FFT and PSD - normalize values

70 vues (au cours des 30 derniers jours)
John Navarro
John Navarro le 6 Mai 2020
Hello.
I have noticed many posts regarding frequency spectrum and many answers to the same questions. How to normalize the values of FFT and PSD from a periodic signal?
Most of us when we refeer to "normalize the frequency domain" we mean that it is required to obtain the frequency domain of the signals where the amplitude of the signal, or the time window, do not affect the obtained plot.
For example in my case I use a periodic sine signal with a frequency of 160HZ, and observed that the amplitude of the foundamental frequency in both figures changes according to the length of the signal. So the questions are:
How to obtain a FFT and PSD where the amplitude values are not affected by:
  1. The time window or the number of points.
  2. The peak amplitude of the signal. I mean for example if the signal amplitude is reduced by the half, still getting the corresponding peak value, or even better the corresponding percentage.
Notice that even the amplitud of the signal is 15. The obtained values in FFT and PSD are not 15.
For FFT the value varies between 10 and 14 depending on the time window, but never the 15 that it should be, even it the time window is very high.
For PSD the amplitude value increases as a function of the time window. From 2 up to 800.
Also notice that PSD shows the foundamental frequency at 160 Hz but in the one-side FFT varies 145 - 160 Hz. How to solve it? Why is this happening?
I attached the code and in the PDF the results with diverse time windows.
Do you have any comments about how to fix my code in order to get the same amplitude for PSD and FFT despite of the time window. Thanks
%% R is the signal data points extracted from the timetable Datos03
%% tiempo is the signal time points.
%% L is the number of data points (Length of the signal)
%Time domain
plot(Datos3.Time,Datos3.VarName1)
R=Datos3.VarName1;
tiempo=seconds(Datos3.Time);
% FFT
Fs=10240; % Sampling frequency
T=1/Fs; % Sampling period
a=size(tiempo) % Number of points
L = a(1,1) % Length of signal
Y = fft(R);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
plot(f,P1)
xlim([0 400])
title('FFT')
%PSD
[SignalSpectrum,SignalFrequencies] = periodogram(R, ...
[],[],Fs);
plot(SignalFrequencies,SignalSpectrum)
xlim([0 400])
title('PSD')
  1 commentaire
dpb
dpb le 6 Mai 2020
Modifié(e) : dpb le 6 Mai 2020
t is undedfined
As noted in the other response to which you posted, the PSD via FFT WILL normalize peaks to the input magnitude as I showed there IFF the input frequency matches identically to a bin and is noise-free.
If the bin doesn't match that frequency exactly, then it is inevitable and unavoidable that some of the energy will be "smeared" into adjacent bins. You can minimize this by having more bins, but unless and until the two match identically there will always be some leakage.
Noise introduces similar issues -- it has energy content besides the input signal that will go wherever it needs go...again, unavoidable.
I haven't used the periodogram function enough to know precisely what it does internally, sorry.
But, the doc says...
"The units of the PSD estimate are in squared magnitude units of the time series data per unit frequency. For example, if the input data is in volts, the PSD estimate is in units of squared volts per unit frequency. For a time series in volts, if you assume a resistance of 1 Ω and specify the sample rate in hertz, the PSD estimate is in watts per hertz."
Pretty straightforward. So the magnitude in input units has to be integrated over the frequency bin and root taken to match the time series peak magnitude per the example PSD scaling. Which scaling to use is dependent upon the problem space and usage being made of the result.
But the same problems re: sampling and where the frequencies are computed are present; it's just a fact of discrete sampling.

Connectez-vous pour commenter.

Réponse acceptée

David Goodmanson
David Goodmanson le 6 Mai 2020
Modifié(e) : David Goodmanson le 6 Mai 2020
Hi John,
This concerns fft only; psd is a different animal.
For a continuous oscillation, most people do not seem to be concerned about fft output that spills into adjacent frequencies, but it's worth looking at as you are doing. With the fft you can get sharp one-point frequency peaks with a varying number of points and a varying amplitude, but not for time windows of arbitrary width. With fft It's required that there be an exact integral number of oscillations in the time window. And, using a single oscillation of a sine wave as an example, the last point in the array can't repeat the first point. The last point has to be "one point short" of a full oscillation
So
y = sin(2*pi*(0:9)/10) % correct but y = sin(2*pi*(0:10)/10) % incorrect
The code below gives two peaks, one at +160 kHz and one at -160 kHz for any number of fft points (above a minimum number), any number of complete oscillations in the time window and any amplitude. The abs(amplitude) of each peak is A/2, consistent with cos(w*t) = (1/2)(exp(i*w*t)+exp(-i*w*t)). If you want, as is common for plotting purposes, you could ignore the negative freqencies and double the height of the positive frequency peak.
f0 = 160e3;
ncyc = 23; % number of cycles in time window
A = 6; % amplitude
psi = pi/6; % arbitrary phase angle
N = 1299; % need N >2*ncyc to avoid aliasing
t0 = 1/f0;
T = ncyc*t0; % time window
delt = T/N;
t = (0:N-1)*delt;
y = A*cos(2*pi*f0*t+psi);
figure(1)
plot(t,y);grid on
delf = 1/(delt*N); % golden rule for fft
if rem(N,2) ==0
f = (-N/2:N/2-1)*delf;
else
f = (-(N-1)/2:(N-1)/2)*delf;
end
z = fftshift(fft(y))/N; % fftshift to put f=0 in the middle
figure(2)
stem(f,abs(z),'o-');grid on

Plus de réponses (0)

Catégories

En savoir plus sur Parametric Spectral Estimation dans Help Center et File Exchange

Produits


Version

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by