gonglvpu.m

用功率谱的方法求系统频率特性

Fs=400;%采样频率
N=1000;%采样点数
N1=400;

dfs=Fs/N1;%频率分辨率
n=0:1:N1;
t=n/Fs; %采样时刻
F=([1:N1]-1)*Fs/N1; %换算成实际的频率值
F=F(1:N1/2);%取N/2个实际频率点

%sys=tf([-0.5],[2.0e-6 1.0e-10 1]);
sys=tf([3.355e7],[1 1.504e3 5.394e5 3.291e7]);
sysd=c2d(sys,1/Fs,'tustin');
[num,den]=tfdata(sysd,'v');
[h0,ff0]=freqz(num,den,N/2,Fs);
mag=abs(h0);
ph=angle(h0);
ph=ph*180/pi;

A=1;%频率F1信号的幅度
f0=1;%起始频率(Hz)
df=1;%频率间隔
f1=4000;%结束频率

M=4000;
DELf=(f1-f0)/M;
AA=exp(j*2*pi*f0/Fs);
W=exp(-j*2*pi*DELf/Fs);

S=0;
for i=0:1:((f1-f0)/df)
f=f0+df*i;
S=A*cos(2*pi*f*t)+S;  %组合正弦波
end
x=10*chirp(t,1,N/Fs,100);%chirp
for i=N1:N %补零消除栅栏效应
    S(i)=0;
    x(i)=0;
end;
S=S(1:N);
x=x(1:N);
y=filter(num,den,x);
%ham=hamming(N);
%x=x.*ham';
%y=y.*ham';
cx=xcorr(x);
cx=fft(cx,N);
cy=xcorr(x,y);
cy=fft(cy,N);
for i=1:100
    cc(i)=cy(i)/cx(i);
end;
%p=0.95;
%pxx=psd(S,N,Fs,[],0,p);
%cc=1;
ha=abs(cc);
pa=angle(cc);
pa=-pa*180/pi;
w=F(1:100);
%for i=1:700
%if pa(i)<0
    %pa(i)=0;
%end;
%end

figure(1)
subplot(211);
semilogx(w,20*log10(ha(1:100)),'r',w,20*log(mag(1:100)),'g');
grid;
title('chirp信号');
xlabel('频率（HZ）');
axis([0 100 -40 20]);
subplot(212);
semilogx(w,pa(1:100),'r',w,ph(1:100),'g');
grid;
title('chirp');
xlabel('频率（HZ）');
axis([0 100 -150 10]);

y=filter(num,den,S);
%S=S.*ham';
%y=y.*ham';
cx=xcorr(S);
cx=fft(cx,N);
%cx=czt(S,M,W,AA); %做FFT
cy=xcorr(S,y);
cy=fft(cy,N);
%cy=czt(cy,M,W,AA); %做FFT
for i=1:100
cc(i)=cy(i)/cx(i);
end;
%p=0.95;
%pxx=psd(S,N,Fs,[],0,p);
%cc=1;
ha=abs(cc);
pa=angle(cc);
pa=-pa*180/pi;

Ff=f0:DELf:f0+(M-1)*DELf;
w=F(1:100);
figure(2)
subplot(211);
semilogx(w,20*log10(ha(1:100)),'r',w,20*log10(mag(1:100)),'b');
grid;
title('组合正弦波');
axis([0 100 -20 20]);
xlabel('频率（HZ）');
subplot(212);
semilogx(w,pa(1:100),'r',w,ph(1:100),'b');
grid;
title('组合正弦波');
xlabel('频率（HZ）');
axis([0 100 -150 10]);

dA=[1:100];
dP=[1:100];
for i=1:100
    dA(i)=ha(i)-mag(i);
    dP(i)=pa(i)-ph(i);
end;
figure(3)
subplot(211);
plot(w,dA);
grid;
title('组合正弦波');
%axis([0 100 -0.1 0.1]);
xlabel('频率（HZ）');
subplot(212);
semilogx(w,dP);
grid;
title('组合正弦波');
xlabel('频率（HZ）');
%axis([0 100 -2 2]);