用 Matlab 进行傅里叶变换

Question

我的问题与this post 类似，但比this post 更笼统，而且我认为在标准化方面存在错误，无论如何都是最新版本的Matlab（2015）。我犹豫是否将其发布在 CodeReview SE 上，如果您认为在 cmets 中告诉我更合适的话。

我想使用 Matlab 的 fft 验证傅立叶变换的以下代码，因为我在网络上发现了相互冲突的信息来源，包括在 Matlab 帮助本身中，我有无法使用某些此类“配方”验证 Parseval 定理（包括来自 MathWorks 团队的answers，见下文），尤其是那些为真实输入提取单面光谱的。

例如，仅在提取正频率时，通常在网上发现的用于解释实值信号的对称频谱的幅度加倍似乎是错误的（Parseval 定理失败），相反，似乎有必要使用平方- Matlab中两个系数的根（我不知道为什么）。有些人似乎也像Y = fft(X)/L 一样直接对DFT 系数进行归一化，但我认为这很混乱，应该不鼓励；幅度是defined，作为复数DFT系数除以信号长度的模，系数本身不应该被除。一旦通过验证，我打算将此代码作为要点发布在 GitHub 上。

function [frq,amp,phi] = fourier_transform( time, vals )
% FOURIER_TRANSFORM computes the Fast Fourier Transform of a given time-series.
% 
% [freq,amp,phi] = fourier_transform(time,vals)
%
% Inputs:
%
%   time  - Nx1 column vector of equally-spaced timepoints (if not, input will be resampled).
%   vals  - NxM matrix of real- or complex-valued observations at previous timepoints.
%
% Outputs:
%
%   frq   - column vector of frequencies in Hz
%   amp   - corresponding matrix of amplitudes for each frequency (row) and signal (column)
%   phi   - corresponding unwrapped phase for each frequency (row) and signal (column)
%
% Note:
%   To compute the power from the output amplitude, you need to multiply by the number of timepoints:
%       power = numel(time) * amp.^2;
%
% References:
%   https://en.wikipedia.org/wiki/Discrete_Fourier_transform

    % make sure input time-series is uniformly sampled
    assert( iscolumn(time), 'Input time should be a column vector.' );
    assert( ismatrix(vals) && size(vals,1) == numel(time), 'Input values should be a matrix with same number of rows than of timepoints.' );

    if std(diff(time)) > 1e-6
        warning('Input time-course does not appear to be uniformly sampled. Resampling before continuing.');
        [vals,time] = resample(vals,time);
    end

    % sampling information
    nt = numel(time);
    dt = time(2)-time(1);
    fs = 1/dt;
    df = fs/nt;

    % complex spectrum coefficients
    coef = fft(vals);

    % real input
    if isreal(vals)

        % extract one-sided spectrum (spectrum is symmetric)
        nfft = floor( nt/2 + 1 ); % eg 8 -> 5, and 7 -> 4
        coef = coef( 1:nfft, : );
        frq  = (0:nfft-1)*df;

        % correct amplitude values to account for previous extraction
        fac = sqrt(2);
        amp = fac*abs(coef)/nt;
        amp(1,:) = amp(1,:)/fac; % .. except for the DC component
        if mod(nt,2) == 0
            amp(end,:) = amp(end,:)/fac; % .. and for the Nyquist frequency (only if nt is even)
        end

    % complex input
    else

        % shift the spectrum to center frequencies around 0
        coef = fftshift( coef );
        frq  = fftshift( (0:nt-1)*df );
        amp  = abs(coef)/nt;

    end

    % make sure frq is a column vector and compute phases
    frq = frq(:);
    phi = unwrap(angle(coef));

end

---

使用示例

>> fs=1e3; t=transpose(0:1/fs:10); nt=numel(t); X=rand(nt,1);
>> [frq,amp,phi] = fourier_transform( t, X ); 
>> sum( abs(X).^2 ) - nt*sum( amp.^2 ) % Parseval's theorem

ans =

  -2.7285e-11

错误示例 1

来自 Matlab 的 own example 和 SO 上的 posted：

>> Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));     % Sinusoids plus noise
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
>> sum(abs(y).^2) - NFFT*sum(abs(Y).^2) % Parseval's theorem

ans =

 -220.4804

问题及解决方案：

来自Y = fft(y,NFFT)/L 行中的规范化。 这应该是：

>> Y = fft(y,NFFT);
Ya = abs(Y)/NFFT; % correctly normalised amplitudes
sum(abs(y).^2) - NFFT*sum(Ya.^2) % Parseval's theorem

错误示例 2

来自 MathWorks 团队自己的clarification post：

>> Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));     % Sinusoids plus noise

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

NumUniquePts = ceil((NFFT+1)/2);
Y=Y(1:NumUniquePts);
MX=2*abs(Y);
MX(1)=MX(1)/2; % DC component
if ~rem(NFFT,2)  % when NFFT is even, Y(1+Y/2) is the Nyquist frequency component of x, and needs to make sure it is unique.
    MX(length(MX))=MX(length(MX))/2;
end

>> sum( abs(y).^2 ) - NFFT*sum( MX.^2 )

ans =

  -5.3812e+03

问题及解决方案：

再次标准化。将Y = fft(y,NFFT)/L; 替换为Y = fft(y,NFFT)，并将MX=2*abs(Y); 替换为MX=2*abs(Y)/NFFT;。但是这里出现了倍频问题；校正因子似乎是 sqrt(2) 而不是 2。

错误示例 3

在 MatlabCentral 上以 answer 的形式找到：

>> Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
x = 0.7*sin(2*pi*Fs/8*t) + sin(2*pi*Fs/4*t); 
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(x,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
>> sum( abs(x).^2 ) - NFFT*sum( abs(Y).^2 )

ans =

  -36.1891

问题及解决方案：

和第一个例子一样，归一化问题。改为写：

Y  = fft(x,NFFT);
Ya = abs(Y)/NFFT;
sum( abs(x).^2 ) - NFFT*sum( abs(Ya).^2 )

Answer 1

TL;DR（摘要）

很难找到使用 Matlab 正确标准化幅度/功率值的 fft 使用的在线示例（例如，可以使用 Parseval's theorem 进行验证）。如果您想比较不同长度的信号之间的光谱，这一点至关重要。实值信号还有一个额外的问题，因为在这种情况下，频谱通常仅针对正频率计算，因此需要缩放幅度或功率值以考虑频率折叠。 根据下面的帖子和答案，here is a gist 我认为它可以正确且一致地为实值和复值输入调整系数。

带回家的信息是：

• 不要直接对 DFT 系数进行归一化（例如不要写 Y = fft(x)/L）；
• 如果您使用 n 点变换（例如 fft(x,nfft)），则归一化器是 nfft 而不是 numel(x)；
• 如果提取单面频谱，则需要根据共轭对 DFT 系数对应的幅度/功率值进行调整；
• 如果您提取单面频谱，则应分别计算幅度和功率（即不要根据幅度计算功率）。

---

幅度、功率和单面光谱

定义和解释on Wikipedia:

• DFT 系数是复数且未归一化，而逆 DFT 的公式在总和前面带有 1/N 因子。这在某种意义上是自然的，因为在时间-频率方向上的移动可以看作是在不同频率的（正交）波的基础上的投影，而在频率-时间方向上的移动可以看作是加权叠加波浪。
• 在那个叠加中，整体震级应该是原始时间点的震级（即，它是一个反演），而加权平均中每个波的 幅值就是对应的 DFT 系数除以波数|Xk| / N。同样，每一波的power为|Xk|^2 / N。 Matlab 也使用这种归一化（嗯，FFTW 确实如此）。
• 对于real-valued inputs，DFT 系数是对应正/负频率的共轭对，除了直流分量（平均值，频率 0），以及当时间点是均匀的。在实践中，大多数应用程序通过仅提取正频率的 DFT 系数来避免这种冗余。这会导致幅度和功率的离散值变得复杂。
• 对应于成对 DFT 系数的幅度（除了第一个和奈奎斯特频率存在时除外）可以简单地加倍，然后丢弃负频率。功率也一样。
• 与功率类似，但请注意，在这种情况下使用调整后的幅度值计算离散功率值是不正确的。实际上 N * amp_adjusted[k]^2 = N * (2*|Xk|/N)^2 不等于 2*|Xk|^2 / N （这是 OP 中两个平方根的来源）。因此，有必要独立地从 DFT 系数中计算幅度和功率值（另一个很好的理由不对其进行缩放）。

N 点变换

许多在线示例使用显式 N 点变换：Y = fft(x,NFFT) 其中NFFT 通常是 2 的幂，从而使 FFTW 的计算效率更高。

这种情况下的有效区别（假设 NFFT >= N）是 x 在其末尾填充 0 直到它达到 NFFT 时间点的长度。这意味着分解中的频率数量发生了变化，因此应该相对于NFFT 波分量进行归一化，而不是原始的N 时间点。

因此，在网上找到的几乎所有示例在标准化系数的方式上都是错误的。它不应该是Y = fft(x,NFFT)/N，而应该是Y = fft(x,NFFT)/NFFT——这又是一个放松对复数系数进行归一化习惯的好理由。

请注意，这对 Parseval 的相等性没有影响，因为时域中的相加项都为零，因此它们对现在更大的和的贡献也为零。但在频域中，添加的离散频率通常会对原始信号产生响应，这直观地说明了为什么在填充和未填充变换之间获得的系数确实可能完全不同。

代码

因此，OP 中的代码不正确，相反，似乎有必要同时输出幅度和功率，因为​​没有通用的归一化系数可以适应复杂和真实的情况，包括偶数或奇数时间-点。您可以找到 Gist here。

Answer 2

我发现这个在线 fft 示例 (https://habr.com/post/112068/) 非常有用。看看这个：

%% Parameters
Tm=5;% Length of signal (s)
Fd=512;% Sampling frequency (Hz)
Ak=0.5;% Constant component (Unit)
A1=1;% The amplitude of the first sinusoid (Unit)
A2=0.7;% Amplitude of the second sinusoid (Unit)
F1=13;% Frequency of the first sinusoid (Hz)
F2=42;% Frequency of the second sinusoid (Hz)
Phi1=0;% Initial phase of the first sinusoid (Degrees)
Phi2=37;% The initial phase of the second sinusoid (Degrees)
An=3*A1;% Noise Dispersion (Unit)
FftL=1024;% Number of Fourier Spectrum Lines
%% Generating work arrays
T=0:1/Fd:Tm;% Time Arrays
Noise=An*randn(1,length(T));% An array of random noise with a length equal to the array of time
Signal=Ak+A1*sind((F1*360).*T+Phi1)+A2*sind((F2*360).*T+Phi2);% Signal array (a mixture of 2x sinusoids and a constant component)
%% Spectral representation of the signal
FftS=abs(fft(Signal,FftL));% The amplitudes of the Fourier transform of the signal
FftS=2*FftS./FftL;% Spectrum normalization by amplitude
FftS(1)=FftS(1)/2;% The normalization of the constant component in the spectrum
FftSh=abs(fft(Signal+Noise,FftL));% The amplitudes of the Fourier transform of the signal + noise mixture
FftSh=2*FftSh./FftL;% Spectrum normalization by amplitude
FftSh(1)=FftSh(1)/2;% The normalization of the constant component in the spectrum
%% Plotting
subplot(2,1,1);
plot(T,Signal);
title('Signal');
xlabel('Time (s)');
ylabel('Amplitude (Unit)');
subplot(2,1,2);
plot(T,Signal+Noise);
title('Signal+Noise');
xlabel('Time (s)');
ylabel('Amplitude (Unit)');

F=0:Fd/FftL:Fd/2-1/FftL;% The frequency array of the calculated Fourier spectrum
figure();
subplot(2,1,1);
plot(F,FftS(1:length(F)));%  Plotting of the spectrum of the Fourier signal
title('Signal spectrum');
xlabel('Frequency (Hz)');
ylabel('Amplitude (Unit)');
subplot(2,1,2);
plot(F,FftSh(1:length(F)));% Plotting of the Fourier signal spectrum
title('Signal spectrum');
xlabel('Frequency (Hz)');
ylabel('Amplitude (Unit)');