谐波检测 matlab程序

figure（1）
%载波
t=0:0.0025:60;
u1=5*1.414*sin（2*pi*50*t）;
subplot（211）;
plot（tu1）;
title（‘载波‘）;
grid;
%t调幅波
u2=0.176*sin（2*pi*8.8*t）;
subplot（212）;
plot（tu2）;
title（‘调幅波‘）;
grid;
figure（2）
%闪变信号
w1=（1+u2）.*u1;
subplot（211）;
plot（tw1）;
axis（[00.6-1010]）;
title（‘闪变信号‘）;
grid;
%平方后闪变信号
w2=w1.*w1
subplot（212）;
plot（tw2）;
axis（[00.6080]）;
title（‘平方后闪变信号‘）;
grid;
figure（3）
%35Hz低通滤波器
b10=[ 0 0 0 0 0 0 219.91^6]
a10=[1 848.85 360768.64 9.14*219.91^3 7.46*219.91^4 3.86*219.91^5 219.91^6]
[b1a1]=bilinear（b10a10400）;
y10=filter（b1a1w2）;
subplot（211）;
plot（ty10）;
axis（[05040]）;
title（‘35Hz低通滤波器‘）;
grid;
y1=y10-mean（y10）;
%一阶高通滤波器
b20=[1 0];
a20=[1 0.314];
[b2a2]=bilinear（b20a20400）;
y20=filter（b2a2y1）;
subplot（212）;
plot（ty20）;
axis（[05-2020]）;
title（‘一阶高通滤波器‘）;
grid;
y2=y20/50;
figure（4）
%视感度加权滤波器
num1=[7.0195 100.5500 0];
den1=[9.4392*10^（-4） 1.8531*10^（-1） 1.1121*10 5.0483*10^2 3.3088*10^3];
fs=400;
[n1d1]=bilinear（num1den1fs）;
y30=filter（n1d1y2）;
subplot（211）;
plot（ty30）;
axis（[05-0.40.4]）;
title（‘视感度加权滤波器‘）;
grid;
y3=y30.*y30;
%0.05Hz低通滤波器
num2=[63.7864];
den2=[0.3 1];
fs=400;
[n2d2]=bilinear（num2den2fs）;
y4=filter（n2d2y3）;
subplot（212）;
plot（