Synthesize СДП в MatLab

Question

Я хочу, чтобы восстановить сигналы времени от заданной спектральной плотности мощности, предполагая нормальное распределение исходного сигнала:

PSD;         % [(m/s)^2/Hz] given spectrum 
T = 60;        % [s] length of original signal 
dt = 0.005;       % [s] time step of original signal 
N = T/dt;       % [-] number of samples 
NFFT = 2^nextpow2(N);    % [-] number of bins for FFT 
fs = 1/dt;       % [Hz] sampling frequency 
ASD = sqrt(PSD);      % [(m/s)/sqrt(Hz)] get amplitude spectrum 
omega = 2*pi*rand(NFFT/2,1);   % [rad] generate phase vector 
Z = ASD.*exp(1i*omega);    % create complex amplitude vector 
Z = [0;Z;flipud(conj(Z))];   % extend to satisfy symmetry 
Y = real(ifft(Z));     % inverse FFT 
[PSDY,f] = pwelch(Y,[],[],NFFT,fs); % generate PSD from Y to compare 

Результаты показывают спектр мощности на несколько порядков ниже, чем оригинал, но форма подходит очень хорошо. Я думаю, что что-то не так с единицами, или может отсутствовать коэффициент масштабирования. Я не уверен в единицах сигнала времени после ifft, так как амплитуда имеет [(m/s)/sqrt (Hz)].

Answer 1

Я считаю, что здесь есть две проблемы.
Во-первых, я думаю, что PSD, как вы его определяете (точнее, как вы его используете), находится в неправильных единицах.
При определении сигнала в качестве

Z = ASD.*exp(1i*omega); 

затем ASD должны быть в m/s и не в (m/s)/Hz.
Так что вы должны сделать что-то вроде этого:

ASD = sqrt(PSD*fs/2) 

Теперь, так как PSD в (m/s)^2/Hz, ASD в единицах m/s.

Далее следует нормализовать ifft. То есть, вы должны определить Y как

Y = ifft(Z)*sqrt(NFFT); 

еще одна вещь, я не уверен, если это нарочно, но следующая строка

[PSDY,f] = pwelch(Y,[],[],NFFT,fs); 

приводит Y разделена на 8 частей (с длиной <NFFT) с перекрытием 50%. Каждая часть равна нулю от длины до NFFT.
Лучше практикой было бы использовать что-то вроде

[PSDY,f] = pwelch(Y,L,L/2,L,fs); 

для некоторого L или

[PSDY,f] = pwelch(Y,NFFT,[],NFFT,fs); 

, если вы настаиваете. Чтобы узнать больше, перейдите на http://www.mathworks.com/help/signal/ref/pwelch.html

В заключение, это ваш (модифицированный) код:

PSD = 5;         % [(m/s)^2/Hz] given spectrum 
T = 60;        % [s] length of original signal 
dt = 0.005;       % [s] time step of original signal 
N = T/dt;       % [-] number of samples 
NFFT = 2^nextpow2(N);    % [-] number of bins for FFT 
fs = 1/dt;       % [Hz] sampling frequency 
ASD = sqrt(PSD*fs/2);      % [(m/s)] get amplitude spectrum 
omega = 2*pi*rand(NFFT/2,1);   % [rad] generate phase vector 
Z = ASD.*exp(1i*omega);    % create complex amplitude vector 
Z = [0;Z;flipud(conj(Z))];   % extend to satisfy symmetry 
Y = ifft(Z)*sqrt(NFFT);     % inverse FFT 
[PSDY,f] = pwelch(Y,256,128,256,fs); % generate PSD from Y to compare 

что приводит к

, где синяя линия оценка PSD ,