请教Peakpicking 程序问题

huazi071783

这是我写的一个峰值拾取法的模态识别程序，算出来的参数和理论相差很大，而且得到的振型不对，不知道哪个地方编错了，请高手指点。或者哪位已经编好了可否提供参考。谢谢！！

clear;clc
load('strainsig.mat');                                               %信号是15通道的长度为1024的动态应变信号
fs=2000;                                                                  %采样频率
n=length(strainsig);                                                  %信号的长度
r=fft(strainsig);r=fftshift(r);                                       %对strainsig进行傅立叶变换
amp=abs(r(513:1024,:))*2/n;                                    %幅值谱
f=linspace(0,fs/2,n/2);                                               %生成一个n点的向量(从-fs/2到fs)
t=(0:1023)/fs;
figure(1)
plot(f,amp);      
set(gca,'FontSize',13)
xlabel('Frequence(Hz)');ylabel('Amplitude');
title('Spectrogram')
for k=1:15
        P(k,:)=cpsd(strainsig(:,k),strainsig(:,k),1024);               %求自功率谱密度
end
figure(2)
PSD=log(P);
ff=linspace(0,fs/2,n/2+1);
plot(ff,PSD);
title('PSD')
xlabel('frequence(Hz)');ylabel('PSD(dB)')
fig=0;
f=1;
rr=r(513:1024,:);
while f~=0
            f=input('input a frequency for modal shape or 0 to stop : ');
                if f==0,
                    break
                end
                if fig==3,
                   close(3)
                end
            Ph=180*(abs(angle(rr(round(f),:)))-pi/2)/pi                            %计算相位角
            disp(Ph');
            for i=1:15
               H(i)=(P(i,f)/P(5,f))*sign(Ph(i));                                    %选取5节点作为参考点求振型
            end
            figure(3)
            fig=3;
            hold on
            whitebg(3,'w')                                                               %Change axes background color
            set(3,'position',[200,200,600,400])
            plot([0:16],[0,H,0],'k','linewidth',2)
            plot([0,16],[0,0],'k')
            axis([0,16,-1.2,1.2])                                               %轴线显示大小
end

ChaChing

感觉没原始数据(strainsig.mat), 别人不好测试!?