Спектр, вычисленный с помощью Matlab FFT, не дает последовательного результата для разных длин выборки (одинаковое количество точек, но Fs различны)?


Я хотел бы построить профиль rugosity (из измерений АСМ), но все еще есть то, что я неправильно понимаю относительно FFT (особенно в документации Matlab).

Я хочу сравнить два измерения, то есть два профиля ковровости. Они были сделаны на одной и той же поверхности, они расходятся только на том, что одна сделана на меньшей длине, чем другая. Однако для каждого профиля у меня есть одинаковое количество выборочных измерений (здесь N=512). Допустим, это мои профили, t10 и t100-это X-абсцисса, вдоль которой производится измерение, и d10 и d100 - вертикальная координата, т. е. высота измерения в профиле ковровости.

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 2

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 периода позже!

Построение d100 и d10 вместе и увеличение масштаба на первых 10 секундах сигнала выявляет проблему:

Введите описание изображения здесь

Следовательно, более низкая частотная составляющая оценивается правильно (широкий пик для d10 обусловлен несколькими периодами, а не целым числом их), но более высокочастотный компонент, алиасированный для d100, появляется на более низкой частоте.

Кстати: что касается изменения длины преобразования в степень два: это ускорит вычисления в целом,но в этом случае у вас уже есть мощность сигнала двух длин. Результаты не будут лучше, это просто скорость вычислений.