请教个白噪声合成的程序问题,望高手们赐教!
最近在学习功率谱合成,编程过程中发现了一些问题,诚恳高手们赐教,小女不胜感激!问题:
(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 ( '生成时程数据的功率谱','目标功率谱')
同求啊
另程序中迭代计算那的k都没用到 不知楼主问题解决没有 请教个问题,白噪声的频域表达式怎么样的? 回复 4 # dw04116 的帖子
上边程序中运行出来的图figure(2)就是白噪声的频谱图 回复 2 # huakong 的帖子
恩,目前还没解决呢
页:
[1]