很难找到fftMatlab使用的在线示例,正确地规范化幅度/功率值.如果您要在不同长度的不同信号上比较这些值,这一点至关重要.这通常是实值输入的问题,因为在这种情况下通常会提取单侧频谱,因此在计算幅度或功率值时应手动应用幅度变化.你可以在这里找到关于GitHub的要点(请告诉我任何错误).
带回家的消息是:
Y = fft(X)/L);Y = fft(y,NFFT)/L),那么规范器是Y = fft(y,NFFT)/L;和不是Y = fft(y,NFFT),我的问题类似于,但比这篇文章更为通用,我认为关于规范化存在错误,最新版本的Matlab(2015)无论如何.我对在CodeReview SE上发布此内容犹豫不决,如果您认为更合适,请在评论中告诉我.
我想用Matlab 来验证以下傅立叶变换代码MX=2*abs(Y);,因为我在网上找到了相互矛盾的信息来源,包括Matlab帮助本身,我无法用某些这样的"食谱"来验证Parseval定理(包括来自MathWorks团队的答案,见下文),特别是那些提取实际输入的单面光谱的答案.
例如,在提取正频率时,通常在网上发现的用于解释实值信号的对称频谱的幅度加倍似乎是错误的(Parseval定理失败),而似乎有必要使用平方根Matlab中的两个系数(我不知道为什么).有些人似乎也直接将DFT系数标准化MX=2*abs(Y)/NFFT;,但我认为这是令人困惑的,应该不鼓励; 幅度定义为复数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
>> 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
问题和解决方案:
来自线路的规范化sqrt(2).这应该是:
>> Y = fft(y,NFFT);
Ya = abs(Y)/NFFT; % correctly normalised amplitudes
sum(abs(y).^2) - NFFT*sum(Ya.^2) % Parseval's theorem
来自MathWorks团队自己的澄清帖子:
>> 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
问题和解决方案:
再次标准化.替换2为fft,并且据说Y = fft(X)/L由Y = fft(y,NFFT)/L.但是这里出现了幅度倍增问题; 修正系数似乎是,Y = fft(y,NFFT)/L;而不是Y = fft(y,NFFT).
在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 )
很难找到fft使用 Matlab 正确标准化幅度/功率值的在线示例（例如可以通过Parseval 定理进行验证）。如果您想比较不同长度信号之间的频谱，这一点至关重要。实值信号还存在一个额外的问题，因为在这种情况下，通常仅针对正频率计算频谱，因此需要缩放幅度或功率值以考虑频率折叠。在下面的帖子和答案之后，我认为这里的要点是针对实值和复值输入正确且一致地缩放系数。
带回家的信息是：
Y = fft(x)/L）；fft(x,nfft)），则归一化器是nfft，而不是numel(x)；正如维基百科上的定义和解释：
1/N在和前面带有一个因子。这在某种意义上是自然的，因为在时间到频率方向上的移动可以看作是在不同频率的（正交）波的基础上的投影，而在频率到时间方向上的移动可以看作是加权叠加的波浪。|Xk| / N。同样，每个波的功率|Xk|^2 / N为。Matlab 也使用这种标准化（嗯，FFTW 就是这样做的）。N * amp_adjusted[k]^2 = N * (2*|Xk|/N)^2不等于2*|Xk|^2 / N（这是OP中二的平方根的来源）。因此，有必要根据 DFT 系数独立计算幅度和功率值（不缩放它们的另一个充分理由）。许多在线示例都使用显式 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——这是放弃标准化复系数的习惯的另一个好理由。
请注意，这与帕塞瓦尔等式没有区别，因为时域中添加的项全部为零，因此它们对现在更大的总和的贡献也为零。但在频域中，添加的离散频率通常会对原始信号产生响应，这直观地说明了为什么在填充和未填充变换之间所获得的系数实际上可能有很大不同。
因此，OP中的代码是不正确的，相反，似乎有必要输出幅度和功率，因为没有通用的归一化系数可以适应具有偶数或奇数时间点的复杂和真实情况。您可以在这里找到要点。