Why the amplitude of theFFT trasform is not equal as the amplitude of the input signal?

close all
clear all
clc
sine = [];
rate = 100000;
fbase=1300;
freq = fbase+100; 
sinesine1 = 0.25*sin(linspace(0, 2*pi,rate/freq)); 
freq = fbase+200;
sinesine2 = 0.25*sin(linspace(0,2*pi,rate/freq)); 
freq = fbase+300; 
sinesine3 =0.25*sin(linspace(0, 2*pi,rate/freq)); 
freq = fbase+400;
sinesine4 = 0.25*sin(linspace(0, 2*pi,rate/freq)); 
freq =fbase+500; 
sinesine5 = 0.25*sin(linspace(0, 2*pi,rate/freq));
n=200; 
for i = 1: n
    sine = [sine sinesine1(1:end-1)];
end
for i = 1: n
    sine = [sine sinesine2(1:end-1)];
end
for i = 1: n
    sine = [sine sinesine3(1:end-1)];
end
for i = 1: n
    sine = [sine sinesine4(1:end-1)];
end
for i = 1: n
    sine = [sine sinesine5(1:end-1)];
end
Fs = rate;
L = length(sine);     % Length of signal
T = 1/Fs;             % Sampling period 
t = (0:L-1)*T;        % Time vector
Z = fft(sine);
P2= abs(Z)./L;
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
figure('Name','Spettro 1400-1800','NumberTitle','off');
plot(f,P1)
axis([fbase fbase+1100 0 max(P5(2:end))+0.01*max(P5(2:end))])

[Walter Roberson]

You use L before you define it.

You do not define L2 but you use it.

[Arianna Senesi]

Now it's ok but it doesn't work correctly, the amplitude it's wrong

[Matt J]

Now it's ok

No, it's not. I still get errors when I run what you've posted.

[Arianna Senesi]

Where does it give you errors? because it works for me. thanks anyway

Answer 1

[Mathieu NOE]

hi again

I was not happy to see that your code needed a very high sampling rate to display correctly the 5 frequencies

this is by no mean required by Shannon's theorem

I modified the way the data are generated and now we can have the right display even with Fs = 5000;

also one single for loop needed vs five

close all
clearvars
clc
sine = [];
% rate = 100000;
% fbase= 1300;
% freq = fbase+100; 
% sinesine1 = 0.25*sin(linspace(0, 2*pi,rate/freq)); 
% freq = fbase+200;
% sinesine2 = 0.25*sin(linspace(0,2*pi,rate/freq)); 
% freq = fbase+300; 
% sinesine3 =0.25*sin(linspace(0, 2*pi,rate/freq)); 
% freq = fbase+400;
% sinesine4 = 0.25*sin(linspace(0, 2*pi,rate/freq)); 
% freq =fbase+500; 
% sinesine5 = 0.25*sin(linspace(0, 2*pi,rate/freq));
% n=200; 
% for i = 1: n
%     sine = [sine sinesine1(1:end-1)];
% end
% for i = 1: n
%     sine = [sine sinesine2(1:end-1)];
% end
% for i = 1: n
%     sine = [sine sinesine3(1:end-1)];
% end
% for i = 1: n
%     sine = [sine sinesine4(1:end-1)];
% end
% for i = 1: n
%     sine = [sine sinesine5(1:end-1)];
% end
% Fs = rate;

% my suggestion
periods = 200;
Fs = 5e3;
dt = 1/Fs;
fbase= 1300;
freq = fbase + [100 200 300 400 500];
for ci = 1:length(freq)
    omega = 2*pi*freq(ci);
    t = (0:dt:periods/freq(ci)-dt);
    sine = [sine 0.25*sin(omega*t)];
end


% T = 1/Fs;             % Sampling period 
% t = (0:L-1)*T;        % Time vector
% L = length(sine);     % Length of signal
% Z = fft(sine);
% P2= abs(Z)./L;
% P1 = P2(1:L/2+1);
% P1(2:end-1) = 2*P1(2:end-1);
% f = Fs*(0:(L/2))/L;
% figure('Name','Spettro 1400-1800','NumberTitle','off');
% plot(f,P1)
% axis([fbase fbase+1100 0 max(P5(2:end))+0.01*max(P5(2:end))])

%%%%%%%%%% PEAK HOLD AVERAGING FFT SPECTRUM DEMO %%%%%%%%%%%%%%%%%
NFFT = 250;    
Overlap = 0.75;

[freq,fft_spectrum] = myfft_peak(sine(:), Fs, NFFT, Overlap); % data must be column oriented (samples x channels)

figure(1),plot(freq,fft_spectrum,'b');grid
title(['Peak Hold FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(' Amplitude')

%%%%%%%%%% SPECTROGRAM DEMO %%%%%%%%%%%%%%%%%

[S,F,T] = myspecgram(sine, Fs, NFFT, Overlap); 

figure(2); 
imagesc(T,F,(abs(S))) 
colorbar('vert');
set(gca,'YDir','Normal') 
xlabel('Time (secs)') 
ylabel('Freq (Hz)') 
title('Short-time Fourier Transform spectrum')
colormap('jet');

function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

if channels > samples
    signal = signal'; % flip dimensions
end

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = max(fft_spectrum,(abs(fft(sw))*4/nfft));     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
%     fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

function  [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal  (example sinus amplitude 1   = 0 dB after fft).
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_specgram = [];
    for ci=1:spectnum
        start = (ci-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_specgram = [fft_specgram abs(fft(sw))*4/nfft];     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;

% time vector 
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(offset/2))/Fs;

end

Answer 2

[Mathieu NOE]

hello

FYI see code below

it iplement a peak hold fft averaging technique; also the spectrogram principle is coded (so you can see what matlab is doing)

both methods show the right amplitude

I reduced your sampling freq that was too high for the demo purpose

close all
clearvars
clc
sine = [];
rate = 10000;
fbase= 1300;
freq = fbase+100; 
sinesine1 = 0.25*sin(linspace(0, 2*pi,rate/freq)); 
freq = fbase+200;
sinesine2 = 0.25*sin(linspace(0,2*pi,rate/freq)); 
freq = fbase+300; 
sinesine3 =0.25*sin(linspace(0, 2*pi,rate/freq)); 
freq = fbase+400;
sinesine4 = 0.25*sin(linspace(0, 2*pi,rate/freq)); 
freq =fbase+500; 
sinesine5 = 0.25*sin(linspace(0, 2*pi,rate/freq));
n=200; 
for i = 1: n
    sine = [sine sinesine1(1:end-1)];
end
for i = 1: n
    sine = [sine sinesine2(1:end-1)];
end
for i = 1: n
    sine = [sine sinesine3(1:end-1)];
end
for i = 1: n
    sine = [sine sinesine4(1:end-1)];
end
for i = 1: n
    sine = [sine sinesine5(1:end-1)];
end
Fs = rate;
% T = 1/Fs;             % Sampling period 
% t = (0:L-1)*T;        % Time vector
% L = length(sine);     % Length of signal
% Z = fft(sine);
% P2= abs(Z)./L;
% P1 = P2(1:L/2+1);
% P1(2:end-1) = 2*P1(2:end-1);
% f = Fs*(0:(L/2))/L;
% figure('Name','Spettro 1400-1800','NumberTitle','off');
% plot(f,P1)
% axis([fbase fbase+1100 0 max(P5(2:end))+0.01*max(P5(2:end))])

%%%%%%%%%% PEAK HOLD AVERAGING FFT SPECTRUM DEMO %%%%%%%%%%%%%%%%%
NFFT = 500;    
Overlap = 0.75;

[freq,fft_spectrum] = myfft_peak(sine(:), Fs, NFFT, Overlap); % data must be column oriented (samples x channels)

figure(1),plot(freq,fft_spectrum,'b');grid
title(['Peak Hold FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(' Amplitude')

%%%%%%%%%% SPECTROGRAM DEMO %%%%%%%%%%%%%%%%%

[S,F,T] = myspecgram(sine, Fs, NFFT, Overlap); 

figure(2); 
imagesc(T,F,(abs(S))) 
colorbar('vert');
set(gca,'YDir','Normal') 
xlabel('Time (secs)') 
ylabel('Freq (Hz)') 
title('Short-time Fourier Transform spectrum')
colormap('jet');

function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

if channels > samples
    signal = signal'; % flip dimensions
end

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = max(fft_spectrum,(abs(fft(sw))*4/nfft));     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
%     fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

function  [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal  (example sinus amplitude 1   = 0 dB after fft).
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_specgram = [];
    for ci=1:spectnum
        start = (ci-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_specgram = [fft_specgram abs(fft(sw))*4/nfft];     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;

% time vector 
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(offset/2))/Fs;

end

[Arianna Senesi]

Thanks for your help. But for what I have to do I need to use a sampling frequency=100k Hz. Moreover in the figure I should see 5 peaks because I go from 1400 to 1800 hz. If I send this code I only get 3 peaks which are at 1600, 2000 and 2500 Hz.

[Mathieu NOE]

Ok

so let's put back the rate at 100 kHz, then the 5 peaks will show up but you'll have to increase the NFFT to 10 times higher if you'd like to distinguish the relatively close frequencies