使用 Matlab FFT 计算的频谱对于不同长度的样本(相同数量的点但 Fs 不同)没有给出一致的结果?
Spectrum computed with Matlab FFT does not give a consistent result for different lengths of sample (same number of points but Fs different)?
我想绘制粗糙度曲线(来自 AFM 测量),但仍然存在我对 FFT 的误解(尤其是在 Matlab 文档中)。
我想比较两个测量值,a.k.a。两个粗糙度配置文件。它们是在同一表面上完成的,它们的不同之处仅在于一个的长度比另一个短。不过,对于每个配置文件,我都有相同数量的样本度量(此处 N=512)。假设这些是我的配置文件,t10
和 t100
是进行测量的 x 横坐标,d10
和 d100
是垂直坐标,a.k.粗糙度剖面中测量的高度。
N=512;
t10 = linspace(0,10, N);
t100= linspace(0,100, N);
d10 = sin(2*pi*0.23 .*t10)+cos(2*pi*12 .*t10);
d100 = sin(2*pi*0.23 .*t100)+cos(2*pi*12 .*t100);
由于我测量的是同一个表面,但空间分辨率不同,即采样周期不同,所以这些粗糙度剖面的单边振幅谱应该重叠,不是吗?
与我应该获得的不同,我有以下图表:
和
使用以下函数:
function [f,P1,S1] = FFT_PowerSpectrumDensity(time,signal,flagfig)
H=signal;
X=time;
ell=length(X);
L = ell;% 2^(nextpow2(ell)-1) % Next power of 2 from length of the signal
deltaTime = mean(diff(X));
Fs=1/deltaTime; %% mean sampling frequency
%% Compute the Fourier transform of the signal.
Y = fft(H);
%% Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L.
P2 = abs(Y/L); % abs(fft(signal Y)) / Length_of_signal
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
if flagfig~=0
figure(flagfig)
loglog(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)','FontSize',18)
xlabel('Spatial frequency f=1/\lambda (m^{-1})','FontSize',14)
ylabel('|P1(f)| (m)','FontSize',14)
end
S = (Y.*conj(Y)).*(2/L).^2; % power spectral density
S1 = S(1:L/2+1);
S1(2:end-1) = S1(2:end-1);
%% Power spectrum (amplitude = a^2+b^2), in length^2
if flagfig~=0
figure(flagfig+1)
loglog(f,S1)
title('Power spectrum','FontSize',18)
xlabel('Spatial frequency f=1/\lambda (m^{-1})','FontSize',14)
ylabel('(Y*2/L)^2 (m^2)','FontSize',14)
end
end
例如,我使用以下命令调用此函数:
[f10, S10]= FFT_PowerSpectrumDensity(t10, d10, 10);
我应该使用 L=2^pow2(ell)-1)
吗?我知道它为 FFT 函数提供了更好的输入?另外,我不太确定我应该找到的大多数单位和值。
感谢您的帮助、指正和建议。
你的问题出在你的输入信号上:
N=512;
t100 = linspace(0,100, N);
d100 = sin(2*pi*0.23 .*t100)+cos(2*pi*12 .*t100);
d100
欠采样。您的第一个样本中有 cos(0)
,然后第二个样本中有 cos(2*pi*12*0.1953
= cos(2*pi*2.3436)
。那是 2.3 个周期之后!
一起绘制 d100
和 d10
,并放大信号的前 10 秒,揭示了问题:
因此,低频分量被正确估计(d10
的宽峰是由于周期少而不是整数),但是高频分量,别名 d100
,出现频率较低。
顺便说一句:关于将变换长度更改为 2 的幂:这通常会加快计算速度,但在这种情况下,您已经拥有了 2 倍长度信号的幂。结果不会更好,只是计算速度。
我想绘制粗糙度曲线(来自 AFM 测量),但仍然存在我对 FFT 的误解(尤其是在 Matlab 文档中)。
我想比较两个测量值,a.k.a。两个粗糙度配置文件。它们是在同一表面上完成的,它们的不同之处仅在于一个的长度比另一个短。不过,对于每个配置文件,我都有相同数量的样本度量(此处 N=512)。假设这些是我的配置文件,t10
和 t100
是进行测量的 x 横坐标,d10
和 d100
是垂直坐标,a.k.粗糙度剖面中测量的高度。
N=512;
t10 = linspace(0,10, N);
t100= linspace(0,100, N);
d10 = sin(2*pi*0.23 .*t10)+cos(2*pi*12 .*t10);
d100 = sin(2*pi*0.23 .*t100)+cos(2*pi*12 .*t100);
由于我测量的是同一个表面,但空间分辨率不同,即采样周期不同,所以这些粗糙度剖面的单边振幅谱应该重叠,不是吗?
与我应该获得的不同,我有以下图表:
使用以下函数:
function [f,P1,S1] = FFT_PowerSpectrumDensity(time,signal,flagfig)
H=signal;
X=time;
ell=length(X);
L = ell;% 2^(nextpow2(ell)-1) % Next power of 2 from length of the signal
deltaTime = mean(diff(X));
Fs=1/deltaTime; %% mean sampling frequency
%% Compute the Fourier transform of the signal.
Y = fft(H);
%% Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L.
P2 = abs(Y/L); % abs(fft(signal Y)) / Length_of_signal
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
if flagfig~=0
figure(flagfig)
loglog(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)','FontSize',18)
xlabel('Spatial frequency f=1/\lambda (m^{-1})','FontSize',14)
ylabel('|P1(f)| (m)','FontSize',14)
end
S = (Y.*conj(Y)).*(2/L).^2; % power spectral density
S1 = S(1:L/2+1);
S1(2:end-1) = S1(2:end-1);
%% Power spectrum (amplitude = a^2+b^2), in length^2
if flagfig~=0
figure(flagfig+1)
loglog(f,S1)
title('Power spectrum','FontSize',18)
xlabel('Spatial frequency f=1/\lambda (m^{-1})','FontSize',14)
ylabel('(Y*2/L)^2 (m^2)','FontSize',14)
end
end
例如,我使用以下命令调用此函数:
[f10, S10]= FFT_PowerSpectrumDensity(t10, d10, 10);
我应该使用 L=2^pow2(ell)-1)
吗?我知道它为 FFT 函数提供了更好的输入?另外,我不太确定我应该找到的大多数单位和值。
感谢您的帮助、指正和建议。
你的问题出在你的输入信号上:
N=512;
t100 = linspace(0,100, N);
d100 = sin(2*pi*0.23 .*t100)+cos(2*pi*12 .*t100);
d100
欠采样。您的第一个样本中有 cos(0)
,然后第二个样本中有 cos(2*pi*12*0.1953
= cos(2*pi*2.3436)
。那是 2.3 个周期之后!
一起绘制 d100
和 d10
,并放大信号的前 10 秒,揭示了问题:
因此,低频分量被正确估计(d10
的宽峰是由于周期少而不是整数),但是高频分量,别名 d100
,出现频率较低。
顺便说一句:关于将变换长度更改为 2 的幂:这通常会加快计算速度,但在这种情况下,您已经拥有了 2 倍长度信号的幂。结果不会更好,只是计算速度。