nfft

nfft

关注
发信
关注(28)粉丝(399)

matlab - 用Matlab进行傅立叶变换

我的问题类似于this post，但比answers更笼统，我认为在归一化方面存在错误，无论如何，最新版本的Matlab(2015)仍然存在。如果您认为更合适的话，请在评论中告诉我，我犹豫将其发布在CodeReview SE上。

我想验证使用Matlab的fft进行傅立叶变换的以下代码，因为我发现网络上的信息源相互矛盾，包括在Matlab帮助本身中，并且我无法通过某些“配方”(包括来自MathWorks团队的defined，请参见下文)，尤其是那些提取单面光谱以用于实际输入的配方。

例如，通常在网上找到用于在提取正频率时解决实值信号对称频谱的加倍幅度似乎是错误的(Parseval定理失败)，取而代之的是必须使用Matlab中的两个系数(我不知道为什么)。有些人似乎也像Y = fft(X)/L一样直接对DFT系数进行归一化，但是我认为这很令人困惑，因此不建议使用。振幅为own example，因为复数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的posted和clarification post上:
>> 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团队自己的answer:
>> 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上以ojit_a的形式找到:
>> 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 )

最佳答案

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。同样，每个wave的功率是|Xk|^2 / N。 Matlab也使用该归一化(FFTW确实如此)。
  • 对于real-valued inputs，DFT系数是对应的正/负频率，除ott_rstrong DC分量(平均项，频率为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中的代码是不正确的，取而代之的是必须同时输出幅度和功率，因为​​没有通用的归一化系数可以容纳复杂或真实的情况，而偶数或奇数个时间点。您可以找到要点here。

    关于matlab - 用Matlab进行傅立叶变换，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/35177274/

    10-11 20:10
    Powered by LMLPHP ©2024 nfft 0.033242