wyuxuenvni 发表于 2011-10-8 20:55

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

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

问题:
(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 发表于 2012-2-28 01:54

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

huakong 发表于 2012-2-28 01:54

不知楼主问题解决没有

dw04116 发表于 2012-3-2 16:33

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

wyuxuenvni 发表于 2012-3-14 17:32

回复 4 # dw04116 的帖子

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

wyuxuenvni 发表于 2012-3-14 17:33

回复 2 # huakong 的帖子

恩,目前还没解决呢
页: [1]
查看完整版本: 请教个白噪声合成的程序问题,望高手们赐教!