Поиск спектра мощности сигнала с использованием двух подходов

Question

Я пытаюсь найти спектр мощности сигнала. Длина сигнала 100000, частота выборки 1000Hz, а количество баллов - 100000. Я нашел спектр мощности, используя два подхода. Первый состоит в том, чтобы взять всю длину как одну часть и найти спектр мощности для нее, а второй подход - делить сигнал на 100*1000 и найти спектр для каждой строки, затем получить среднее значение для всех строк. Моя проблема в том, что я должен получить тот же ответ в обоих подходах, но у меня есть разные ответы. Я не знаю, в чем ошибка в моем коде.

N=100000; 
SF=1000;  
a=0.1; 
b=0.3; 
amplitude1=1; 
amplitude2=0.5; 
t=0:1/SF:100; 
f1=SF*a; 
f2=SF*b; 
A=amplitude1*sin(2*pi*f1*t)+amplitude2*sin(2*pi*f2*t); 
Y=2*randn(1,length(A))+A; 
bin=[0 :N/2]; 
fax_Hz=(bin*SF)/N; 
FFT=fft(Y); 
spectra=2/(SF*length(Y))*(FFT.*conj(FFT)); 
plot(fax_Hz,spectra(1,1:50001)); 
D=reshape(Y(1,1:100000),[100,1000]); 
M=length(D(1,:)); 
for i=1:100 
    FFT_1(i,:)=fft(D(i,:)); 
    S(i,:)=(2/(SF*M))*(FFT_1(i,:).*conj(FFT_1(i,:))); 
end 
S_f=mean(S); 
figure 
plot (S_f); 

Я просто обновляю код. Я не знаю, но когда я добавил шум, чтобы сигнализировать, что два сюжета выглядят смещенными.

Answer 1

Основная проблема: reshape Вы работаете с каждой строкой, являющейся отдельной последовательностью. Однако Reshape заполняет первый столбец, прежде чем переходить ко второму.

Вместо этого вы можете использовать следующее.

D=reshape(A(1,1:100000),[1000,100]).'; 

Нормализация - еще одна проблема. Вы можете использовать ifft вместо fft, поскольку он нормализуется по умолчанию (не уверен, почему). Или, альтернативно, сохраните свою нормализацию и вместо использования mean вы можете использовать sum, возможно, это связано с ошибкой, которую вы, возможно, сделали. По-прежнему существует небольшое расхождение в амплитудах, не уверен, откуда это происходит.

В конце концов участок использовать следующее:

bin=[0 :N]; 
fax_Hz=(bin*SF)/N; 
FFT=ifft(A); 
spectra=FFT.*conj(FFT); 
plot(fax_Hz,spectra); hold on 
D=reshape(A(1,1:100000),[1000,100]).'; 
M=length(D(1,:)); 
for i=1:100 
    FFT_1(i,:)=ifft(D(i,:)); 
    S(i,:)=FFT_1(i,:).*conj(FFT_1(i,:)); 
end 
S_f=mean(S); 
plot(fax_Hz(1:100:end-1), S_f); 

Замечание: при fax_Hz(1:100:end-1) является Hacky способ получить длину векторов быть одинаковыми.