Wiener Filter + Inverse Filter, Contd.

In the last post i derived the formula for wiener filter. I am not getting enough time to write something here. This night i decided to write something.  If you carefully examine the wiener filter formula it can be seen that when the K is zero ( that is no noise),it act just an inverse filter.
     This means we can just divide the Fourier transform of the input signal(degraded) with the Fourier transform of the degrade function. This will produce the original data. You may be wondering how this works. It is based on convolution theorem.

In my pc mathematica stopped working. So i am switching to matlab from now on

Following mat lab snippet proves that convolution in spatial domain and frequency domain is same

%% in this example convolution signal is not centered ( so the source)
%% by this way we can avoid the one shifting(ifftshift) at the end.

data1 = [1,2,3,2,1,0,0];
kernel = [1,1,1];
kernel_padded = [1,1,1,0,0,0,0];
'convolution spatial'
conv(data1,kernel,'valid')
d = fft(data1) .* fft(kernel_padded);
'fourier'
ifft(d)

When you run it this will produce the result  1  3    6   7  6  3  1. (conv function will out two more zeros, but it is not a part of input data)

Using this idea inverse filter works. Since inverse filter uses a division problem can happen when the denominator is zero or say too high. because of this it is not stable as itself also with noise.

Here is the results and code

a = imread('f:\\ip\\imtests\\exotic.jpg');
d = rgb2gray(a);
degrade = [.1 .1 .1;.1 .5 .1;.1 .1 .1];
dim = conv2(d,degrade,'valid');
figure('Name','degraded image with some filter')
imshow(dim,[0,255])

img_size = size(d) ;
img_size = size(dim); 
conv_mat = zeros(img_size(1),img_size(2));

% create conv matrix for inverse.
pw = floor(img_size(1))-3
ph = floor(img_size(2))-3
conv_mat = padarray(degrade,[pw,ph],'post');

fo = fft2(dim);
fc = fft2(conv_mat);
% the final division.
myimg = fo ./ (fc )

corrected = ifft2(myimg);
figure('Name','corrected with inverse filter')
imshow(real(corrected),[0 255])

Oringal Image

Corrected with Inverse

Don't expect this kind of correction all cases , because of that division things can go totally wrong . So one idea is to add some constant with denominator before division and play with it until you find some satisfactory results Or you need to check it against zero or non-numeric value.

Ok its time for me to stop. I just made this because of me :) . I was worried about not writing anything for sometime. That's why this quick post.  More interesting this are coming(to my mind)

Frequency Identification

In this post i will try to explain the basic things to identify major frequency components from a data.This helpful filtering certain frequencies from the source or for doing some analysis. Frequency analysis a must learn thing if you are learning image processing, signal processing.
 Basic idea is to convert the data to frequency domain using discrete fourier transformation, and using that we can find the major frequency in the data.Before going to the details lets look the equation of a basic Sin wave

it is sin( 2*Pi/T * t )
where T is the time period of wave

            t is the instantaneous  time.
            Pi Ofcourse 22/7 , 
Lets plot this wave. I am using mathematica to plot the wave. 
Here is the wave form
T = 20
Plot[ Sin[2*Pi/T*t], {t, 0, 50}] . (We can verify T =20 from the graph.)

So the frequency of this wave is 1/T . that is 1/20.

Ok.  this is no big deal unless you are very weak in mathematics. Now lets add some random noise to it

n = 1500;

T = 20;

SampledData= Table[Sin[2 Pi* x/T] + RandomReal[.8], {x, n}] ;

ListLinePlot[SampledData]

I hope now you cannot identify the frequency from this ;) .Lets find the T period from this , it must approximate to 20.Before that lest look the equation of one dimensional DFT

Where N is the total number of elements and k/N will give the frequency (because Sin (2*Pi*f * n ) .Ok , now i am going to do DFT with data generated with equation Sin[2 Pi* x/T] and I am  going to find the power spectrum(magnitude of complex numbers)

DftData = Abs[Fourier[SampledData]];

lets plot the power spectrum graph

In graph you can see two spikes. The rightmost spike represent the nyquest frequency( i will explain it later). lets find the first position where magnitude is highest.

pos = Position[f, Max[f]][[1, 1]]

it will be at 76.

Now we can find the frequency by substituting this value for k in equation k/N .

that is 76/1500 = .0506. which approximates to 1/20.

If you want to find the timer period just find the reciprocal. 

T = 1500/76 = 19.73!! approximates to 20.

Hey you just learnt a great thing!!.