Why the amplitude of theFFT trasform is not equal as the amplitude of the input signal?
5 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Arianna Senesi
el 26 de Abr. de 2021
Respondida: Mathieu NOE
el 27 de Abr. de 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 comentarios
Respuesta aceptada
Mathieu NOE
el 27 de Abr. de 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
0 comentarios
Más respuestas (1)
Mathieu NOE
el 26 de Abr. de 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 comentarios
Mathieu NOE
el 27 de Abr. de 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
Ver también
Categorías
Más información sobre Spectral Measurements en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!