使用 Matlab FFT 计算的频谱对于不同长度的样本(点数相同但 Fs 不同)给出的结果不一致?

Posted

技术标签:

【中文标题】使用 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)? 【发布时间】:2018-07-14 22:52:12 【问题描述】:

我想绘制粗糙度的轮廓(来自 AFM 测量),但我仍然对 FFT 有误解(尤其是在 Matlab 文档中)。

我想比较两个测量值,也就是两个粗糙度曲线。它们是在同一个表面上完成的,它们只是在一个比另一个的长度短这一事实上有所不同。不过,对于每个配置文件,我都有相同数量的样本度量(这里 N=512)。假设这些是我的配置文件,t10t100 是进行测量的 x 横坐标,d10d100 是垂直坐标,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 函数提供了更好的输入?另外,我不太确定我应该找到的大多数单位和值。

感谢您的帮助、更正和建议。

【问题讨论】:

【参考方案1】:

你的问题在于你的输入信号:

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 个周期之后!

d100d10 一起绘制,并放大信号的前 10 秒,可以发现问题:

因此,正确估计了较低频率分量(d10 的宽峰值是由于几个周期而不是它们的整数),但较高频率分量(别名为 d100)出现在较低的频率。

顺便说一句:关于将变换长度更改为 2 的幂:这通常会加快计算速度,但在这种情况下,您已经有了 2 的幂长度信号。结果不会更好,只是计算速度。

【讨论】:

好的,谢谢。这是有道理的,我增加了采样,幅度谱中的峰值似乎确实向不变值收敛......顺便说一下,你会对这段代码中计算的功率谱的计算有一些见解吗?我想我想计算功率密度谱,显然(例如根据dsp.stackexchange.com/questions/33957/…),似乎是别的东西......?非常感谢您的帮助和见解。 有趣的是,您链接的帖子又链接到***,这似乎说没有区别:“功率谱密度(或简称功率谱)”。 en.wikipedia.org/wiki/Spectral_density -- 我当然不知道有任何这样的区别。我认为计算将是S = (Y.*conj(Y))/L,但除以L 是微不足道的,无论您使用L 还是L^2 都符合惯例。毕竟,有趣的是函数的形状。但老实说,我已经很多年没有计算过功率谱了,所以我有点生疏了。 :) 是的,这很奇怪,在我链接的评论中,据我了解,这家伙说有区别(“在第一种情况下,每个分立元件都有功率单位(W,mW等)。在第二种情况下,连续频谱中的每个点都有每个频率的功率单位(W/Hz、mW/Hz 等)。在第二种情况下,必须在一个频带上积分以获得功率单位。”)无论如何,它似乎代表的是同一件事,并且不会改变图形的形状,正如您所说,这是我们感兴趣的。 我主要想知道研究表面皱纹的人从图表的角度来看是什么,他们从他们的形状中得到了什么......? @LeChat:我绝对不是问这个问题的人,但我猜峰值是相关的,因为它们表明了一些重要的周期性;高频功率表明皱纹的精细程度;非常低频率的功率指示表面是否是直的。或类似的东西。我不知道你的领域! :)

以上是关于使用 Matlab FFT 计算的频谱对于不同长度的样本(点数相同但 Fs 不同)给出的结果不一致?的主要内容,如果未能解决你的问题,请参考以下文章

对12bit的AD建模后matlab写程序对其进行FFT,计算SNR并输出频谱图,程序该怎么写???急

信号处理常用matlab之频谱

MATLAB画信号频谱的子函数

matlab怎样进行频谱分析

利用matlab怎样进行频谱分析

matlab 作出信号频谱图