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

ap0108220

上面的我的系统模型，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 = [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

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

ap0108220

怎么没有人回复呢

songzy41

原帖由 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

谢谢你的回复。

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

songzy41

原帖由 ap0108220 于 2008-10-15 11:44 发表
谢谢你的回复。

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

对我上边谈到的笫2个问题：“2，在计算H时，是计算0-Wmax之间的频谱，是单边谱（即只取正频率部分），而X是双边谱（即有正频率部分，也有负频率部分），这两者怎么能一起处理（尽管两者长度一样，概念上说不通）“，是可以单边谱处理，即两者都是单边谱，又等长：Y=H.*X
但在做IFFT之前，要把Y变成双边谱，再做IFFT。

ap0108220

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

songzy41

原帖由 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']);

kyu16866

学习了。。。。

yghit08

大神，学习了。。