0

我想从给定的功率谱密度恢复时间信号，假设原始信号的正态分布：

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)]。

4

1 回答 1

1

I believe there are two problems here.
First, I think that the PSD as you define it (or rather, as you use it) is in the wrong units.
When you define the signal as 

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

then ASD should be in m/s and not in (m/s)/Hz.
So you should do something like that: 

ASD = sqrt(PSD*fs/2)

Now, since PSD is in (m/s)^2/Hz, ASD is in units of m/s.

Next, the ifft should be normalised. That is, you should define the Y as

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

One more thing, I am not sure if this is on purpose, but the following line

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

results in Y being divided into 8 parts (with length <NFFT) with 50% overlap. Each part is zero padded to length of NFFT.
A better practice would be to use something like

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

for some L or 

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

if you insist. To find out more, go to http://www.mathworks.com/help/signal/ref/pwelch.html

In conclusion, this is your (modified) code:

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

which results in

enter image description here

where the blue line is the estimate PSD.

于 2014-09-11T19:56:16.027 回答