Wrong FFT for an audio file

I want to sketch the power spectrum of an audio file but i get wrong answer.human speech should be in range of 50 to 300 hz

[x Fs] = audioread('v0.mp3');
nf=length(x);
Y = fft(x,nf);
Y = Y-mean(Y);
f = Fs/2*linspace(0,1,nf/2+1);
plot(f,abs(Y(1:nf/2+1)));

i should get this:

but instead i get this:

[Walter Roberson]

Why are you subtracting the mean() of the fft from the fft results? It would make a lot more sense to have subtracted the mean of x from x

Answer 1

[Star Strider]

The plot is correct. You are not considering the exponential multiplication at the right end of the frequency axis.

This makes it a bit more obvious:

[x Fs] = audioread('Amirhosein Khanlari v0.mp3');
nf=size(x,1);
Y = fft(x-mean(x))/nf;
f = Fs/2*linspace(0,1,fix(nf/2)+1);

figure
plot(f,abs(Y(1:nf/2+1))*2)
xlim([0 1.5E+4])
set(gca, 'XTick', (0:2500:15000))

producing:

[Amirhosein Khanlari]

thanks man!and can you tell me how can i determine the frequency that has the max power?

[Star Strider]

As always, my pleasure!

Yes. The code changes to:

[x Fs] = audioread('Amirhosein Khanlari v0.mp3');
nf=size(x,1);
Y = fft(x-mean(x))/nf;
f = Fs/2*linspace(0,1,fix(nf/2)+1);

[maxAmp,maxAmpidx] = max(abs(Y(1:nf/2+1))*2);           % Maximum Amplitude & Index

figure
plot(f,abs(Y(1:nf/2+1))*2)
xlim([0 1.5E+4])
set(gca, 'XTick', (0:2500:15000))
text(f(maxAmpidx), maxAmp, sprintf('\\leftarrow %.1f Hz = %.5f Amplitude', f(maxAmpidx), maxAmp), 'HorizontalAlignment','left')

and produces:

Here since you only want to maximum amplitude, use the max function (with two outputs). To find more than one, use the findpeaks function, or the islocalmax function.

Also note that this code will give you the amplitude of the signal at each frequency, not the power. To calculate and plot power, square the amplitude:

Y = fft((x-mean(x)).^2)/nf;

[Amirhosein Khanlari]

if i use the power the maximum code wont work:(

[x Fs] = audioread('Amirhosein Khanlari v0.mp3');
nf=size(x,1);
Y = fft((x-mean(x)).^2)/nf;
f = Fs/2*linspace(0,1,fix(nf/2)+1);
[maxAmp,maxAmpidx] = max(abs(Y(1:nf/2+1))*2);           % Maximum Amplitude & Index
figure
plot(f,abs(Y(1:nf/2+1))*2)
xlim([0 1.5E+4])
set(gca, 'XTick', (0:2500:15000))
text(f(maxAmpidx), maxAmp, sprintf('\\leftarrow %.1f Hz = %.5f Amplitude', f(maxAmpidx), maxAmp), 'HorizontalAlignment','left')

[Star Strider]

That is interesting. Subtracting the mean should produce a zero D-C offset. It does in the Fourier transform, however it does not when the time-domain signal is squared first.

According to Parseval’s Theorem, it is correct to square the Fourier transform or the original time-domain signal. Squaring the Fourier transformed signal to create ‘Psd’:

[x Fs] = audioread('Amirhosein Khanlari v0.mp3');
nf=size(x,1);
Y = fft(x-mean(x))/nf;
f = Fs/2*linspace(0,1,fix(nf/2)+1);
Psd = (abs(Y(1:nf/2+1))*2).^2;
[maxAmp,maxAmpidx] = max(Psd);              % Maximum Amplitude & Index

figure
plot(f,Psd)
xlim([0 1.5E+4])
set(gca, 'XTick', (0:2500:15000))
text(f(maxAmpidx), maxAmp, sprintf('\\leftarrow %.1f Hz = %.2E Amplitude', f(maxAmpidx), maxAmp), 'HorizontalAlignment','left')

This produces the correct result.