在 -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 范围内重复此操作。
%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 是信号的采样频率，因此 fs/2 是奈奎斯特频率。请注意，仅当信号是真实的时，考虑半频谱才有效。这意味着负光谱与正光谱对称，因此您实际上不需要在这里考虑负光谱。
但是，您想绘制幅度和相位的完整频谱。请注意，在使用 MATLAB 计算 fft 时，它使用 Cooley-Tukey 算法，因此在计算 N 点 FFT 时，结果的一半是针对从 0 Hz 到 fs/2 Hz 独占，另一半用于从 -fs/2 Hz 到 0 Hz 独占的频率。
因此，要绘制完整频谱，只需对完整信号执行 fftshift，以便交换频谱的右半部分和左半部分，从而使 0 Hz 频率位于中心的信号。此外，您必须生成 -fs/2 到 fs/2 之间的频率以覆盖整个频谱。具体来说，您需要生成 N 个线性间隔在 -fs/2 到 fs/2 之间的点。但是，请注意，fs/2 Hz 的奈奎斯特频率在最后被排除在外，因此您需要在 -fs/2 到 fs/2 之间生成 N+1 点并删除最后一个点为了使每个频率仓之间的正确步长正确。生成此线性点阵列的最简单方法是使用 linspace 命令，其中起始频率为 -fs/2，结束频率为 fs/2，并且您希望 N+1 点之间这个范围并删除最后一点：
freq = linspace(-fs/2, fs/2, n+1);
freq(end) = [];
因此，借用您的代码的某些部分，这就是修改后的代码绘制幅度和相位的完整频谱的样子：
%// 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 秒）。使用该文件，这是我得到的幅度和相位：
对于离散傅立叶变换 (DFT) 及其快速实现 (FFT)，频率使用采样频率 fs 进行归一化，即原始范围 -fs/2:fs/2更改为 -pi:pi.
此外，DFT/FFT总是从0开始，你可以使用fftshift()将0频率移到中心。因此，在fftshift()之后，范围是-pi:pi，那么，你可以缩放到-fs/2:fs/2.
看看下面的Matlab函数，它可以非常精确地计算相位谱和幅度谱：
该程序以可接受的精度计算输入信号的幅度和相位谱，尤其是在相位计算方面 spectrum.The 代码在计算幅度和相位谱方面执行三个主要工作。首先，它将输入信号扩展到无穷大；因为为了计算傅立叶变换 (FT)（Matlab 中的 fft 函数），我们认为我们的信号是周期性的，具有无限波长，代码通过将原始信号放在自身旁边直到 [=31= 的长度来创建 super_signal ] 大约是1000000个样本，为什么我选择了1000000个样本？实际上，它只是基于尝试和错误！对于我尝试过的大多数信号，具有 1000000 个样本的超级信号具有最佳输出。
其次，在Matlab中计算fft可以选择不同的分辨率，Mathwork文档和帮助使用NFFT=2^nextpow2(length(signal))，这对于想要高精度输出的人来说肯定是不够的。在这里，我选择适用于大多数信号的分辨率 NFFT=100000。
三、代码对FT的结果进行阈值过滤，非常重要的一步！对于相位谱的计算，其结果由于浮动舍入误差而非常嘈杂，导致在计算过程中 "arctan" 即使是很小的舍入误差也会在相位谱的结果中产生显着的噪声，为了抑制这种噪声，您可以定义一个阈值。这意味着如果特定频率的幅度小于预定义的阈值（您必须定义它），它会用零代替它。
这三个步骤有助于显着改善幅度和相位谱的结果。
如果您在研究中使用此程序，请引用以下论文：
Afshin Aghayan、Priyank Jaiswal 和 Hamid Reza Siahkoohi (2016)。 "Seismic denoising using the redundant lifting scheme." 地球物理学，81(3)，V249-V260。 https://doi.org/10.1190/geo2015-0601.1
我在绘制 wav 文件的 FFT 时遇到问题。我设法绘制了信号的幅度和相位频谱，但是我需要在 -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 是信号的采样频率，因此 fs/2 是奈奎斯特频率。请注意，仅当信号是真实的时，考虑半频谱才有效。这意味着负光谱与正光谱对称，因此您实际上不需要在这里考虑负光谱。
但是，您想绘制幅度和相位的完整频谱。请注意，在使用 MATLAB 计算 fft 时，它使用 Cooley-Tukey 算法，因此在计算 N 点 FFT 时，结果的一半是针对从 0 Hz 到 fs/2 Hz 独占，另一半用于从 -fs/2 Hz 到 0 Hz 独占的频率。
因此，要绘制完整频谱，只需对完整信号执行 fftshift，以便交换频谱的右半部分和左半部分，从而使 0 Hz 频率位于中心的信号。此外，您必须生成 -fs/2 到 fs/2 之间的频率以覆盖整个频谱。具体来说，您需要生成 N 个线性间隔在 -fs/2 到 fs/2 之间的点。但是，请注意，fs/2 Hz 的奈奎斯特频率在最后被排除在外，因此您需要在 -fs/2 到 fs/2 之间生成 N+1 点并删除最后一个点为了使每个频率仓之间的正确步长正确。生成此线性点阵列的最简单方法是使用 linspace 命令，其中起始频率为 -fs/2，结束频率为 fs/2，并且您希望 N+1 点之间这个范围并删除最后一点：
freq = linspace(-fs/2, fs/2, n+1);
freq(end) = [];
因此，借用您的代码的某些部分，这就是修改后的代码绘制幅度和相位的完整频谱的样子：
%// 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 秒）。使用该文件，这是我得到的幅度和相位：
对于离散傅立叶变换 (DFT) 及其快速实现 (FFT)，频率使用采样频率 fs 进行归一化，即原始范围 -fs/2:fs/2更改为 -pi:pi.
此外，DFT/FFT总是从0开始，你可以使用fftshift()将0频率移到中心。因此，在fftshift()之后，范围是-pi:pi，那么，你可以缩放到-fs/2:fs/2.
看看下面的Matlab函数，它可以非常精确地计算相位谱和幅度谱：
该程序以可接受的精度计算输入信号的幅度和相位谱，尤其是在相位计算方面 spectrum.The 代码在计算幅度和相位谱方面执行三个主要工作。首先，它将输入信号扩展到无穷大；因为为了计算傅立叶变换 (FT)（Matlab 中的 fft 函数），我们认为我们的信号是周期性的，具有无限波长，代码通过将原始信号放在自身旁边直到 [=31= 的长度来创建 super_signal ] 大约是1000000个样本，为什么我选择了1000000个样本？实际上，它只是基于尝试和错误！对于我尝试过的大多数信号，具有 1000000 个样本的超级信号具有最佳输出。
其次，在Matlab中计算fft可以选择不同的分辨率，Mathwork文档和帮助使用NFFT=2^nextpow2(length(signal))，这对于想要高精度输出的人来说肯定是不够的。在这里，我选择适用于大多数信号的分辨率 NFFT=100000。
三、代码对FT的结果进行阈值过滤，非常重要的一步！对于相位谱的计算，其结果由于浮动舍入误差而非常嘈杂，导致在计算过程中 "arctan" 即使是很小的舍入误差也会在相位谱的结果中产生显着的噪声，为了抑制这种噪声，您可以定义一个阈值。这意味着如果特定频率的幅度小于预定义的阈值（您必须定义它），它会用零代替它。
这三个步骤有助于显着改善幅度和相位谱的结果。
如果您在研究中使用此程序，请引用以下论文：
Afshin Aghayan、Priyank Jaiswal 和 Hamid Reza Siahkoohi (2016)。 "Seismic denoising using the redundant lifting scheme." 地球物理学，81(3)，V249-V260。 https://doi.org/10.1190/geo2015-0601.1