matlab编写的时间序列ARMA模型的信号处理程序

clear all clc; load 'outsjjl'; mn=31; ig=1; %load '';%ig=1，时域振动数据 b=y(1,1:6000); clear y; f=512;%ig=1，f为采样频率，ig=2时为频率间隔 nm=2*mn; if ig==1 sf=f; n=fix(length(b)/2);%n=500 h=b(1,1:2*n)'; clear b; dt=1/sf; t=0:dt:(2*n-1)*dt; else df=f; n=length(b(1,:)); f=0:df:(n-1)*df; H=b(1,:)'; H(n+1)=real(H(n)); H(n+2:2*n)=conj(H(n:-1:2)); h=real(ifft(H)); t=linspace(0,1/df,2*n); dt=t(2)-t(1); end h1=zeros(n,1); [A B]=prony(h,mn,nm); V=roots(B); F1=abs(log(V))/(2*pi*dt); D1=sqrt(1./(((imag(log(V))./real(log(V))).^2)+1)); for k=0:(2*n-1) Va(k+1,:)=[conj(V').^k]; end S1=(pinv(conj(Va')*Va)*conj(Va')*h); vms=Va*S1; h1=real(vms); nn=1:length(t); figure(1); plot(t(nn),h(nn),':',t(nn),h1(nn)); xlabel('时间(s)'); ylabel('幅值'); legend('实测','拟合'); grid on; if ig>1 H1=fft(Va*Sl); nn=1:n; figure(2); subplot(2,1,1); plot(f,real(H(nn)),':',f,real(H1(nn))); xlabel('频率(Hz)'); ylabel('实部'); legend('实测','拟合'); grid on； subplot(2,1,2); plot(f(nn),b(2:nn),':',f(nn),imag(H1(nn))); xlabel('频率(Hz)'); ylabel('虚部'); legend('实测','拟合'); grid on； end [F2,I]=sort(F1); m=0; for k=1:nm-1 if F2(k)~=F2(k+1); continue end m=m+1; l=I(k); F(m)=F1(l); D(m)=D1(l); S(m)=S1(l); end