Calculating fft spectrum of a signal and obtain correct amplitudes

Question

I'm trying something rather simple which is getting the correct values in a spectrum of a signal with a fixed length. I don't know the containing frequencies and therefor I can't chose the number of fft nodes accordingly. This leaves me with the problem of leakage. I try to avoid that by using a window to fade out the beginning and end of my time signal. With that I'm having problems to get the right amplitudes. Here is what I've done so far:

• compensate the windowing by multiplying the spectrum by a factor
• two-sided to single-sided spectrum -> multiplying spectrum amplitude by two
• normalice spectrum amplitude by dividing by the number of frequency nodes 'N'

As long as the number of fft nodes equals the number in my time signal it works well. But when the number of fft nodes are less than the time signal nodes the amplitudes of the spectrum are getting wrong.

I built a little example which calculates three spectrums:

1. Using fft() with windowing
2. Using fft() without windowing
3. Using spectrogram()

By changing the SampleTime one can change the length of the time signal and by changing N one can change the number of FFT nodes. With the following setting all amplitudes are correct, but for the spectrum with the Hanning window. Why is it that much off? Thank you for any help in advance.

SampleTime=0.001;
Fs=1/SampleTime;
t=0:SampleTime:10;

y=2*sin(2*pi*10*t)+3*cos(2*pi*30*t); % Time Signal

fn = Fs/2; % Nyquistfrequency
N = 1001; % FFT-length (use N=2^x to use DFT-Algorithm)
df = Fs/N; % frequency resolution

%% Generating spectrum using fft()
win = hann(length(y))';  % generate Hanning window for Time Signal

ywin = y.*win*length(win)/sum(win); % amplitude correction to get correct values with window

Hwin = fft(ywin, N); % spectrum of windowed signal
H = fft(y,N);  % complex spectrum of signal without window

amplHwin = abs(Hwin); % absolute values of spectrum of windowed signal
amplH = abs(H); % absolute values of spectrum without window

% Double sided spectrum -> single sided spectrum of windowed signal
amplitudengang = fftshift(amplH/N);
amplitudengang = amplitudengang(ceil(length(amplitudengang)/2+1):end);
amplitudengang(2:end) = amplitudengang(2:end)*2;

% Double sided spectrum -> single sided spectrum of unwindowed signal
amplitudengangwin = fftshift(amplHwin/N);
amplitudengangwin = amplitudengangwin(ceil(length(amplitudengangwin)/2+1):end);
amplitudengangwin(2:end) = amplitudengangwin(2:end)*2;

%% Generating spectrum using spectrogram()
noverlap = floor(0.5*N);
window = hann(N);
[B,F,T] = spectrogram(y,window,noverlap,N,Fs);
% Scaling
FFTScalingFactor = 2 * 1 / N; % 2 for Hanning window, 2^(-1/2) for RMS, 1 for peak
A = FFTScalingFactor * B;
B = FFTScalingFactor * abs(B);
B(2:end-1,:) = 2.*B(2:end-1,:); % single side spectrum with even block length

%% Plots

x_fn = 0 : df : fn-df; % Frequency vector

figure
subplot(2,2,1)
plot(x_fn,amplitudengangwin,'r')
hold on
plot(x_fn,amplitudengang)
legend('Hanning window','no window')
title('FFT()')
subplot(2,2,3)
plot(t,y)
title('Time signal')
subplot(2,2,2)
plot(F,B(:,1))
title('Spectrogram()')

Answer 1

You should change the number of FFT.

Hwin = fft(ywin, N); % spectrum of windowed signal
H = fft(y,N);  % complex spectrum of signal without window

Here N is 1001.

It means when you performed FFT of windowed signal, you used just first 1001 samples in all data samples(10001 samples).

The windowed signal you used would be like below.

Change your code like below

Hwin = fft(ywin, length(y)); % spectrum of windowed signal
H = fft(y,length(y));  % complex spectrum of signal without window
    amplHwin = abs(Hwin); % absolute values of spectrum of windowed signal
amplH = abs(H); % absolute values of spectrum without window

% Double sided spectrum -> single sided spectrum of windowed signal
amplitudengang = fftshift(amplH/length(y));
amplitudengang = amplitudengang(ceil(length(amplitudengang)/2+1):end);
amplitudengang(2:end) = amplitudengang(2:end)*2;

% Double sided spectrum -> single sided spectrum of unwindowed signal
amplitudengangwin = fftshift(amplHwin/length(y));
amplitudengangwin = amplitudengangwin(ceil(length(amplitudengangwin)/2+1):end);
amplitudengangwin(2:end) = amplitudengangwin(2:end)*2;

fftAll=[amplitudengang' amplitudengangwin'];

figure(1)
    plot(fftAll)
    xlim([0 500])

You will get below spectrum.

You defined N=1001, but your data length is 10001. If you want to use just 1001 samples, you should make your window length as 1001.

If you just want use 1001 samples in your data, please change your code like below.

SampleTime=0.001;
Fs=1/SampleTime;
t=0:SampleTime:10;

y=2*sin(2*pi*10*t)+3*cos(2*pi*30*t); % Time Signal

fn = Fs/2; % Nyquistfrequency
N = 1001; % FFT-length (use N=2^x to use DFT-Algorithm)
df = Fs/N; % frequency resolution

%% Generating spectrum using fft()
win = hann(N)';  % generate Hanning window for Time Signal

figure(10)
plot(win)

ywin = y(1:N).*win*length(win)/sum(win); % amplitude correction to get correct values with window
figure(11)
plot(ywin)
xlim([0 1001])
   Hwin = fft(ywin, N); % spectrum of windowed signal
H = fft(y,N);  % complex spectrum of signal without window
    amplHwin = abs(Hwin); % absolute values of spectrum of windowed signal
amplH = abs(H); % absolute values of spectrum without window

% Double sided spectrum -> single sided spectrum of windowed signal
amplitudengang = fftshift(amplH/N);
amplitudengang = amplitudengang(ceil(length(amplitudengang)/2+1):end);
amplitudengang(2:end) = amplitudengang(2:end)*2;

% Double sided spectrum -> single sided spectrum of unwindowed signal
amplitudengangwin = fftshift(amplHwin/N);
amplitudengangwin = amplitudengangwin(ceil(length(amplitudengangwin)/2+1):end);
amplitudengangwin(2:end) = amplitudengangwin(2:end)*2;

fftAll=[amplitudengang' amplitudengangwin'];

figure(1)
    plot(fftAll)
    xlim([0 50])

In this code, I changed all length(y) as N. And I changed ywin like below.

ywin = y(1:N).*win*length(win)/sum(win);

The result will be like below.