Why do I receive results mismatch in the FFT signal while using MATLAB built in FFT function? How can I able to rectify this?
Show older comments
Hi
I am imposing a signal such as 60cos(21t) + 40cos(42t) + 20cos(63t). So I am supposed to get the following amplitude and frequency from the FFT signal.
1) Frequency (21), Amplitude (60)
2) Frequency (42), Amplitude (40)
3) Frequency (63), Amplitude (20).
The following script used to obtain the above FFT results:
%% FFT Macro
Start_time = 0;
End_time = 10;
dt=0.001;
Fs=(1/dt)/2; % Sampling frequency/2
range = End_time-Start_time;
FFT_zone = 1/range; % frequency resolution
NPL = ((range)/dt)+1; % Number of points for local run
IS = (Start_time/dt)+1;
Ns = IS+NPL-1;
Ns1 = (((range)/dt)-150);
Frequency_vector=0:FFT_zone:Fs;
Rin.Frequency_vector=Frequency_vector';
columns = 2;
for i=1:columns-1
eval(['Vector_selected1=Rin.signal(:,' num2str(i) ');']); % selection of the run #i
eval(['Vector_selected2=Vector_selected1(' num2str(IS) ':' num2str(Ns) ');']); % selection of the time range for FFT
eval(['FFT_wzd_Run_signal=abs(fft(Vector_selected2))/Ns1;']); % FFT calculated in magnitude ( amplitude Y - axis)
end
[r,c] = size(Rin.Frequency_vector);
Rin.xx11=FFT_wzd_Run_signal(1:r);
Rin.xx11(1,1) = 0;
%% Just for reference ( i have not given complete steps of plotting macro, error is not in this point)
plot(Rin.t,Rin.signal,'r','LineWidth',2); % subplot 1
plot(Rin.Frequency_vector,(Rin.xx11*2),'r','LineWidth',2); % subplot 2
The following signals are obtained by using the macro.
SIGNAL 1 SIGNAL 2
In the first signal, i am getting the FFT amplitude are almost closer to my expectation. But in the second signal, the FFT amplitudes are deviates from my expectation. In Inputs, The only deviation between both signals are start time and end time.
For signal 1 --> Start time = 0, End time = 10
For signal 2 --> Start time = 4, End time = 6.
Total available signal time = 10 s.
Eventhough both siganl are same, due to the zoomed windowing time leads to the deviation in FFT amplitudes. So i need to avoid this issue.
Kindly help me out on this.
Thanks in advance.
Answers (2)
dpb
on 4 Jun 2021
0 votes
If the fundamental frequency doesn't precisely match the frequency bins of the FFT, then the energy will be spread over multiple bins around that center frequency. To get the full power of the peak in that case you must integrate the peak, not just take the value of the peak at the peak frequency.
A few prior Answers may be of interest as well---
Somewhere, sometime, I illustrated the effect of changing the sampling frequency so the output frequency doesn't quite match the input frequencies as does in the example at doc fft but I don't seem to find that particular Answer at the moment. It's not hard to do however, just change the input frequency in the generated signal to not be an exact match of one of the frequency bins with the same sampling parameters or change the sample rate a little to shift where the frequency bins lie and observe the change from a single bin spike to one that has energy in several bins.
18 Comments
Gopinath Karuppannan
on 15 Jun 2021
dpb
on 15 Jun 2021
Only by having noise-free input signal that is at an exact center frequency of a frequency bin will the single-bin point value match identically -- that's the whole point.
You have to integrate the peak otherwise to get the whole energy contained across multiple bins and there is inevitable leakage that is unavoidable owing to the introduction of noise.
Gopinath Karuppannan
on 15 Jun 2021
As I noted before and in the above answers, refer to the doc for FFT to see how to more neatly get the needed frequency vector.
And, you left out the the 2pi factor in generating your signals so they're not really at those frequencies...
Fs = 1000; % the sampling frequency
T = 1/Fs; % the time sample dt
L = 10*Fs; % generate 10 seconds of data
t = (0:L-1)*T; % the time vector to match
f=Fs*(0:(L/2))/L; % and the frequency vector for baseband
S=60*cos(2*pi*21*t)+40*cos(2*pi*42*t)+20*cos(2*pi*63*t);
Y=fft(S);
P2=abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
plot(f,P1) % gives plot of three spikes of mag 60,40,20 @ 21,42,63 Hz
ix4=find(t==4); % which point is the 4 second one...
Y=fft(S(ix4:end)); % use that point on
L2=numel(t(ix4:end)); % get the length of the signal for normalization
P2=abs(Y/L2); % do same as above but with new L
P1 = P22(1:L2/2+1);
P12=P1;P12(2:end-1) = 2*P1(2:end-1);
f2=Fs*(0:(L2/2))/L2; % have to recompute f to go along, too...
hold on
plot(f2,P12);
The result is a plot with only one set of peaks, the two side-by-side look like
The key is there is a bin that is identically equal to the expected frequency and the number of elements in the signal is used to normalize and the frequency vector and signal are consistenty in generation.
Paul
on 17 Jun 2021
Regarding this line:
P1(2:end-1) = 2*P1(2:end-1);
Why does the indexing start at 2 and not at 1?
dpb
on 17 Jun 2021
The DC frequency is present in the returned FFT only once and contains the total DC magnitude already. (Ditto for Nyquist for the -1 on the end as well).
But the area under the (interpolated) curve between P1(1) and P1(2) will be too low, won't it? I thought the point of the multiply-by-2 was to try to make the area under the (interpolated) curve of a half period of the magnitude of the scaled DFT equal to the area under the curve of the full period of the magnitude of the DTFT, at least to within a factor of L (which is another scaling on which I'm unclear why it's used, in general).
Nope...as Goldilocks says, this one is just right!
NB: the peaks above are identical to within machine-precision of the magnitude of the input time signal for each peak. Of course that was a zero-mean signal, but if add a DC component,
>> Fs=1000;
>> T = 1/Fs; % the time sample dt
L = 10*Fs; % generate 10 seconds of data
t = (0:L-1)*T; % the time vector to match
f=Fs*(0:(L/2))/L; % and the frequency vector for baseband
S=pi+60*cos(2*pi*21*t)+40*cos(2*pi*42*t)+20*cos(2*pi*63*t);
Y=fft(S);
P2=abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
>> P1(1)
ans =
3.14
>>
is exactly (well, again, to machine precision) pi --
>> P1(1)-pi
ans =
7.1054e-15
>>
Since the DC component is only in the output vector one time, the internal algorithm takes care of it for you automagically; you just don't want to double it.
As for the divisor L, that's again because of the definition of the DFT as implemented in FFTW that MATLAB uses.
For the above note:
>> sum(abs(Y))
ans =
1.2314e+06
>> sum(abs(Y))/L
ans =
123.1416
>> 60+40+20+pi
ans =
123.1416
>>
Or, the magnitude of each peak is L times the peak value of the time signal frequency...
One mistake one sees is that folks use NFFT as the divisor when augmenting the time signal with nextpow2 points where the NFFT-L added points are nothing but a zero-valued augmentation of the signal and so contain no energy. Thus, this will ratio down the observed amplitude by that factor -- the normalization should always only use L, the length of the actual signal.
I understand all of this, but I'm still not convinced (thought still open to being so) that it's a compelling argument. I wouldn't scale by L or by 2.
Let me use a different sequence to illustrate why I say that.
n = 0:5;
x = 10 - n; % a six point, finite duration sequence
L = numel(n);
All of the information in x is contained in one period of its DTFT:
w = 0:.01:2*pi;
DTFTX = sum(x.*exp(-1j*w(:).*n),2);
figure;plot(w,abs(DTFTX)),grid
Now I can take the 6-point DFT:
nfft = 6; X6 = fft(x,nfft); w6 = (0:nfft-1)/nfft*2*pi;
figure;plot(w,abs(DTFTX));
hold on
stem(w6,abs(X6));
As expected, the DFT sequence is sampling the DTFT. Though X6 does tell me everything I need to know about x (after all, I can reconstruct x using ifft() ), it's kind of coarse with six points. So we zero pad to interpolate in the frequency domain.
nfft = 30; X30 = fft(x,nfft); w30 = (0:nfft-1)/nfft*2*pi;
figure;plot(w,abs(DTFTX));
hold on
stem(w30,abs(X30));
So far so good. X30 is giving me a nice visual representation of X. And of course I can recover x from X30 using ifft() and taking the first L points. Furthermore, I can do some operations on X30 (like use it to approximate the energy in x after I add a point at 2*pi).
Now I'll do the procedure outlined above:
nfft = 30;
P2 = abs(fft(x,nfft)/L);
P1 = P2(1:nfft/2+1);
P1(2:end-1) = 2*P1(2:end-1);
figure;plot(w,abs(DTFTX));
hold on
stem(2*pi*(0:(nfft/2))/nfft,P1);
I look at that last stem plot and fail to see how I would use it in comparison to the others. Not only does it not approximate the DTFT of x, it's actually misleading because values of P1 at 0 and pi are too small relative to their neighbors.
One can normalize/interpret in many different ways depending on the need. Here OP wanted the FFT spectra to be normalized to match the input time wave magnitude at the input frequencies -- the normalization given does that -- no more, no less. Nothing else is claimed.
How would you answer the OP's desire without either the factor of L or 2?
ADDENDUM:
I'd have to dig deeper into why your above didn't return the expected DC; with the normalization as I would do it works just as well (albeit, of course, if one doesn't match up the number of points precisely, the nearest frequency bin won't be identically integer-valued so you'll get some peak-spreading, but...
NFFT=16384; % something other than NFFT==L
Y=fft(S,NFFT); % FFT over N
P2=abs(Y/L); % but the power is only in the L signal points
P1 = P2(1:NFFT/2+1); %
P1(2:end-1) = 2*P1(2:end-1);
P1(1)
f=Fs*(0:(NFFT/2))/NFFT; % and the frequency vector for baseband
[pk,fr]=findpeaks(P1,f,'MinPeakHeight',10)
ans =
3.1416
pk =
59.8487 39.5715 19.4981
fr =
20.9961 41.9922 62.9883
>>
shows the above normalization still works precisely as would want it to -- the peaks are no longer identical to the time signal to machine precision because the frequencies no longer match the input signal frequencies identically, but they're close because there's enough detail to get almost where the bin centers need to be.
NB: the DC component still comes out at the DC level of the signal.
NB2: When I was a half-mile down in the bowels of the paper mill, if there were a DC component to the signal that would have indicated I was in serious trouble as the whole line would be translating on its way towards the Tennessee River! :) Needless to say, DC was of essentially no interest to us other than ensuring didn't have some sort of bias signal on the inputs...
But, for those applications where it might possibly have some meaning/use, this will get it right.
The DC didn't show up well so I stuck an asterisk there just so it would show...
Paul
on 18 Jun 2021
I think this statement is key: "One can normalize/interpret in many different ways depending on the need."
If I wanted to manipulate the output of fft() if the need is the plot desired by the OP for this particular question, I'm sure I would have done the same thing (or at least very similar).
But I feel like I see take-the-first-half-divide-by-L-and-multiply-by-2-except-at-DC-and-pi for cases where I don't think it's needed and might even be misleading.
I fail to see what cases that might be for someone to mistake the use -- the application is clear to make the one-sided PSD from the MATLAB FFT that is two-sided and normalize, one does NOT want to double either the DC nor the Nyquist frequency bins. That's true whether it is to normalize to the amplitude of the time signal here or to the power instead; the energy of the signal is all in those two bins already by the FFTW algorithm (as is the scaling by L, the length of the input signal containing the real time signal).
If you generate the DFT in some other fashion, then obviously you'll need to treat it commensurately with whatever it is that is going to be consistent with that development; but that's hardly ever the Q? raised in Answers and certainly wasn't here.
It's a red herring to raise the point in my view as it has no bearing on the actual use case, nor does the example documentation show other uses where it would not be appropriate.
dpb
on 19 Jun 2021
That's ok, I think I couldn't see where you were coming from and still don't fully understand how there can be any confusion simply using the results of FFT() on a sampled signal although it is a point not specifically documented about the unique nature of DC and Nyquist and so one often sees the blanket factor of 2 on the full half-spectrum.
But, often one will have used dtrend first, anyway, and so the DC component is zero anyway and since 2*0 ==> 0 it never shows up...or it does get doubled but nobody cares because they don't ever look at the DC component, anyway.
Similarly for the Nyquist frequency bin, only extremely rarely would one ever look at the actual magnitude there; if there's significant energy content showing up there, the signal was badly undersampled for any real analysis and needs fixing for that problem, anyways...
Those are practical lessons; but doesn't conflate the issue of the representation of the DTFT with the FFT that your example does. If one applies in an apples:apples comparison, your case works just as well although you've not defined a sampling frequency to go along with the sampled signal that woud be need to get actual frequency components;
n = 0:5;
x = 10 - n; % a six point, finite duration sequence
L = numel(n);
plot(x,'*-')
Y=fft(x);
P2=abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
>> plot(P1,'*-')
>> mean(x) ans =
7.5000
>> P1(1)
ans =
7.5000
>>
It's pretty-much meaningless, of course, but the FFT() in MATLAB applied to the signal and converted into a one-sided PSD returns the expected results irrespective of the shape or length of the input signal. Whether it makes any sense to do so or not is a totally different question, of course! :)
Your comparison is just looking at two completely different things and so should not be surprising as they two don't compare.
Paul
on 19 Jun 2021
For NFFT an even number, I can see why one might want to do the one sided thing and scale P1(2:end-1) by 2, or better better yet, sqrt(2). I'm still not sure of the utility of scaling everything by L in general.
With regards to your last statement "Your comparison ..." can you be more specific about which comparison of mine you're referring to and what are the two completely different things in that comparison?
dpb
on 19 Jun 2021
" utility of scaling everything by L in general."
That's precisely what I've pointed out you seem to keep ignoring -- I've never said anything about "everything in general". This is applies to the specific normalization of a PSD to reproduce the time signal amplitude from the FFT as implemented in FFTW that MATLAB uses.
For odd NFFT, use 1:fix(NFFT/2); the DC component is still in the first location and still contains the mean but the tail end is reflected so use 2*P1(2:end). (Or not, your choice depending on the output scaling desired--but, if you do use one-sided instead of two, then the total magnitude is half the integral of the two-sided and there is still only the one DC component so you'll be off by that much if you double it. You'll still see that mean(S)==P1(1) identically even or odd NFFT after dividing by L).
As for scaling by L or not, again it all depends on what interpretation one wants for the output -- in my world, the actual acceleration or displacements are wanted in g/Hz or sqrt(Hz) so if is normalized to the magnitude of the input signal, then scaling by the accelerometer sensitvity in mV/g returns g. Returning to the DC component, if one doesn't divide by L, then the DC component will be L*mean(S) which doesn't reflect that magnitude. But, agreed, it is totally arbitrary; it simply is a choice of what scaling/units one wants the result to be at/in.
Paul
on 19 Jun 2021
Apologies for the poor wording on my part. I diidn't mean to imply that you said to always scale by L in general. Rather, I was referring to my anecdotal recollection that a lot of other on Answers do include that scaling in a variety of general problems where it's unclear to me how that scaling aids in the analysis of those problems (unlike the specific problem in this Question where that scaling makes the amplitude plot in the frequency domain match the signal amplitudes in the time domain).
While it may not help, I doubt it hurts, either... :)
I don't otomh have any idea how to interpret the results without, though...what does it mean unscaled other than purely relative when the magnitudes are dependent upon the length of the input series? What specific analysis would be better for it, you think?
For example, say the input signal of OP's were the input from an accelerometer from a perfectly noise-free machine and if we had sampled for 1 second or, wanting more averaging-out of any noise that might be present, 2 seconds. If we computed the PSD for those two measurements and didn't normalize by the sample length, we would get something like:
Now, while both of those are the same in general (although we'll note the peak spreading in the first owing to the frequency binning not being really well matched with the lower frequency so the actual peak itself is quite under-estimated, if we do integrate we'll be able to estimate the total energy reasonably well), the second tells us we've got something much more dastardly going on in the machine than does the first--but that's not so, it's the same signal.
How can that result be desirable by not normalizing the FFT by the sample length to compensate for the fact the FFT algorithm itself returns the values scaled by the number of samples?
I grant there being alternatives of whether to then scale to the peak or rms value, etc., etc., as to what specifically is the actual scaling/units, but I just can't see the objection to putting both on the same footing with the factor L.
your signal does not have any noise. The way to correct the problem is to correct your time vector. You seem to have:
dt=0.001; t=[0 : dt : 10];
This means that you have 10001 miliseconds. Each time point represents itself and the interval in front of it. The correct way to have 10 seconds (and make your waves have an exact integer number of periods in your signal) is this:
dt=0.001; t=[0 : dt : 10-dt];
Edit your code accordingly.
Categories
Find more on Spectral Measurements in Help Center and File Exchange
on 4 Jun 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
