组合正弦波曲线拟和.m

基于曲线拟合的频率特性测试方法研究

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%   曲线拟和求频率特性    %%%%%%%%%%%%%%%%%
close all; 
clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%系统离散化及初始值设置%%%%%%%%%%%%%%%%%%%%%
sys=tf([3.355e7],[1 1.504e3 5.394e5 3.291e7]);
Fs=600;
h=1/Fs;
N=6000;
F0=0.1;
F1=60;
df=0.1;
t=[0:h:N*h]; 
%sys=tf([2500],[1 20 2500]);
sysd=c2d(sys,1/Fs,'tustin');
[num,den]=tfdata(sysd,'v');
[h0,ff0]=freqz(num,den,N/2,Fs);
mag=20*log10(abs(h0));
ph=angle(h0);
ph=ph*180/pi;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%输入信号%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=0.5;%频率F1信号的幅度
M=round((F1-F0)/df+1);
ff=[1:M];
AA=[1:M];
PP=[1:M];
dA=[1:M];
dP=[1:M];
x=0;
for i=1:M
  f=F0+df*(i-1);
  ee(i)=i*(i-1)*pi/M;
  x=A*sin(2*pi*f*t+ee(i))+x; 
end;
%x=x+randn(1,N+1);
y=filter(num,den,x);
%y=awgn(y,20);
%y=y+normrnd(0,0.06,1,N+1);
U=zeros(N+1,2*M);
for i=1:M
  f=F0+df*(i-1);
  for j=1:N+1
  U(j,2*i-1)=sin(2*pi*f*t(j)+ee(i));
  U(j,2*i)=cos(2*pi*f*t(j)+ee(i));
 end;
end;
V=zeros(N+1,1);
  for j=1:N+1
      V(j,1)=y(j);
  end;
 C=zeros(2*M,1);
 C=inv(U'*U)*U'*V;
 for i=1:M
 f=F0+df*(i-1);
 k=round(f/(Fs/N));
 c1=C(2*i-1,1);
 c2=C(2*i,1);
 AA(i)=20*log10(sqrt(c1*c1+c2*c2)/A);
 PP(i)=atan(c2/c1)*180/pi;
  if c1<0
   if c2>0
         PP(i)=PP(i)+180;
     else
         PP(i)=PP(i)-180;
     end;
 end;
 dA(i)=AA(i)-mag(k+1);
 dP(i)=PP(i)-ph(k+1);
 ff(i)=f;
 end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%绘制频率特性曲线%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
w0=ff;
semilogx(w0,AA,'b.',ff0(1:N/2),mag,'k');%幅频特性 
grid;ylabel('增益（dB）');xlabel('频率（HZ）');
legend('组合正弦波','理论波形');
title('幅频特性');
axis([0.5 60 -20 1]);
figure(2);
semilogx(w0,PP,'b.',ff0(1:N/2),ph,'k');   %相频特性
grid;ylabel('相角（deg)');
title('相频特性');
legend('组合正弦波','理论波形');xlabel('频率（HZ）');
axis([0.5 60 -180 10]);
figure(3);
w0=ff;
plot(w0,dA);
grid;ylabel('增益误差（dB）');xlabel('频率（HZ）');
%axis([0.1 60 -0.02 0.01]);
figure(4);
plot(w0,dP);
grid;ylabel('相位误差（deg）');xlabel('频率（HZ）');
%axis([0.1 60 -0.15 0.1]);