在-fs/2至fs/2的范围内绘制WAV文件的幅度和相位谱 [英] Plotting the magnitude and phase spectra of a wav file in the range of -fs/2 to fs/2

问题描述

我在绘制wav文件的FFT时遇到问题.我设法绘制了信号的幅度和相位频谱，但是我需要在范围-fs/2:fs/2中重复此过程.

I'm having problems plotting the FFT of a wav file. I managed to plot the magnitude and phase spectrums of the signal, however I need to repeat this in range -fs/2:fs/2.

%read sound files 
%'y' is the vector holding the original samples & 'fs' refers to the sampling frequency 
[y,fs] = wavread('handel.wav'); 
ydft = fft(y); %fft to transform the original signal into frequency domain 
n = length (y); %length of the original signal 
%  y has even length 
ydft = ydft(1:length(y)/2+1); 
% create a frequency vector 
freq = 0:fs/length(y):fs/2; 
shiftfreq = fftshift(freq); 

%plot original signal in time domain; 
figure; 
plot ((1:n)/fs, y); 
title('handel.wav in time domain'); 
xlabel ('second'); 
grid on; 

% plot magnitude in frequency domain 
figure; 
plot(freq,abs(ydft)); 
title('handel.wav in frequency domain'); 
xlabel ('Hz'); 
ylabel('Magnitude'); 
grid on; 

% plot phase in frequency domain 
figure; 
plot(freq,unwrap(angle(ydft))); 
title ('handel.wav in frequency domain'); 
xlabel ('Hz'); 
ylabel ('Phase'); 
grid on; 

推荐答案

您当前正在做的是绘制半频谱，所以从0 <= f < fs/2开始，其中fs是信号的采样频率，所以是奈奎斯特频率.请注意，仅当信号为实数时才考虑半频谱.这意味着负光谱与正光谱是对称的，因此您实际上不需要在这里考虑负光谱.

What you are currently doing now is plotting the half spectrum, so from 0 <= f < fs/2 where fs is the sampling frequency of your signal, and so fs/2 is the Nyquist frequency. Take note that considering the half spectrum is only valid if the signal is real. This means that the negative spectra is symmetric to the positive spectra and so you don't really need to consider the negative spectra here.

但是，您想要绘制幅度和相位的整个频谱.请注意，使用MATLAB计算fft时，它会使用 Cooley- Tukey 算法，因此在计算N点FFT时，结果的一半是从0 Hz(包括0Hz)到fs/2 Hz的频率，另一半是从-fs/2 Hz(包括两端)的频率到0 Hz独占.

However, you would like to plot the full spectrum of the magnitude and phase. Take note that when calculating the fft using MATLAB, it uses the Cooley-Tukey algorithm so when computing the N point FFT, half of result is for the frequencies from 0 Hz inclusive up to fs/2 Hz exclusive and the other half is for the frequencies from -fs/2 Hz inclusive up to 0 Hz exclusive.

因此，要绘制完整频谱，只需对完整信号执行fftshift，以便频谱的右半部分和左半部分互换，以使0 Hz频率位于信号的中心.另外，您必须生成-fs/2至fs/2之间的频率才能覆盖整个频谱.具体来说，您需要生成在-fs/2至fs/2之间线性间隔的N点.但是，请注意，最后排除了fs/2 Hz处的奈奎斯特频率，因此您需要在-fs/2至fs/2之间生成N+1点，并删除最后一个点，以便在每个频点都正确.生成点的线性数组的最简单方法是使用 linspace 命令，其中开始频率为-fs/2，结束频率为fs/2，并且您希望N+1点在此范围之间并删除最后一个点:

As such, to plot the full spectrum, simply perform a fftshift on the full signal so that the right half and left half of the spectrum is swapped so that the 0 Hz frequency is located in the centre of the signal. Also, you must generate frequencies between -fs/2 to fs/2 to cover the full spectrum. Specifically, you need to generate N points linearly spaced between -fs/2 to fs/2. However, take note that the Nyquist frequency at fs/2 Hz is being excluded at the end, so you need to generate N+1 points between -fs/2 to fs/2 and remove the last point in order for the right step size between each frequency bin to be correct. The easiest way to generate this linear array of points is by using the linspace command where the start frequency is -fs/2, the ending frequency is fs/2 and you want N+1 points between this range and remove the last point:

freq = linspace(-fs/2, fs/2, n+1);
freq(end) = [];

这样，借用了您的代码的某些部分，这就是修改后的代码的样子，以绘制幅度和相位的全部频谱:

As such, borrowing some parts of your code, this is what the modified code looks like to plot the full spectrum of the magnitude and phase:

%// Read in sound file
[y,fs] = wavread('handel.wav'); 

%// Take N-point FFT where N is the length of the signal
ydft = fft(y);
n = numel(y); %// Get N - length of signal

%// Create frequency vector - make sure you remove last point
freq = linspace(-fs/2, fs/2, n+1);
freq(end) = [];

%// Shift the spectrum
shiftSpectrum = fftshift(ydft);

%//plot original signal in time domain; 
figure; 
plot ((0:n-1)/fs, y); %// Note you should start from time = 0, not time = 1/fs
title('handel.wav in time domain'); 
xlabel ('second'); 
grid on; 

%// plot magnitude in frequency domain 
figure; 
plot(freq,abs(shiftSpectrum)); 
title('handel.wav in frequency domain'); 
xlabel ('Hz'); 
ylabel('Magnitude'); 
grid on; 

%// plot phase in frequency domain 
figure; 
plot(freq,unwrap(angle(shiftSpectrum))); 
title('handel.wav in frequency domain'); 
xlabel('Hz'); 
ylabel('Phase'); 
grid on; 

我无权访问您的handel.wav文件，但我将使用MATLAB随附的文件.您可以使用load handel;将其加载.采样频率存储在名为Fs的变量中，因此在我上面编写的代码起作用之前，我必须做fs = Fs;.此特定文件的采样频率为8192 Hz，这大约是9秒长的文件(numel(y) / fs = 8.9249秒).有了这个文件，这就是我得到的幅度和阶段:

I don't have access to your handel.wav file, but I'll be using the one provided with MATLAB. You can load this in with load handel;. The sampling frequency is stored in a variable called Fs, so I had to do fs = Fs; before the code I wrote above could work. The sampling frequency for this particular file is 8192 Hz, and this is approximately a 9 second long file (numel(y) / fs = 8.9249 seconds). With that file, this is the magnitude and phase that I get:

这篇关于在-fs/2至fs/2的范围内绘制WAV文件的幅度和相位谱的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！