CODEWITHSUNDEEP
matlabfft
Matthias Pospiech
Matthias Pospiech

Reputation: 3494

implement low pass filter in matlab

I want to implement a filter that eliminates higher frequencies. In this example I want to eliminate the sin curve and keep the linear curve.

EDIT:

I corrected my code, however the function that I implemented for filtering changes the edges of the data significantly which is not acceptable.

clc; clear all;
xaxis = linspace(1, 10, 1000);
data = xaxis + sin(xaxis*3);
Nf = 2^12;
xAxisf = linspace(-10,10,Nf);
% plot(xaxis, data);

% FFT

xsize = numel(data);    
Xf = zeros([1 Nf]);
indices = Nf/2-floor(xsize/2):Nf/2-floor(xsize/2)+xsize - 1;
Xf(indices) = data;

% Xf = fftshift(Xf);
Xf = fft(Xf);
Xf = fftshift(Xf);

% plot 
Xfa = abs(Xf); plot(xAxisf, Xfa);

% generate super-gaussian filter function
Nf = numel(Xf);    
widthfilter = 0.12;
filterpower = 2;
filter = exp(-(xAxisf.^2./widthfilter^2).^filterpower);

% filter
filtertimes = 20;
Xf = Xf .* filter.^filtertimes;

% plot 
Xfa = abs(Xf); plot(Xfa);

% iFFt
Xfs = ifftshift(Xf);
Xif = ifft(Xfs);
% Xif = ifftshift(Xif);
result = abs(Xif);

plot(result(indices))

Upvotes: 0

Views: 13284

Answers (1)

thang
thang

Reputation: 3466

First issue:

    Xf = fftshift(data);     % NOT NEEDED
    Xf = fft(Xf);
    Xf = fftshift(Xf);

Do not fftshift the data before fft. The shift is only needed AFTER fft. This is because the radix-n(probably 2) fft "decimates" the data in the process. You don't need to fix it before because it hasn't been decimated.

Second issue: 

    Xfs = ifftshift(Xf);
    Xif = ifft2(Xfs);            
    Xif = ifftshift(Xif);   % NOT NEEDED

ifftshift re-decimates the data (undoes fftshift), which ifft requires as input. The output of ifft reconstructs the original signal if the input is already decimated. DO NOT ifftshift after.

Third issue: 

    Xfs = ifftshift(Xf);
    Xif = ifft2(Xfs);       % USE ifft INSTEAD OF ifft2     
    Xif = ifftshift(Xiff);

Why in the world did you switch to 2D ifft all of a sudden?

I didn't look at your filter code in detail, but I would like to remark that if you want a low pass filter, it needs to be symmetric around the mid point. Otherwise your frequency response is not symmetric and you're gonna end up with at bunch of imaginaries.

And please change your title. This isn't a "Fourier filter". It is a low pass filter using the window method and fft. Window in that you're applying a window in frequency space.

Ok, it's late and I am getting cranky from the back and forth... just trying to help. Faster for me to just write the code for you.

If you're looking for an effect of the filter in your code, you're not gonna be able to because the cut off frequency of your filter is too high and/or the frequency of the sine wave in your data is too low. Here is a version where I increased the frequency of oscillation of the input sine wave:

clc; clear all;
xaxis = linspace(1, 10, 1000);
data = xaxis + sin(xaxis*10);
% plot(xaxis, data);

% FFT
Xf = data;
Xf = fft(Xf);
Xf = fftshift(Xf);

% generate super-gaussian filter function
Nf = numel(data);
xAxisf = linspace(-5,5,Nf);
widthfilter = 0.1;
filterpower = 2;
filter = exp(-(xAxisf.^2./widthfilter^2).^filterpower);

% filter
filtertimes = 1;
Xf = Xf .* filter.^filtertimes;

% plot
Xfa = abs(Xf); plot(Xfa);

% iFFt
Xfs = ifftshift(Xf);
Xif = ifft(Xfs);
result = abs(Xif);

plot(result); hold on; plot(data,'r');
    legend('filtered','data');

going to bed. good night! did my public service :p

Upvotes: 4

Related Questions