Help needed about Matlab's windowed fft precision!!

EDAboard.com Forum Index -> Digital Signal Processing -> Help needed about Matlab's windowed fft precision!!

Author

Message

qslazio

Joined: 23 May 2004
Posts: 194
Helped: 9

15 Jul 2007 10:18 fft precision

Hi everybody,

I've met a weird problem really makes me desperate.

The problem is as below,

As we know, hann window's spectrum is defined as[1]:
W(phi)=0.5*D(phi)+0.25(D(phi-2*pi/N)+D(phi+2*pi/N))
where, D(phi)=exp(j*phi)*sin(N/2*phi)/sin(1/2*phi) is the Dirichlet kernel, and N is number of DFT points
the amplitude of spectrum at the first fft bin (phi=0) should be N/2
the amplitude of spectrum at the second and last fft bin (phi=1*2pi/N and phi=(N-1)/N*2pi) should be N/4
the amplitude of spectrum at the rest of fft bins should be 0
hence, in theory, there should be 3 non-zero bins on the spectrum of hann window.

But when I do this in Matlab, I found that the rest of fft bins' amplitude is not so small that it cannot be simply neglected.
As can be seen from the picture in the attachment

for N=1024, the first and second bin's amplitude is 0dB and -6dB as expected, but the rest bins' amplitude are not

zero as theory says (they are around -70dB~-80dB).
The consequence of this problem is that when I do fft with any sine wave plus noise of which the noise floor is below these non-signal bins' amplitude (for example the noise floor is around -100dB), Matlab simply gives me the wrong SNR!!

I can not find any answers or related information from any textbook. I guess it is caused by Matlab fft's internal precision.

Anyone has idea about this?

Thanks a lot!

[1] "On the use of windows for harmonic analysis with the discrete fourier transform", F. J. Harris, Proceedings of the IEEE VOL.66 NO.1, Jan. 1978

the matlab script I used (I'm using Matlab R2006b):

close all;
nbpts=1024;
w=hann(nbpts); % generates hanning window vectors
W=fft(w)/sum(w); % generates normalized fft spectrum
WM=abs(W); % gets amplitude spectrum
figure(1);
plot([1:nbpts],20*log10(WM),'bo-');
xlabel('fft bin');
ylabel('dB');
title('Hann window DFT with N=1024 normalized spectrum');
grid on;

Sorry, but you need login in to view this attachment

lena

Joined: 01 Jul 2006
Posts: 18

31 Jul 2007 8:09 fft matlab hanning window

in matlab follwing function used for

hanning window

w(n) = 0.5(1-cos(2*pi*n/N)) 0<=n<=N-1

bulx

Joined: 07 Aug 2004
Posts: 171
Helped: 24

01 Aug 2007 14:47 matlab fft hanning window example

Strange that you are looking at the frequency spectrum of the hannng window vectors itself. perhpas you already know this: The hanning window is to compensate for distortion caused by the limited number of samples on which FFT is done in a particular application. ie; if you had the oppurtunity to take 1 million sets, each of 1024 samples of the of your signal, in sequence, then take 1 million FFTs and average them, a hanning widow is unnecessary.

To look at how good your windowing is, sample a pure sine (or some signal with limited and known frequncy components) and see how the noise floor is. There is of course some finite precision effect due to of matlab, but it is should be more or less ignorable as compared to what an application implemented with hardware or a dsp processor would produce. In fact, most applications are more concerned about how close they can get to matlab.
- b

ILv2Xlr8

Joined: 04 Mar 2009
Posts: 12

05 Mar 2009 22:02 energy correction hann window

What is the coherent gain factor for a hanning window?
could this cause error in calculating the magnitudes of the FFT bins?

ILv2Xlr8

Joined: 04 Mar 2009
Posts: 12

06 Mar 2009 16:52 correction spectrum for hann window

I have determined that my window and CGF is not causing the error. I am getting errror when I integrate the data either in the time domain or frequency domain. I am discussing this in another post.

http://www.edaboard.com/viewtopic.php?t=342852

rramya

Joined: 14 Dec 2008
Posts: 76
Helped: 19

06 Mar 2009 18:21 fft precision

Hi ILv2Xlr8,
ILv2Xlr8 wrote:
What is the coherent gain factor for a hanning window?
ans:
coherent gain factor for a hanning window is 2

Multiply by 2 for amplitude, sqrt(8/3) for energy.

If you integrate the hanning window, you'll get the amplitude correction.
If you integrate the square of the window, you'll get the energy correction.

Which correction you use is your choice and depends on your application.

if possible, i will tell u the root cause of error due to integration.

happy learning

ILv2Xlr8

Joined: 04 Mar 2009
Posts: 12

06 Mar 2009 18:34 x =0:1:50; in matlab??

Thanks for the clarification.

I am getting this magnitude error only when integrating the windowed data in the time domain, or integrating the windowed FFT data in the frequency domain.

rramya

Joined: 14 Dec 2008
Posts: 76
Helped: 19

06 Mar 2009 20:44 precision on figures matlab

Hi ILv2Xlr8,
can i ask you a silly question ( suddenly mind is very blank)
Y do u want to
integrate the windowed data in the time domain, or integrate the windowed FFT data in the frequency domain.

I cant get U .

I will send the spectrum analysis usinf fft. It contains the code as well.

If possible post your code.

Happy learning.

Added after 1 minutes:

I forgot to attach the file

Google
AdSense
Google Adsense

06 Mar 2009 20:44 Ads

Sorry, but you need login in to view this attachment

aarr

Joined: 06 Mar 2009
Posts: 5

06 Mar 2009 20:47 hanning window fft

thanks
i want some simple examples using MATLAB functions to desigen digital filters or any simple projects by MATLAB to see results to understand DSP.
thanks .

Aya2002

Joined: 12 Dec 2006
Posts: 1409
Helped: 254
Location: Iraq

07 Mar 2009 0:51 matlab fft bode

Hi aarr,

see this example: to estimate the PSD using Welch's method (see matlab help).

t=0:0.001::1;
y=sin(2*pi*10*t);
subplot(211);
plot(t,y),grid;
subplot(212);
pwelch(y);
------------------------------------------------------------------
example:

Fs=40e3; % sampling frequency
W=[8e3, 12e3]/(Fs/2); % passband spec
[b,a]=butter(4,W,'z'); % design Butterworth
% note the order is 2*4=8
[d,c]=cheby1(4,1,W,'z'); % design Chebyshev
figure(1) % make Bode plots
freqz(b,a,60,Fs) % Butterworth in 60 points
hold on % both in the same diagram
freqz(d,c,60,Fs) % Chebyshev
figure(2) % make pole/zero plot
pzmap(b,a) % for Butterworth
figure(3) % make pole/zero plot
pzmap(d,c) % for Chebyshev

----------------------------------------------------------------

example:

% Generate the sine wave sequence
fs = 8000; T=1/fs; % Sampling rate and sampling period
x = 2*sin(2000*pi*[0:1:50]*T); %Generate the 51 2000-Hz samples
% Apply the FFT algorithm
N = length(x);
index_t=[0:1:N-1];
f = [0:1:50]*8000/N; %Map the frequency bin to the frequency (Hz)
length(f)
xf=abs(fft(x))/N; %Calculate the amplitude spectrum
length(xf)
figure(1)
%Using the Bartlett window
x_b=x.*bartlett(N)'; %Apply the triangular window function
xf_b=abs(fft(x_b))/N; %Calculate the amplitude spectrum
subplot(2,2,1);plot(index_t,x);grid
xlabel('Time index n'); ylabel('x(n)');
subplot(2,2,3); plot(index_t,x_b);grid
xlabel('Time index n'); ylabel('Triangular windowed x(n)');
subplot(2,2,2);plot(f,xf);grid;axis([0 8000 0 1]);
xlabel('Frequency (Hz)'); ylabel('Ak (no window)');
subplot(2,2,4); plot(f,xf_b);grid; axis([0 8000 0 1]);
xlabel('Frequency (Hz)'); ylabel('Triangular windowed Ak');

enjoy

rramya

Joined: 14 Dec 2008
Posts: 76
Helped: 19

07 Mar 2009 6:53 fft coherent gain bartlett

Hi ILv2Xlr8,
From math works.com site i got the following code: ( for spectrum analysis using hanning window)

function [vFrequency, vAmplitude] = fastfft(vData, SampleRate, Plot)

%FASTFFT Create useful data from an FFT operation.
% Usage: [vFrequency, vAmplitude] = fastfft(vData, SampleRate, [Plot])
%
% (no plot will be shown if the last input == 0 or is not included)
%
% This function inputs 'vData' as a vector (row or column),
% 'SampleRate' as a number (samples/sec), 'Plot' as anything,
% and does the following:
%
% 1: Removes the DC offset of the data
% 2: Puts the data through a hanning window
% 3: Calculates the Fast Fourier Transform (FFT)
% 4: Calculates the amplitude from the FFT
% 5: Calculates the frequency scale
% 6: Optionally creates a Bode plot
%
% Created 7/22/03, Rick Auch, mekaneck(at)campbellsville.com

%Make vData a row vector
if size(vData,2)==1
vData = vData';
end

%Calculate number of data points in data
n = length(vData);

%Remove DC Offset
vData = vData - mean(vData);

%Put data through hanning window using hanning subfunction
vData = hanning(vData);

%Calculate FFT
vData = fft(vData);

%Calculate amplitude from FFT (multply by sqrt(8/3) because of effects of hanning window)
vAmplitude = abs(vData)*sqrt(8/3);

%Calculate frequency scale
vFrequency = linspace(0,n-1,n)*(SampleRate/n);

%Limit both output vectors due to Nyquist criterion
DataLimit = ceil(n/2);
vAmplitude = vAmplitude(1:DataLimit);
vFrequency = vFrequency(1:DataLimit);

if exist('Plot', 'var')==1 & Plot~=0
plot(vFrequency, vAmplitude);
title('Bode Plot');
xlabel('Frequency (Hz)');
ylabel('Amplitude');
end

%------------------------------------------------------------------------------------------
%Hanning Subfunction
function vOutput = hanning(vInput)
% This function takes a vector input and outputs the same vector,
% multiplied by the hanning window function

%Determine the number of input data points
n = length(vInput);

%Initialize the vector
vHanningFunc = linspace(0,n-1,n);

%Calculate the hanning funtion
vHanningFunc = .5*(1-cos(2*pi*vHanningFunc/(n-1)));

%Output the result
vOutput = vInput.*vHanningFunc;

Added after 3 minutes:

If you further examples: just visit this site. it will be very much useful to u to proceed further.
http://blinkdagger.com/matlab/matlab-windowing-part-2
http://blinkdagger.com/matlab/matlab-windowing-part-1
http://blinkdagger.com/matlab/matlab-windowing-part-3

Happy learning