MATLAB Answers

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

5 views (last 30 days)
Arianna Senesi
Arianna Senesi on 26 Apr 2021
Answered: Mathieu NOE on 27 Apr 2021
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))])
  4 Comments

Sign in to comment.

Accepted Answer

Mathieu NOE
Mathieu NOE on 27 Apr 2021
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

More Answers (1)

Mathieu NOE
Mathieu NOE on 26 Apr 2021
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
  2 Comments
Mathieu NOE
Mathieu NOE on 27 Apr 2021
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

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!