Why do I receive results mismatch in the FFT signal while using MATLAB built in FFT function? How can I able to rectify this?

Question

Gopinath Karuppannan 2021-6-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/847665-why-do-i-receive-results-mismatch-in-the-fft-signal-while-using-matlab-built-in-fft-function-how-ca

编辑： pepe 2021-6-22

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

dpb 2021-6-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/847665-why-do-i-receive-results-mismatch-in-the-fft-signal-while-using-matlab-built-in-fft-function-how-ca#answer_716825

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---

https://www.mathworks.com/matlabcentral/answers/356376-can-someone-explaing-the-computation-for-double-sided-and-single-sided-spectrum-in-fft-example?s_tid=srchtitle#answer_281307

https://www.mathworks.com/matlabcentral/answers/429086-how-to-calculate-psd-estimate-of-emg-data-in-1hz-frequency-bins-for-100ms-time-windows?s_tid=srchtitle#answer_346358

https://www.mathworks.com/matlabcentral/answers/395585-normalize-fft-signal-of-ecg-data?s_tid=srchtitle#answer_315674

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 个评论
显示 16更早的评论隐藏 16更早的评论

dpb 2021-6-15

编辑：dpb 2021-6-15

在 MATLAB Online 中打开

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.

dpb 2021-6-17

编辑：dpb 2021-6-17

在 MATLAB Online 中打开

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.

Paul 2021-6-18

编辑：Paul 2021-6-18

在 MATLAB Online 中打开

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.

dpb 2021-6-18

编辑：dpb 2021-6-18

在 MATLAB Online 中打开

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...

dpb 2021-6-19

在 MATLAB Online 中打开

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.

dpb 2021-6-19

" 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 2021-6-19

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).

dpb 2021-6-19

编辑：dpb 2021-6-21

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.

请先登录，再进行评论。

Answer 2

pepe 2021-6-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/847665-why-do-i-receive-results-mismatch-in-the-fft-signal-while-using-matlab-built-in-fft-function-how-ca#answer_730275

编辑：pepe 2021-6-22

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.

https://www.mathworks.com/matlabcentral/answers/847665-why-do-i-receive-results-mismatch-in-the-fft-signal-while-using-matlab-built-in-fft-function-how-ca#comment_1584380

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Why do I receive results mismatch in the FFT signal while using MATLAB built in FFT function? How can I able to rectify this?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

18 个评论
显示 16更早的评论隐藏 16更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Why do I receive results mismatch in the FFT signal while using MATLAB built in FFT function? How can I able to rectify this?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

18 个评论 显示 16更早的评论隐藏 16更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

18 个评论
显示 16更早的评论隐藏 16更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论