这是我用来测试MATLAB基本声音分析函数的脚本,包括
fft()
的输出使用
fftshift()
进行显示。
if ~exist('inputFile', 'var')
inputFile = 'vibe.wav';
end
[inputBuffer, Fs] = audioread(inputFile);
fileSize = length(inputBuffer);
numSamples = 2.^(ceil(log2(fileSize)));
x = zeros(numSamples, 1);
x(1:fileSize) = inputBuffer(:,1);
clear inputBuffer;
clear fileSize;
t = linspace(0, (numSamples-1)/Fs, numSamples)';
f = linspace(-Fs/2, Fs/2 - Fs/numSamples, numSamples)';
X = fft(x);
plot(t, x);
xlabel('time (seconds)');
ylabel('amplitude');
title(['time-domain plot of ' inputFile]);
sound(x, Fs);
pause;
plot(f, real(fftshift(X)));
xlabel('frequency (Hz)');
ylabel('real part');
title(['real part frequency-domain plot of ' inputFile]);
pause;
plot(f, imag(fftshift(X)));
xlabel('frequency (Hz)');
ylabel('imag part');
title(['imag part frequency-domain plot of ' inputFile]);
pause;
plot(f, abs(fftshift(X)));
xlabel('frequency (Hz)');
ylabel('amplitude');
title(['abs frequency-domain plot of ' inputFile]);
pause;
plot(f, 20*log10(abs(fftshift(X))+1.0e-10));
xlabel('frequency (Hz)');
ylabel('amplitude (dB)');
title(['dB frequency-domain plot of ' inputFile]);
pause;
semilogx(f(numSamples/2+2:numSamples), 20*log10(abs(X(2:numSamples/2))));
xlabel('frequency (Hz), log scale');
ylabel('amplitude (dB)');
title(['dB vs. log freq, frequency-domain plot of ' inputFile]);
pause;
semilogx(f(numSamples/2+2:numSamples), (180/pi)*angle(X(2:numSamples/2)));
xlabel('frequency (Hz), log scale');
ylabel('phase (degrees)');
title(['phase vs. log freq, frequency-domain plot of ' inputFile]);
pause;
semilogx(f(numSamples/2+2:numSamples), unwrap(angle(X(2:numSamples/2))));
xlabel('frequency (Hz), log scale');
ylabel('unwrapped phase (radians)');
title(['unwrapped phase vs. log freq, frequency-domain plot of ' inputFile]);
如果您正在将音频分段窗口处理并将其传递给FFT,则应在FFT的输入上使用fftshift()
,以将窗口化段的中心定义为t = 0点。
[x_input, Fs] = audioread('vibe.wav');
N = 2*floor(length(x_input)/2);
x = x_input(1:N);
t = linspace(-N/2, N/2-1, N);
omega = linspace(-pi, pi*(1-2/N), N);
X = fftshift( fft( fftshift( x.*hamming(length(x)) ) ) );
[X_max k_max] = max( abs(X) );
figure(1);
plot(t, x, 'g');
figure(2);
plot(omega, abs(X), 'b');
hold on;
plot(omega(k_max), X_max, 'or');
hold off;