ap0108220 发表于 2008-10-9 11:19

Matlab仿真随机激励下系统的响应




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


我的思路是


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



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



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



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


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


出现的问题:


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







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


激励位置1,响应1





激励1,响应2





激励1,响应3







下面是我的程序:




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 = ;                         % mass matrix
damp = ;    % damp matrix
stiff = ;      % 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

希望有朋友可以帮忙解决一下。

ap0108220 发表于 2008-10-10 15:01

怎么没有人回复呢

songzy41 发表于 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?

ap0108220 发表于 2008-10-15 11:44

回复 板凳 songzy41 的帖子

谢谢你的回复。

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

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。

ap0108220 发表于 2008-10-16 16:23

回复 5楼 songzy41 的帖子

那应该如何实现呢?可以把程序写给我看看吧

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 = ;
damp = ;
stiff = ;
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();
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=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();
axis();
subplot(212)
plot(t,real(y))
xlabel('time(sec)'),ylabel('Amp. y'),
title();

kyu16866 发表于 2012-11-5 17:21

学习了。。。。

yghit08 发表于 2012-12-12 16:44

大神,学习了。。
页: [1]
查看完整版本: Matlab仿真随机激励下系统的响应