请教个白噪声合成的程序问题，望高手们赐教！

wyuxuenvni

最近在学习功率谱合成，编程过程中发现了一些问题，诚恳高手们赐教，小女不胜感激！

问题：
（1）为什么最后要求两次功率谱，第一次  p_sd = pwelch ( y , [] , [] , nfft , 10.24);     第二次[ p_sd_t , ff ] = pwelch ( y , [] , [] , 256 , 10.24 ); 做第二幅图时采取时采取plot ( ff , p_sd_t , 'b' , w , S , 'r' )，而非 plot ( w , p_sd , 'b' , w , S , 'r' )呢？

（2）合成功率谱时必须对傅里叶幅值谱进行修正吗？

为说明白问题，先将原程序贴出来：


%按照功率谱生成平稳过滤白噪声时程数据
clear
clc

% 设定目标功率谱相关参数

ksg = 0.63;           % 功率谱衰减系数，已知
wg = 2.5;              %场地土的卓越频率，已知
S0 = 0.002;            %反应地震动强弱程度的谱参数，已知
w = 0 : 0.005 : 5.12;     % 定义离散频率向量   【 w=0 ：fs/nfft ： （nfft/2-1）*fs/nfft 】

% 设定时程曲线的相关参数

fs = 10.24;            %目标函数的采样频率
tl = 100;               % 目标函数的总采样时间
iter_num = 20;            %迭代次数

% 计算目标功率谱

[ h , l ] = size(w);      % 计算目标功率谱的频率点数
S = [];                   % 定义目标功率谱向量
for ii = 1 : l
    SS = ( 1 + 4 * ksg ^2 * w(ii)^2 / wg^2 ) * S0 / (( 1 - w(ii)^2 / wg^2) ^2 + 4 * ksg^2 * w(ii)^2 / wg^2 );  % 式（5-10）
    S = [ S SS ];
end

%按照目标功率谱生成相应时程
nt = fs * tl + 1;                 % 计算生成时程的采样点数
nfft = 2 ^ nextpow2( nt );        % 计算傅里叶变换点数
dt = 1 / fs ;                     % 采样间隔
df = fs / nfft ;                  % 频率间隔
t = 0 : dt : ( nt - 1 ) * dt;         % 生成时间离散向量
g = randn ( 1 , nfft / 2 + 1 ) * 2 * pi ;        % 生成0~2pi之间按照高斯分布的随机相位
bl = sqrt ( 4 * df * S ) * nfft / 2;              %将目标功率谱转换成对应的傅里叶幅值谱
% bl = sqrt ( 4 * S * df * 2 * pi  ) ;     %公式（5-1）

% 迭代计算
for k = 1 : iter_num
    c = bl .* exp ( i * g );                    % 将幅值谱和相位分别作为一复数序列的实部和虚部
    d = [ c , c( nfft/2 :  -1 : 1 )];                  % 将实部和虚频放在同一个向量中
    e = ifft ( d , nfft );                      % 对向量d进行逆傅里叶变换
    y = real ( e ( 1 : nt ));                   % 取逆傅里叶变换的实部结果作为生成的时程数据
%   [ p_sd,ff ] = pwelch ( y , [] , [] , 256 , 10.24);     % 计算生成时程数据的功率谱
    p_sd = pwelch ( y , [] , [] , nfft , 10.24);     % 计算生成时程数据的功率谱
    dd = bl;
    bl = dd .* sqrt ( S ./ p_sd' );                  % 计算目标傅里叶幅值谱与新生成时程数据的傅里叶幅值谱的比值，并用该比值来修正傅里叶幅值谱
%  按照设定的迭代次数进行迭代计算
end

% 后处理绘图部分

figure(1)                   % 绘制生成的时程曲线图
plot( t,y );
xlabel ('time ( s ) ');
ylabel ('acceleration ( g ) ');
[ p_sd_t , ff ] = pwelch ( y , [] , [] , 256 , 10.24 );


figure(2)                 %绘制目标功率谱和生成的功率谱进行比较
plot ( ff , p_sd_t , 'b' , w , S , 'r' );
%  plot ( w , p_sd , 'b' , w , S , 'r' );
xlabel ( 'Frequency ( Hz ) ');
ylabel ( 'PSD ( m^2 / s^3 )');
axis( [ 0 5.5 0.0005 0.004 ] )
set( gca,'xtick',0:0.5:5.5 )
legend ( '生成时程数据的功率谱','目标功率谱')

huakong

同求啊
另程序中迭代计算那的k都没用到

huakong

不知楼主问题解决没有

dw04116

请教个问题，白噪声的频域表达式怎么样的？

wyuxuenvni

回复 4 # dw04116 的帖子

上边程序中运行出来的图figure（2）就是白噪声的频谱图

wyuxuenvni

回复 2 # huakong 的帖子

恩，目前还没解决呢