Why MATLAB's abs(fft(X)) returns maximum amplitude first harmonic's while in signal maximum is second's?

Question

Please let me ask you about my problem with Matlab code executing in R2022a release.

A signal has 3 harmonics with amplitudes 170, 220, and 150 in that order. Matlab's function abs(fft(X)) returns maximum amplitude for the first harmonic. What is the reason? Code is following:

SampFreq = 16000;
Segm = 1:2048;
Pitch = 45;
FirstHarmAngles = Pitch*2*pi/SampFreq*Segm+1.9*pi;
SinFirstHarmAngles = sin(FirstHarmAngles);
SecondHarmAngles = Pitch*2*2*pi/SampFreq*Segm+2.9*pi;
SinSecondHarmAngles = sin(SecondHarmAngles);
ThirdHarmAngles = Pitch*3*2*pi/SampFreq*Segm+0.3*pi;
SinThirdHarmAngles = sin(ThirdHarmAngles);
Xn = 170*SinFirstHarmAngles+220*SinSecondHarmAngles+150*...
SinThirdHarmAngles;
ABS = (abs(fft(Xn)))
plot(ABS(1:50))

I have asked similar question at Matlab but no answer yet. Observe that in the plot, the second amplitude is not conical as are the first and third.

Just when I posted this question, Mr. Jeffrey Clark posted his answer to my question at Matlab. Tough it's instructive I well appreciate any answer from you. Regards. — George Theodosiou, Sep 14 '22 at 14:16
His answer is the correct one, you don't need another one. Your frequency resolution is 16000 / 2048 ~= 8Hz so the FFT can't resolve 45Hz, 90Hz, and 135Hz. When increasing the fft length to 16384, your resolution is now ~= 1Hz, so your fft can resolve these frequencies better. — Jdip, Sep 14 '22 at 14:29
Jdip, please accept my many thanks for your comment. It's instructive too. Regards. — George Theodosiou, Sep 14 '22 at 15:06

Hilmar · Answer 1 · 2022-09-18T11:40:09.870

4

Spectral leakage.

In order for FFT of a sine wave to be a single line, the sine wave must have an integer multiple of periods in the FFT window.

Try the following

SampFreq = 16000;
Segm = 1:2048;
Pitch = 45;
% adjust to an integer period
freqResolution = SampFreq/length(Segm);
Pitch = round(Pitch/freqResolution)*freqResolution;
FirstHarmAngles = Pitch2pi/SampFreqSegm+1.9pi;
SinFirstHarmAngles = sin(FirstHarmAngles);
SecondHarmAngles = Pitch22pi/SampFreqSegm+2.9pi;
SinSecondHarmAngles = sin(SecondHarmAngles);
ThirdHarmAngles = Pitch32pi/SampFreqSegm+0.3pi;
SinThirdHarmAngles = sin(ThirdHarmAngles);
Xn = 170SinFirstHarmAngles+220SinSecondHarmAngles+150...
SinThirdHarmAngles;
ABS = (abs(fft(Xn)));
% plot with frequency axis and amplitude perserving scaling
n = length(Segm);
freqAxis = (0:n/2)/nfreqResolution;
plot(freqAxis(1:50),ABS(1:50)/length(Segm)*2);

edited Sep 18 '22 at 11:40

answered Sep 14 '22 at 15:44

Hilmar

44,604
1
32
63

1

I don’t think the OP was so worried about the spectrum consisting of spectral Lines as he was about the spectrum having the wrong amplitudes and inconsistent shapes at the frequencies of interest, but of course this answer addresses both concerns! – Jdip Sep 14 '22 at 15:48
@Jdip, please let me tell you that I'm just beginner in DSP and Matlab. So I can't understand your answer and comment under it. I choose - for the moment - Mr. Jeffrey Clark's suggestion and change fft(Xn) to fft(Xn,16384). I run your code in matlab and it responded doesn't recognize n, however for n = 2048 it responds correctly. Regards. – George Theodosiou Sep 17 '22 at 09:45
It's @Hilmar's code but yes, you want to add the line 'n=2048' on top (that's the length of 'Segm') – Jdip Sep 17 '22 at 16:24
good catch. Fixed – Hilmar Sep 18 '22 at 11:40

John Bofarull Guix · Answer 2 · 2022-09-20T14:56:59.047

The error in spectrum amplitudes of the 3 tones is caused by low resolution caused by the way fft is applied, and if the input signal is undersampled, which is not the case.

To show that there are fft parameters that increase amplitude and resolution, therefore achieving amplitude readings close to expected values, I have written for you the following MATLAB script.

If you find this answer useful would you please consider clicking on the accepted answer.

clear all;close all;clc  % clear the table

1.- Reproducing signal xn(t)

I rewrite your 3-tone signal with compact notation

fs=16000; % [Hz] sampling frequency     
dt=1/fs
at=2048  % amount time samples
nt=[0:at-1]; % time indexes
t=round(nt*dt,4);
f0=45;  % [Hz] base pitch
f1=f0*[1:3]; % all tone frequencies
ph1=round(1.9pi,4) % phase shifts, just 4 decimals
ph2=round(2.9pi,4)
ph3=round(.3*pi,4)
ph=[ph1 ph2 ph3];
a=2pif1'*t+ph'; % constant phases
s_ind=sin(a); % individual signals A=1
A=[170 220 150]  % amplitudes
% combined signal
xn=sum(A'.sin(a),1);
% same as
% xn=A(1)s1+A(2)s2+A(3)s3;

2.- Reproducing Xn(f)=fft(xn)

Xn=fft(xn);
absXn = abs(Xn);
% frequency window, as in question
nw_fft=50
figure(1)
ax1=gca
plot(ax1,absXn([1:nw_fft]))
grid on
xlabel('f')
title(ax1,['FFT(|Xn(f)|)    fs = ' num2str(fs) 'Hz, amount t samples = ' num2str(at)])

3.- Align frequency

when fft not supplied with amount frequencies then fft amount frequencies = amount time samples af : amount frequencies, default, or so it tries.

af=at
af=length(Xn)
f = fs*[0:af-1]/at;

4.- Remove repeated spectrum

real time signals (no quadrature) produce fft with repeated samples, Above fs/2 there on Xn=fft(xn) samples are mirrored, so let's remove them

Xn([floor(af/2):end])=[];
absXn([floor(af/2):end])=[];
f([floor(af/2):end])=[];

5.- Xn(f) amplitude correction

FFT amplitudes should be independent of amount input time samples

Xn=2*Xn/at;   % *2 to preserve signal energy
absXn=2*absXn/at;
figure(2)
ax2=gca
plot(ax2,f(1:nw_fft),absXn(1:nw_fft))
% plot(ax2,f,absXn)
hold(ax2,'on')
grid on
xlabel('f')
title(ax2,['FFT(|Xn(f)|)    fs = ' num2str(fs) ' amount t samples = ' num2str(at)])

detected frequencies and amplitudes

[pks,locs]=findpeaks(absXn(1:nw_fft),'NPeaks',3) 
fpeaks=(locs-1)*fs/at   % detected frequency peaks
plot(ax2,fpeaks,pks,'ro')
err_f=abs(fpeaks-f1)./(f1)100  % frequency error
err_absXn=abs(abs(pks)-abs(A))./abs(A)100 % amplitude error
str1=['f errors f1: ' ...
    num2str(round(err_f(1),2)) '% f2: ' ...
    num2str(round(err_f(2),2)) '% f3: ' ...
    num2str(round(err_f(3),2)) '%.']
text(ax2,max(xlim)/2,max(ylim)/2,str1)
str2=['A errors A1: ' ...
    num2str(round(err_absXn(1),2)) '% A2: ' ...
    num2str(round(err_absXn(2),2)) '% A3: ' ...
    num2str(round(err_absXn(3),2)) '%.']
text(ax2,max(xlim)/2,max(ylim)/2-10,str2)

6.- Improving frequency and amplitude resolution

% what is the current frequency resolution
f(1)
f(2)
% same as 
fs/at

f(2) > 1Hz. As pointed out in one of the posted comments the frequency resolution is too low to resolve these 3 tones with error below 1Hz.

If there are frequency errors, there are also amplitude errors.

Problem : Sampling really fast but not supplying enough time samples to the fft reduces frequency resolution.

Because fs=16000Hz is well above 2*3*f1 = 270Hz there may be overhead to play with.

Solution : Sampling slower, while still above Nyquist limit, may reduce fft frequency steps and therefore incrase frequency resolution.

Nyquist limit : 2 samples per cylce of signal highest frequency is the limit to avoid degrading signal by undersampling:

fs_nyquist=2*max(signal_f)

One must sample above this limit.

Every signal has it's own Nyquist limit, to avoid sampling below it or start losing signal

f_nyquist=3*2*f1(end)
% let's use 16 sample per cycle of signal highest frequency
fs2=8*f_nyquist  % now fs=2160
dt2=1/fs2
t2=nt*dt2;
a2=2pif1't2+ph';
xn2=sum(A'.sin(a2),1);
fft
Xn2=fft(xn2);
absXn2 = abs(Xn2);
new frequency reference
af2=length(Xn2)
f2 = fs2*[0:af2-1]/at;
fft length and amplitude corrections
Xn2([floor(af2/2):end])=[];
absXn2([floor(af2/2):end])=[];
f2([floor(af2/2):end])=[];
Xn2=2Xn2/at;  % Xn amplitude corrections 
absXn2=2absXn2/at;
figure(3)
ax3=gca
plot(ax3,f2([1:3*nw_fft]),absXn2([1:3*nw_fft]))
hold(ax3,'on')
grid on
xlabel('f')
title(ax3,['FFT(|Xn(f)|)    fs = ' num2str(fs2) 'Hz  amount t samples = ' num2str(at)])
detected frequencies and amplitudes
[pks2,locs2]=findpeaks(absXn2([1:3*nw_fft]),'NPeaks',3) % ,'SortStr','ascend'); 
fpeaks2=(locs2-1)*fs2/at   % detected frequency peaks
plot(ax3,fpeaks2,pks2,'ro')
err_f2=abs(fpeaks2-f1)./(f1)100  % frequency error
err_absXn2=abs(abs(pks2)-abs(A))./abs(A)100 % amplitude error
str1=['f1:' ...
    num2str(round(fpeaks2(1),2)) 'Hz  f2:' ...
    num2str(round(fpeaks2(2),2)) 'Hz  f3:' ...
    num2str(round(fpeaks2(3),2)) 'Hz']
text(ax3,max(xlim)/2,max(ylim)/2+20,str1)
str2=['A1:' ...
    num2str(round(pks2(1),2)) '  A2:' ...
    num2str(round(pks2(2),2)) '  A3:' ...
    num2str(round(pks2(3),2))]
text(ax3,max(xlim)/2,max(ylim)/2+10,str2)
str3=['f errors f1:' ...
    num2str(round(err_f2(1),2)) '%  f2:' ...
    num2str(round(err_f2(2),2)) '%  f3:' ...
    num2str(round(err_f2(3),2)) '%']
text(ax3,max(xlim)/2,max(ylim)/2,str3)
str4=['A errors A1:' ...
    num2str(round(err_absXn2(1),2)) '% A2:' ...
    num2str(round(err_absXn2(2),2)) '% A3:' ...
    num2str(round(err_absXn2(3),2)) '%']
text(ax3,max(xlim)/2,max(ylim)/2-10,str4)

now frequency error below 1Hz for all carriers yet amplitude errors are still significant

7.1.- Sampling faster while keeping same amount time samples

fs2=2*fs % 10*f_nyquist  % now fs=2160
% at=12048  % amount time samples
% nt=[0:at-1]; % time indexes
dt2=1/fs2
t2=nt*dt2;
a2=2pif1't2+ph';
xn2=sum(A'.sin(a2),1);
% fft
Xn2=fft(xn2);
absXn2 = abs(Xn2);
% new frequency reference
af2=length(Xn2)
f2 = fs2*[0:af2-1]/at;
% fft length and amplitude corrections
Xn2([floor(af2/2):end])=[];
absXn2([floor(af2/2):end])=[];
f2([floor(af2/2):end])=[];
Xn2=2Xn2/at;  % Xn amplitude corrections 
absXn2=2absXn2/at;
figure(4)
ax4=gca
plot(ax4,f2([1:3*nw_fft]),absXn2([1:3*nw_fft]))
hold(ax4,'on')
grid on
xlabel('f')
title(ax4,['FFT(|Xn(f)|)    fs = ' num2str(fs2) 'Hz  amount t samples = ' num2str(at)])
% detected frequencies and amplitudes
[pks2,locs2]=findpeaks(absXn2,'NPeaks',3) 
fpeaks2=(locs2-1)*fs2/at   % detected frequency peaks
plot(ax4,fpeaks2,pks2,'ro')
err_f2=abs(fpeaks2-f1)./(f1)100  % frequency error
err_absXn2=abs(abs(pks2)-abs(A))./abs(A)100 % amplitude error
str1=['f1:' ...
    num2str(round(fpeaks2(1),2)) 'Hz  f2:' ...
    num2str(round(fpeaks2(2),2)) 'Hz  f3:' ...
    num2str(round(fpeaks2(3),2)) 'Hz']
text(ax4,max(xlim)/2,max(ylim)/2+20,str1)
str2=['A1:' ...
    num2str(round(pks2(1),2)) '  A2:' ...
    num2str(round(pks2(2),2)) '  A3:' ...
    num2str(round(pks2(3),2))]
text(ax4,max(xlim)/2,max(ylim)/2+10,str2)
str3=['f errors f1:' ...
    num2str(round(err_f2(1),2)) '%  f2:' ...
    num2str(round(err_f2(2),2)) '%  f3:' ...
    num2str(round(err_f2(3),2)) '%']
text(ax4,max(xlim)/2,max(ylim)/2,str3)
str4=['A errors A1:' ...
    num2str(round(err_absXn2(1),2)) '% A2:' ...
    num2str(round(err_absXn2(2),2)) '% A3:' ...
    num2str(round(err_absXn2(3),2)) '%']
text(ax4,max(xlim)/2,max(ylim)/2-10,str4)

7.2.- Sampling faster, more time samples = more frequencies

fs2=10*fs
at=12048  % amount time samples
nt=[0:at-1]; % time indexes
dt2=1/fs2
t2=nt*dt2;
a2=2pif1't2+ph';
xn2=sum(A'.sin(a2),1);
% fft
Xn2=fft(xn2,at);
absXn2 = abs(Xn2);
% new frequency reference
af2=length(Xn2)
% f2 = fs2[0:af2-1]/at;
f2 = fs2[0:af2-1]/af2;
% fft length and amplitude corrections
Xn2([floor(af2/2):end])=[];
absXn2([floor(af2/2):end])=[];
f2([floor(af2/2):end])=[];
Xn2=2Xn2/at;  % Xn amplitude corrections 
absXn2=2absXn2/at;
figure(5)
ax5=gca
plot(ax5,f2([1:10*nw_fft]),absXn2([1:10*nw_fft]))
hold(ax5,'on')
grid on
xlabel('f')
title(ax5,['FFT(|Xn(f)|)    fs = ' num2str(fs2) 'Hz  amount t samples = ' num2str(at)])
% detected frequencies and amplitudes
[pks2,locs2]=findpeaks(absXn2, 'NPeaks',3,'MinPeakHeight',10)

fpeaks2=fs2*(locs2-1)/af2   % detected frequency peaks
plot(ax5,fpeaks2,pks2,'ro')
err_f2=abs(fpeaks2-f1)./(f1)100  % frequency error
err_absXn2=abs(abs(pks2)-abs(A))./abs(A)100 % amplitude error
str1=['f1:' ...
    num2str(round(fpeaks2(1),2)) 'Hz  f2:' ...
    num2str(round(fpeaks2(2),2)) 'Hz  f3:' ...
    num2str(round(fpeaks2(3),2)) 'Hz']
text(ax5,max(xlim)/2,max(ylim)/2+20,str1)
str2=['A1:' ...
    num2str(round(pks2(1),2)) '  A2:' ...
    num2str(round(pks2(2),2)) '  A3:' ...
    num2str(round(pks2(3),2))]
text(ax5,max(xlim)/2,max(ylim)/2+10,str2)
str3=['f errors f1:' ...
    num2str(round(err_f2(1),2)) '%  f2:' ...
    num2str(round(err_f2(2),2)) '%  f3:' ...
    num2str(round(err_f2(3),2)) '%']
text(ax5,max(xlim)/2,max(ylim)/2,str3)
str4=['A errors A1:' ...
    num2str(round(err_absXn2(1),2)) '% A2:' ...
    num2str(round(err_absXn2(2),2)) '% A3:' ...
    num2str(round(err_absXn2(3),2)) '%']
text(ax5,max(xlim)/2,max(ylim)/2-10,str4)

Now it looks as if there were again loss frequency/amplitude resolution but actually checking the errors one realises that the sketchy spectrum is accurate, because now both frequency and amplitude resolutions have enough small steps to measure accurately.

8.- Amplitude and frequency errors as function of sampling fs and at

at_range=2.^[6:.1:20];
fs_range=linspace(90,64000,280);
err_f_log=zeros(numel(fs_range),numel(at_range),3);
err_A_log=zeros(numel(fs_range),numel(at_range),3);
for k1=1:1:numel(fs_range)
    for k2=1:1:numel(at_range)
dt2=1/fs_range(k1);
nt=[0:at_range(k2)-1]; % time indexes
t2=nt*dt2;
% sin() inputs
a1=2pif1(1)t2+ph1;a2=2pif1(2)t2+ph2;a3=2pif1(3)*t2+ph3;
% individual signals A=1
s1=sin(a1);s2=sin(a2);s3=sin(a3);
% combined signal
xn=A(1)s1+A(2)s2+A(3)*s3;
% fft
Xn=fft(xn);Xn2=2*Xn2/at;  % Xn amplitude corrections
af=length(Xn);
f = fs_range(k1)*[0:af-1]/at_range(k2);
Xn([floor(af/2):end])=[];f([floor(af/2):end])=[];
absXn=abs(Xn);
% detected frequencies and amplitudes
[pks,locs]=findpeaks(absXn,'NPeaks',3,'MinPeakHeight',10);
fpeaks=(locs-1)*fs_range(k1)/at_range(k2);   % detected frequency peaks
% plot(ax4,fpeaks2,pks2,'ro') % same as findpeaks option 'Annotate',['peaks' | 'extents']
switch numel(pks)
    case 0
        err_f=[NaN NaN NaN];
        err_absXn=[NaN NaN NaN];
case 1  % assume single f detected then it's f0, although could be alias
    err_f(1)=abs(fpeaks-f1(1))./(f1(1))*100;
    err_f(2)=NaN;
    err_f(3)=NaN;
    err_absXn(1)=abs(abs(pks)-abs(A(1)))./abs(A(1))*100;
    err_absXn(2)=NaN;
    err_absXn(3)=NaN;

case 2
    err_f(1)=abs(fpeaks(1)-f1(1))./(f1(1))*100;
    err_f(2)=abs(fpeaks(2)-f1(2))./(f1(2))*100;
    err_f(3)=NaN;
    err_absXn(1)=abs(abs(pks(1))-abs(A(1)))./abs(A(1))*100;
    err_absXn(2)=abs(abs(pks(2))-abs(A(2)))./abs(A(2))*100;
    err_absXn(3)=NaN;

case 3
    err_f=abs(fpeaks-f1)./(f1)*100; 
    err_absXn=abs(abs(pks)-abs(A))./abs(A)*100; 

otherwise
    disp('findpeaks should return up to 3 peaks only')

end
% update logs

err_f_log(k1,k2,1)=err_f(1);
err_f_log(k1,k2,2)=err_f(2);
err_f_log(k1,k2,3)=err_f(3);
err_A_log(k1,k2,1)=err_absXn(1);
err_A_log(k1,k2,2)=err_absXn(2);
err_A_log(k1,k2,3)=err_absXn(3);
end

end

plot error surfaces as function of fs and at look for lowest point

% frequency error
figure(6)
ax6=gca
hs6=surf(err_f_log(:,:,1));
hs6.EdgeColor='none'
xlabel('fs [Hz]');ylabel('amount time samples');zlabel('f error [%]')
title(ax6,['f1 = ' num2str(f1(1)) 'Hz error [%]'])

figure(7)
ax7=gca
hs7=surf(err_f_log(:,:,2));
hs7.EdgeColor='none'
xlabel('fs [Hz]');ylabel('amount time samples');zlabel('f error [%]')
title(ax7,['f1 = ' num2str(f1(2)) 'Hz error [%]'])

figure(8)
ax8=gca
hs8=surf(err_f_log(:,:,3));
hs8.EdgeColor='none'
xlabel('fs [Hz]');ylabel('amount time samples');zlabel('f error [%]')
title(ax8,['f1 = ' num2str(f1(3)) 'Hz error [%]'])

amplitude error

figure(9)
ax9=gca
hs9=surf(err_A_log(:,:,1));
hs9.EdgeColor='none'
xlabel('fs [Hz]');ylabel('amount time samples');zlabel('A error [%]')
title(ax9,['A = ' num2str(A(1)) ' error [%]'])

figure(10)
ax10=gca
hs10=surf(err_A_log(:,:,2));
hs10.EdgeColor='none'
xlabel('fs [Hz]');ylabel('amount time samples');zlabel('A error [%]')
title(ax10,['A = ' num2str(A(2)) ' error [%]'])

figure(11)
ax11=gca
hs11=surf(err_f_log(:,:,3));
hs11.EdgeColor='none'
xlabel('fs [Hz]');ylabel('amount time samples');zlabel('A error [%]')
title(ax11,['A = ' num2str(A(3)) ' error [%]'])

9.1.- af up holding at constant changing frequency correction

nf=2.^[10:15];   % amount fft frequencies
figure
ax=gca
for k=1:1:numel(nf)
Xn=fft(xn,nf(k));
absXn=abs(Xn);
f = fs*[0:nf(k)-1]/at;

absXn([floor(nf(k)/2):end])=[];
f([floor(nf(k)/2):end])=[];
absXn=2*absXn/at;

plot(ax,f,absXn)
hold(ax,'on')
grid on
xlabel('f')


end
title(['|Xn(f)| with nf ' num2str(nf)]);

In this graph I did not update the frequency axis because it would get too crowded.

Make the amount of fft frequencies a power of 2, otherwise again loss of resolution happens shown in next point.

9.2.- af up holding at constant with more frequencies changing frequency correction

nf=floor(linspace(2048,8192,10));   % amount fft frequencies
figure
ax=gca
for k=1:1:numel(nf)
Xn=fft(xn,nf(k));
absXn=abs(Xn);
f = fs*[0:nf(k)-1]/at;

absXn([floor(nf(k)/2):end])=[];
f([floor(nf(k)/2):end])=[];
absXn=2*absXn/at;

plot(ax,f,absXn)
hold(ax,'on')
grid on
xlabel('f')


end
title(['|Xn(f)| with nf ' num2str(nf)]);

Just pushing more frequencies with larger nf in fft(xn,nf) keeping constant the amount of time samples deteriorates frequency resolution substantially, increasing side lobes.

nf same as at.

frequencies span and the frequency peaks with it because amount time samples at is kept constant.

If nf (nf same as at) is increased, so should amount time samples at

10.- Listening xn

fs=16000; % [Hz] sampling frequency     
dt=1/fs
at=2048  % amount time samples
nt=[0:at-1]; % time indexes
t=round(nt*dt,4);
f0=45;  % [Hz] base pitch
f1=f0*[1:3]; % all tone frequencies
ph1=round(1.9pi,4) % phase shifts, just 4 decimals
ph2=round(2.9pi,4)
ph3=round(.3*pi,4)
ph=[ph1 ph2 ph3];
% sin phases
a=2pif1'*t+ph';
% combined signal
xn=sum(A'.*sin(a),1);
% load example
xData=load('chirp.mat')
xObj=audioplayer(xData.y,xData.Fs)
play(xObj)
% replace data
xdata.y=xn
xData.Fs=fs
xObj.SampleRate=fs
xObj.TimerPeriod=t(end)
play(xObj)

the signal is the simluation of a bird tweet :)

Thanks for reading my answer.

John Bofarull, please let me express you my sincere gratitude for this great work you did for me. I'm beginner in DSP and Matlab, but I hope some day understand your answer. For the moment I choose Mr. Jeffrey Clark's suggestion (at Matlab) and just change fft(Xn) to fft(Xn,16384). Regards. — George Theodosiou, Sep 21 '22 at 07:52

Why MATLAB's abs(fft(X)) returns maximum amplitude first harmonic's while in signal maximum is second's?

2 Answers2

Linked