声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3313|回复: 8

[综合] Matlab仿真随机激励下系统的响应

[复制链接]
发表于 2008-10-9 11:19 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
Snap2.jpg


上面的我的系统模型,M,C, K系统。


我的思路是


先假设系统输入白噪声为x(t),对其进行傅立叶变换为X(w),



然后由M, C,K,得到系统的系统的频域响应函数H(w),



则输出为 Y(w)=H(w)X(w),



然后去Y(w)求逆傅立叶变换得到时域响应y(t).


不知道思路对不对?请论坛大侠们帮助。


出现的问题:


第一,白噪FFT变换后,为何第一个值比其他值都大?结果导至输出响应第一个值都比其他的大。


Snap1.jpg




第二,变换响应点的位置,为何响应信号幅值变化那么大?


激励位置1,响应1


11.JPG


激励1,响应2


12.JPG


激励1,响应3


13.JPG




下面是我的程序:




clear all
clc
close all hidden


% Set parameters

% time domain
fs = 512;                 % Sample frequency
N_Timesample =5120;    % #sample in time domain
T=N_Timesample/fs;


% frequency domain
Wmax = 256;              % max. frequency (rad/sec)
N_Sample = 2560;        % number of frequency samples
dw = Wmax/N_Sample;

Noise_Factor =0;         % noise factor

inp = 1;            % input point
out = 3;            % response point (Change 'out' vaule to get different point's response)


F=5;               % Magnitude of Excitation

% degrees of freedom lineal forced system with viscous damping
m1=0.3;m2=0.3;m3=0.3;
c1=0.01;c2=0.01;c3=0.01;c4=0.01;
k1=1000;k2=1000;k3=1000;k4=1000;


mass = [m1,0,0;0,m2,0;0,0,m3];                         % mass matrix
damp = [c1+c2,-c2,0;-c2,c2+c3,-c3;0,-c3,c3+c4];      % damp matrix
stiff = [k1+k2,-k2,0;-k2,k2+k3,-k3;0,-k3,k3+k4];      % stiff matrix

W=linspace(0,Wmax,N_Sample);
for ii=1:N_Sample
  w=W(ii);
  Kd=-w*w*mass+sqrt(-1)*w*damp+stiff;
  Kinv=inv(Kd);
  H(ii)=Kinv(out,inp)*F;
end


figure(1)
plot(W,real(H)),
xlabel('Frequency (rad/sec)'),ylabel('Real(H)')
% Generate random excitation
t=linspace(0,T,N_Timesample);
rand('state',0)
x=rand(1,N_Timesample);

X=fft(x,N_Timesample/2);

figure(2)
subplot(211)
plot(t,x),xlabel('time(sec)'),ylabel('Amp. x'),
subplot(212),plot(W,real(X))
xlabel('Frequency (rad/sec)'),ylabel('Real(X)'),

% Get the response
Y=X.^H;
y=ifft(Y,N_Timesample);

figure(3)
subplot(211)
plot(W,real(Y))
xlabel('Frequency (rad/sec)'),ylabel('Real(Y)')
subplot(212)
plot(t,real(y))
xlabel('time(sec)'),ylabel('Amp. y'),

%end of program

希望有朋友可以帮忙解决一下。
回复
分享到:

使用道具 举报

 楼主| 发表于 2008-10-10 15:01 | 显示全部楼层
怎么没有人回复呢
发表于 2008-10-10 17:00 | 显示全部楼层
本帖最后由 VibInfo 于 2016-11-8 16:15 编辑
原帖由 ap0108220 于 2008-10-9 11:19 发表
3734


上面的我的系统模型,M,C, K系统。
我的思路是
先假设系统输入白噪声为x(t),对其进行傅立叶变换为X(w)
然后由M, C,K,得到系统的系统的频域响应函数H(w),
则输出为 Y(w)= ...

我对LZ在数字处理的部分感觉有些问题:
1,LZ设置fs=512,后设置Wmax=256,是不是应认为是fs/2,而不是圆频率?单位是Hz而不是弧度/秒?
可能把Wmax作为圆频率,那最大圆频率应为256*2*pi,因为信号最大的频率可以到fs/2。
2,在计算H时,是计算0-Wmax之间的频谱,是单边谱(即只取正频率部分),而X是双边谱(即有正频率部分,也有负频率部分),这两者怎么能一起处理(尽管两者长度一样,概念上说不通)?
3,LZ在前面说明部分写Y(w)=H(w)X(w),后面程序中写为Y=X.^H,是否错了,应改为Y=X.*H?
 楼主| 发表于 2008-10-15 11:44 | 显示全部楼层

回复 板凳 songzy41 的帖子

谢谢你的回复。

我已经做修改了。
对于第二个问题,我觉得如果全部都用单边谱去讲算的话,结果应该是一样的。
发表于 2008-10-15 14:48 | 显示全部楼层
本帖最后由 VibInfo 于 2016-11-8 16:15 编辑
原帖由 ap0108220 于 2008-10-15 11:44 发表
谢谢你的回复。

我已经做修改了。
对于第二个问题,我觉得如果全部都用单边谱去讲算的话,结果应该是一样的。

对我上边谈到的笫2个问题:“2,在计算H时,是计算0-Wmax之间的频谱,是单边谱(即只取正频率部分),而X是双边谱(即有正频率部分,也有负频率部分),这两者怎么能一起处理(尽管两者长度一样,概念上说不通)“,是可以单边谱处理,即两者都是单边谱,又等长:Y=H.*X
但在做IFFT之前,要把Y变成双边谱,再做IFFT。
 楼主| 发表于 2008-10-16 16:23 | 显示全部楼层

回复 5楼 songzy41 的帖子

那应该如何实现呢?可以把程序写给我看看吧
发表于 2008-10-17 14:37 | 显示全部楼层
本帖最后由 VibInfo 于 2016-11-8 16:15 编辑
原帖由 ap0108220 于 2008-10-16 16:23 发表
那应该如何实现呢?可以把程序写给我看看吧

我编的程序如下:
fs=512;
N_Timesample=5120;
T=N_Timesample/fs;
dt=1/fs;
Wmax=256;
N_Sample=2560;
dw=Wmax/N_Sample;
Noise_Factor=0;

inp=1;
out=3;
F=5;
m1=0.3;m2=0.3;m3=0.3;
c1=0.01;c2=0.01;c3=0.01;c4=0.01;
k1=1000;k2=1000;k3=1000;k4=1000;
mass = [m1,0,0;0,m2,0;0,0,m3];
damp = [c1+c2,-c2,0;-c2,c2+c3,-c3;0,-c3,c3+c4];
stiff = [k1+k2,-k2,0;-k2,k2+k3,-k3;0,-k3,k3+k4];
W=(0:N_Sample)*dw;
for ii=1:N_Sample+1
  w=W(ii);
  Kd=-w*w*mass+sqrt(-1)*w*damp+stiff;
  Kinv=inv(Kd);
  H(ii)=Kinv(out,inp)*F;
end
figure(1)
subplot 211; plot(W,real(H));
xlabel('Frequency '),ylabel('Real(H)')
title([num2str(inp) num2str(out) ' System Response']);
subplot 212; plot(W,imag(H));
xlabel('Frequency '),ylabel('Imag(H)')
% Generate random excitation
t=(0:N_Timesample-1)*dt;
randn('state',0)
x=randn(1,N_Timesample);
X=fft(x);
figure(2)
subplot(211)
plot(t,x),xlabel('time(sec)'),ylabel('Amp. x'),
subplot(212),plot(W,abs(X(1:N_Timesample/2+1)));
xlabel('Frequency (Hz)'),ylabel('abs(X)'),
% Get the response
Y=X(1:N_Timesample/2+1).*H;
Y=[Y conj(Y(end-1:-1:2))];
y=ifft(Y);
figure(3)
subplot(211)
plot(W,20*log10(abs(Y(1:N_Timesample/2+1)))); grid;
xlabel('Frequency (Hz)'),ylabel('Log(abs(Y))')
title([num2str(inp) num2str(out) ' Output Log-Response']);
axis([0 120 -80 60]);
subplot(212)
plot(t,real(y))
xlabel('time(sec)'),ylabel('Amp. y'),
title([num2str(inp) num2str(out) ' Time Response']);
发表于 2012-11-5 17:21 | 显示全部楼层
学习了。。。。
发表于 2012-12-12 16:44 | 显示全部楼层
大神,学习了。。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 12:57 , Processed in 0.064152 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表