{"id":381,"date":"2011-07-22T14:12:21","date_gmt":"2011-07-22T18:12:21","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2011\/07\/22\/filtering-fun\/"},"modified":"2019-10-29T16:49:09","modified_gmt":"2019-10-29T20:49:09","slug":"filtering-fun","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2011\/07\/22\/filtering-fun\/","title":{"rendered":"Filtering fun"},"content":{"rendered":"
\r\n

Today I want to take the test pattern I created last time<\/a> and subject it to a variety of frequency-based filters.\r\n <\/p>\r\n

In this post I'll be using a variety of frequency design and frequency-response visualization tools in MATLAB, Signal Processing\r\n Toolbox, and Image Processing Toolbox, including fir1<\/tt><\/a>, ftrans2<\/tt><\/a>, freqspace<\/tt><\/a>, fwind1<\/tt><\/a>, and freqz2<\/tt><\/a>. I won't be giving much explanation of these filter design techniques. Perhaps in the future I'll explore filter design here\r\n in more detail. For now, though, I really just want to have some fun with this test pattern.\r\n <\/p>

[x,y] = meshgrid(-200:200);\r\nr = hypot(x,y);\r\nkm = 0.8*pi;\r\nrm = 200;\r\nw = rm\/10;\r\nterm1 = sin( (km * r.^2) \/ (2 * rm) );\r\nterm2 = 0.5*tanh((rm - r)\/w) + 0.5;\r\ng = (term1 .* term2 + 1) \/ 2;\r\n\r\nimshow(g, 'XData'<\/span>, [-.8 .8], 'YData'<\/span>, [-.8 .8])\r\naxis on<\/span>\r\nxlabel('Normalized frequency'<\/span>), ylabel('Normalized frequency'<\/span>)\r\ntitle('Test pattern'<\/span>)<\/pre> 
Example 1:<\/b> Design a one-dimensional lowpass filter and apply it in both directions.\r\n   <\/p>
b = fir1(30,0.3);\r\nh = b' * b;\r\ng1 = imfilter(g, h);\r\n\r\nimshow(g1, 'XData'<\/span>, [-.8 .8], 'YData'<\/span>, [-.8 .8])\r\naxis on<\/span>\r\nxlabel('Normalized frequency'<\/span>), ylabel('Normalized frequency'<\/span>)<\/pre> 
Example 2:<\/b> Design a one-dimensional highpass filter and apply it in both directions.\r\n   <\/p>
b2 = fir1(30, 0.3, 'high'<\/span>);\r\nh2 = b2' * b2;\r\ng2 = imfilter(g, h2);\r\n\r\nimshow(g2, [], 'XData'<\/span>, [-.8 .8], 'YData'<\/span>, [-.8 .8])\r\naxis on<\/span>\r\nxlabel('Normalized frequency'<\/span>), ylabel('Normalized frequency'<\/span>)<\/pre> 
Example 3:<\/b> Design a one-dimensional bandpass filter and apply it on both directions.\r\n   <\/p>
b3 = fir1(30, [0.3 0.4]);\r\nh3 = b3' * b3;\r\ng3 = imfilter(g, h3);\r\n\r\nimshow(g3, [], 'XData'<\/span>, [-.8 .8], 'YData'<\/span>, [-.8 .8])\r\naxis on<\/span>\r\nxlabel('Normalized frequency'<\/span>), ylabel('Normalized frequency'<\/span>)<\/pre> 
Example 4:<\/b> Apply the one-dimensional bandpass filter from the previous step in only one direction.\r\n   <\/p>
g4 = imfilter(g, b3);\r\n\r\nimshow(g4, [], 'XData'<\/span>, [-.8 .8], 'YData'<\/span>, [-.8 .8])\r\naxis on<\/span>\r\nxlabel('Normalized frequency'<\/span>), ylabel('Normalized frequency'<\/span>)<\/pre> 
Example 5:<\/b> Use ftrans2<\/tt> to turn the one-dimensional bandpass filter from the previous example into a two-dimensional bandpass filter that is approximately\r\n      circularly symmetric.\r\n   <\/p>
h5 = ftrans2(b3);\r\ng5 = imfilter(g, h5);\r\n\r\nimshow(g5, [], 'XData'<\/span>, [-.8 .8], 'YData'<\/span>, [-.8 .8])\r\naxis on<\/span>\r\nxlabel('Normalized frequency'<\/span>), ylabel('Normalized frequency'<\/span>)<\/pre> 
Example 6:<\/b> Use fwind1<\/tt> to design a two-dimensional filter that has a wedge-shaped frequency response.\r\n   <\/p>
[f1, f2] = freqspace(101, 'meshgrid'<\/span>);\r\ntheta = cart2pol(f1, f2);\r\nHd = ((theta >= pi\/6) & (theta < pi\/3)) | ...<\/span>\r\n    ((theta >= -5*pi\/6) & (theta < -2*pi\/3));\r\nmesh(f1, f2, double(Hd))\r\ntitle('Ideal frequency response'<\/span>)<\/pre> 
h6 = fwind1(Hd, hamming(31));\r\nfreqz2(h6)\r\ntitle('Designed filter frequency response'<\/span>)<\/pre> 
g6 = imfilter(g, h6);\r\n\r\nimshow(g6, [], 'XData'<\/span>, [-.8 .8], 'YData'<\/span>, [-.8 .8])\r\naxis on<\/span>\r\nxlabel('Normalized frequency'<\/span>), ylabel('Normalized frequency'<\/span>)<\/pre> 
Example 7:<\/b> Compare the use of an averaging filter and a Gaussian filter. This example is inspired by Cris Luengo's comment<\/a> from last week.\r\n   <\/p>
g7 = imfilter(g, fspecial('average'<\/span>, 9));\r\nimshow(g7, [], 'XData'<\/span>, [-.8 .8], 'YData'<\/span>, [-.8 .8])\r\naxis on<\/span>\r\nxlabel('Normalized frequency'<\/span>), ylabel('Normalized frequency'<\/span>)\r\ntitle('13-by-13 averaging filter'<\/span>)<\/pre> 
You can see that some very high frequencies remain in the output. The Gaussian filter performs much better:<\/p>
g9 = imfilter(g, fspecial('gaussian'<\/span>, 9, 2));\r\nimshow(g9, [], 'XData'<\/span>, [-.8 .8], 'YData'<\/span>, [-.8 .8])\r\naxis on<\/span>\r\nxlabel('Normalized frequency'<\/span>), ylabel('Normalized frequency'<\/span>)\r\ntitle('13-by-13 Gaussian filter'<\/span>)<\/pre> 
OK, I'm going to give this test pattern a rest for now. Phew!<\/p>