IFFT of Gaussian-filtered real signal is complex - why?

15 views (last 30 days)
Hello everyone,
I've written a script to apply a Gaussian filter to a time series (signal), however I must have done something wrong, as the ifft of the filtered signal is complex, even though the original signal was real.
I know that when the first and last elements of the imaginary part vector of the FFT are equal to zero, the symmetry condition is met, i.e. the result of the IFFT will be real rather than complex. However I'm not sure if I should just manually make them be equal to zero.
Could the problem lie in how I've defined the Gaussian filter, i.e. along the entire frequency range of the FFT, and symmetrical?
Code is below. Many thanks in advance for any suggestions.
Y = fft(y);
% filter parameters
f0 = 2; % Hz
sigma = 0.2 * f0;
% define filter
f = 1:length(Y);
gaussFilter_firstHalf = normpdf (f, f0, sigma);
gaussFilter_firstHalf = gaussFilter_firstHalf ./ max(gaussFilter_firstHalf); % scale it so it goes from 0-1
for i=1:length(Y)
gaussFilter_secondHalf(i) = gaussFilter_firstHalf(length(Y)-i+1);
end
gaussFilter = gaussFilter_firstHalf + gaussFilter_secondHalf;
% apply filter
Y_filtered = Y .* gaussFilter;
y_filtered = ifft (Y_filtered);

Accepted Answer

Matt J
Matt J on 2 Dec 2012
Edited: Matt J on 3 Dec 2012
You need to be more careful about using fftshift in appropriate places. You also need to be more careful about how you set up the frequency axis in Hz. It won't be f=1:length(Y). Remember that with DFTs, your frequency sampling interval is related to your time sampling interval as follows
frequencysampling = 1/N/timesampling
leading to
N=length(Y);
f= ( (0:N-1) -ceil((N-1)/2) )/N/timesampling;
gaussFilter = normpdf(abs(f),2, 0.4);
plot(f,gaussFilter); %Check
Y_filtered = Y .* ifftshift(gaussFilter);
y_filtered = ifft(Y_filtered,'symmetric');
plot(f, fftshift(abs(Y_filtered))); %Check
  4 Comments
Matt J
Matt J on 4 Dec 2012
See the documentation for the difference between IFFT and IFFTSHIFT.

Sign in to comment.

More Answers (3)

Image Analyst
Image Analyst on 2 Dec 2012
This is not true: "when the first and last elements of the imaginary part vector of the FFT are equal to zero, the symmetry condition is met". That does not guarantee that the spectrum is Hermitian. It would have to be symmetric about the middle. You might want to check out the chart about 1/3 of the way down this page: http://www.cv.nrao.edu/course/astr534/FourierTransforms.html
If you start with a real time domain signal, you will end up with a Hermitian spectrum (= an even real part and an odd imaginary part). If you then filter it with a symmetric Gaussian, the real part should stay even and the imaginary part should stay odd. So transforming them back should get you a mostly real function with only a very small imaginary part due to spatial quantization.
I can't run your code because I don't have whatever toolbox contains normpdf.
  1 Comment
AwedBy Matlab
AwedBy Matlab on 3 Dec 2012
Thanks!
Normpdf is a function in the Stats Toolbox, which can be replaced by the formula of a Gaussian distribution:
gaussFilter_firstHalf = normpdf (f, f0, sigma);
gaussFilter_firstHalf = (1/sigma*sqrt(2*pi)) * exp (-(f-fo)^2/(2*sigma^2)); %equivalent

Sign in to comment.


Vieniava
Vieniava on 2 Dec 2012
Your Gaussian filter is NOT symmetric in frequency domain. Check this out:
stem(fftshift(gaussFilter))
it supposed to be symmetric with respect to the (length(Y)/2 + 1) sample.

AwedBy Matlab
AwedBy Matlab on 2 Dec 2012
Thanks for both your comments. Unfortuantely, even after making those 2 changes to the code, I still get a Y_filtered which, when plotted, doesn\t look as I expected it to, i.e. the same as the Y spectrum but with only the components within the bell curve having survived

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by