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 件のコメント
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 ExchangeFilter Analysis についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by