matlab - Matlab 傅里叶变换

Question

我的问题类似于，但比 this post 更笼统，而且我认为在规范化方面存在错误，无论如何使用最新版本的 Matlab(2015)。我犹豫是否将其发布到 CodeReview SE 上，如果您认为更合适，请在评论中告诉我。

我要验证以下代码 使用 Matlab 的 fft 的傅立叶变换，因为我在网络上发现了相互矛盾的信息来源，包括在 Matlab 帮助本身中，并且我无法用某些此类“食谱”验证 Parseval 定理(包括来自 MathWorks 团队的 answers，见下文) ，尤其是那些为实际输入提取单边光谱的人。

例如，在提取正频率时，通常在网上发现的用于解释实值信号的对称谱的倍幅似乎是错误的(Parseval 定理失败)，而似乎有必要使用Matlab 中的两个系数(我不知道为什么)。有些人似乎也直接对 DFT 系数进行归一化，如 Y = fft(X)/L ，但我认为这是令人困惑的，应该劝阻；振幅为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 , 和 posted关于 SO:

>> 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

发现为 answer在 MatlabCentral 上:

>> 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(摘要)

网上很难找到fft的例子。与 Matlab 一起使用以正确归一化幅度/功率值(例如，可以使用 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 .同样，电源每波是|Xk|^2 / N . Matlab 也使用归一化(好吧，FFTW 确实如此)。

对于 real-valued inputs ，DFT 系数是对应正/负频率的共轭对，除了 DC 分量(平均项，频率 0)，以及当时间点数为偶数时的 Nyquist 频率。实际上，大多数应用程序通过仅提取正频率的 DFT 系数来避免这种冗余。这导致幅度和功率的离散值的复杂化。

对应于成对 DFT 系数的幅度(除了第一个和存在的奈奎斯特频率之外的所有幅度)可以简单地加倍，然后丢弃负频率。权力也是一样。

功率类似，但请注意它是 不正确 在这种情况下，使用调整后的幅度值计算离散功率值。确实N * amp_adjusted[k]^2 = N * (2*|Xk|/N)^2不等于 2*|Xk|^2 / N (这是 OP 中 2 的平方根的来源)。因此需要计算独立 来自 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 中的代码是不正确的，似乎有必要同时输出幅度和功率，因为​​没有通用的归一化系数可以适应具有偶数或奇数时间点的复杂和真实情况。您可以找到要点 here .