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:
1 个评论
Walter Roberson
2020-1-21
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
采纳的回答
Star Strider
2020-1-21
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:
4 个评论
Amirhosein Khanlari
2020-1-22
Star Strider
2020-1-22
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:
%20-%202020%2001%2021.png)
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
2020-1-26
Star Strider
2020-1-26
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.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Spectral Measurements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
