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

8 次查看（过去 30 天）
AwedBy Matlab
AwedBy Matlab 2012-12-1
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);

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2012-12-2
编辑：Matt J 2012-12-3
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 个评论
显示 2更早的评论
AwedBy Matlab
AwedBy Matlab 2012-12-4
编辑：AwedBy Matlab 2012-12-5
Thanks, the filtering script is now working as it should.
Just one more question. How come the FFT of the filtered signal,
Y_filtered = Y .* ifftshift(gaussFilter);
..is obtained by multiplying the unfiltered signal by the IFFT of the filter shape and not with the filter itself as it was defined, i.e. in the frequency domain? I would have thought that this needs to be the product of two frequency-domain signals, i.e. Y.*gaussFilter , whereas ifftshift(gaussFilter) must surely be a time-domain signal!
Matt J
Matt J 2012-12-4
See the documentation for the difference between IFFT and IFFTSHIFT.

请先登录，再进行评论。

更多回答（3 个）

Image Analyst
Image Analyst 2012-12-2
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 个评论
AwedBy Matlab
AwedBy Matlab 2012-12-3
编辑：AwedBy Matlab 2012-12-3
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

请先登录，再进行评论。

Vieniava
Vieniava 2012-12-2
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 2012-12-2
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

类别

在 Help Center 和 File Exchange 中查找有关 Filter Analysis 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by