程序1.1.txt

matlab 计算框架结构模态

%function chengguo11
tf=90;     %<-----------长度加大
delt=0.003;
m=10*[4 4 3 3 3 3 3 3];
M=diag(m);
k=[6 -3 0 0 0 0 0 0
    -3 5 -2 0 0 0 0 0
    0 -2 4 -2 0 0 0 0
    0 0 -2 4 -2 0 0 0
    0 0 0 -2 4 -2 0 0
    0 0 0 0 -2 4 -2 0
    0 0 0 0 0 -2 4 -2
    0 0 0 0 0 0 -2 2];
C=2*k;
K=3000*k;
[E,F]=eig(K,M)
diag (sqrt(F))/2/pi              %<-----------------圆频率除以2pi
y0=rand(8,1);
v0=rand(8,1);
bita=1/6;
Md=inv(M+delt/2*C+bita*delt^2*K);
M1=inv(M);
yn=zeros(8,tf/delt);
xy=zeros(4097,8);
i=1;
for t=0:delt:tf;
    f=[zeros(7,1)
        2*randn(1)];
    if t==0;ydd0=M1*(f-K*y0-C*v0);
    else
        ydd=Md*(f-C*(v0+delt/2*ydd0)-K*(y0+delt*v0+(1/2-bita)*delt^2*ydd0));
        yd=v0+delt/2*(ydd0+ydd);
        y=y0+delt*v0+delt^2/2*ydd0+bita*delt^2*(ydd-ydd0);
        ydd0=ydd;
        v0=yd;y0=y;
    end
    yn(:,i)=y0;
    i=i+1;
end
n=tf/delt;
t=1:n;
nFFT=8192;    % <-----------------FFT 谱线增多
for j=1:8;
    yj=yn(j,t);
   [pxy,f]=csd(yj,yn(8,t),nFFT,1/delt,hanning(nFFT)); 
   figure('name',sprintf('自由度%d的互功率谱密度',j)); %<-------------信号处理工具箱，可以直接获取功率谱密度(psd)图形
    plot(f,10*log10(abs(pxy)))
    xlabel('单位: Hz')
    ylabel('单位: dB')
    xy(:,j)=abs(pxy);   
end                                               %自己用 help psd 命令，看看帮助文
[pxx,f]=pwelch(yn(8,t),nFFT,1/delt,hanning(nFFT)); 
figure('name',sprintf('自由度%d的自功率谱密度',j));
plot(f,10*log10(abs(pxx)))
    xlabel('单位: Hz')
    ylabel('单位: dB')