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
希望有朋友可以帮忙解决一下。
怎么没有人回复呢 本帖最后由 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？
回复 板凳 songzy41 的帖子
谢谢你的回复。我已经做修改了。
对于第二个问题，我觉得如果全部都用单边谱去讲算的话，结果应该是一样的。 本帖最后由 VibInfo 于 2016-11-8 16:15 编辑
原帖由 ap0108220 于 2008-10-15 11:44 发表
谢谢你的回复。
我已经做修改了。
对于第二个问题，我觉得如果全部都用单边谱去讲算的话，结果应该是一样的。
对我上边谈到的笫2个问题：“2，在计算H时，是计算0-Wmax之间的频谱，是单边谱（即只取正频率部分），而X是双边谱（即有正频率部分，也有负频率部分），这两者怎么能一起处理（尽管两者长度一样，概念上说不通）“，是可以单边谱处理，即两者都是单边谱，又等长：Y=H.*X
但在做IFFT之前，要把Y变成双边谱，再做IFFT。
回复 5楼 songzy41 的帖子
那应该如何实现呢？可以把程序写给我看看吧 本帖最后由 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(); 学习了。。。。 大神，学习了。。
页:
[1]