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

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);

Answer 1

[Matt J]

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

[AwedBy Matlab]

Thanks Matt. So it's incorrect to create the Gaussian filter straight in the frequency domain, from two mirror halves concatenated together?

And shouldn't the plot be off abs(fftshift(Y_filtered)), rather than fftshift(abs(Y_filtered))?

Also, could you please remind me by virtue of which property of the DTFT is the frequencysampling equal to 1/(N*timesampling)? I'm getting really confused here..

The plots that I get out of the code below don't look quite right.

timesampling = 1:length(y);
Y = fft(y);
N=length(Y);
f=  (  (0:N-1) -ceil((N-1)/2)  )/N/timesampling;
figure
plot(f, fftshift(abs(Y)));

% filter parameters
f0 = 2; % Hz
sigma = 0.2 * f0; % 20% of f0

% define filter
gaussFilter = normpdf(abs(f),2, 0.4);
figure
plot(f,gaussFilter); %Check

% apply filter
Y_filtered = Y .* ifftshift(gaussFilter);
y_filtered = ifft(Y_filtered,'symmetric');
figure
plot(f, fftshift(abs(Y_filtered))); %Check

[Matt J]

Thanks Matt. So it's incorrect to create the Gaussian filter straight in the frequency domain, from two mirror halves concatenated together?

I don't know that it's incorrect, just more awkward and error-prone as compared to the 1-line construction that I gave you.

And shouldn't the plot be off abs(fftshift(Y_filtered)), rather than fftshift(abs(Y_filtered))?

The two are equivalent.

Also, could you please remind me by virtue of which property of the DTFT is the frequencysampling equal to 1/(N*timesampling)? I'm getting really confused here..

It might be time to dust off your DSP course textbook :-) Basically, though, if you start with the definition of the continuous Fourier transform and you discretize that integral by sampling time at intervals of T and frequency at intervals of 1/(N*T), you will see that it turns into the formula for the DFT.

The plots that I get out of the code below don't look quite right.

We have no way of running your code or knowing how you want the plots to look. However, one error that I see is that you're setting timesampling=1:length(y). timesampling is supposed to be a scalar equal to the physical time in seconds between your y(i) samples.

[AwedBy Matlab]

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]

See the documentation for the difference between IFFT and IFFTSHIFT.

Answer 2

[Image Analyst]

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.

[AwedBy Matlab]

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

Answer 3

[Vieniava]

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.

Answer 4

[AwedBy Matlab]

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