【公开课】数字信号分析理论与实践（基于matlab） - 华中科技大学（数字信
课程中的matlib版本较低，很多函数发生了变化。在这里补充一下课程中的实例程序。
使用到的控件：按钮、编辑字段（数值）、滑块、坐标区。
需要注意的是：最后老师演示的实例，滑动滑块后自动显示了波形，这是因为老师将程序添加在了滑块的回调函数中。所使用的程序不需要额外改变。
%在标尺中可以设置坐标轴的范围来实现对整体信号的截取，程序中整个声音时长为1s，通过坐标轴范围来显示0.01s中的波形
% 生成音频文件需要为全局变量，通过使用global关键词来设置
global player;
F=app.Slider.Value;
Fs=44100;
dt=1.0/Fs;
T=1;
N=T/dt;
t=linspace(0,T,N);
x=0.3*sin(2*pi*F*t);
plot(app.UIAxes,t,x);
app.EditField.Value = F;
player = audioplayer(x,Fs);%生成音频文件
play(player);%播放音频文件
P16 作业-----------------------------------------------
%用11025Hz的采样频率，用Matlab生成时间长度为1s的500Hz的正弦波
%方波等标准信号，并用plot函数绘制信号波形，用sound函数播放信号声音
Fs=11025;
dt=1.0/Fs;
T=1;
N=T/dt;
t=(0:N-1)/N;
%正弦波
%x=0.3*sin(2*pi*500*t);
%方波
%x=0.3*square(2*pi*500*t);
%锯齿波
%x=0.3*sawtooth(2*pi*500*t,1/2);
%白噪声
%x=0.3*randn(1,N);
plot(t,x);
axis([0,0.01,-0.5,0.5]);
player = audioplayer(x,Fs);%生成音频文件
play(player);%播放音频文件
P19 电子琴加泛音-----------------------------------
只以1号键为例
global player
A=linspace(0,0.9,8800);
D=linspace(0.9,0.8,2200);
S=linspace(0.8,0.8,18000);
R=linspace(0.8,0,15100);
global adsr
adsr=[A,D,S,R];
Fs= 44100;
dt=1.0/Fs;
T=1;
N=T/dt;
t=linspace(0,T,N);
x=0.3*sin(2*pi*247*t);
y=x.*adsr;
plot(app.UIAxes,t,y);
player = audioplayer(y,Fs);
play(player);
项目1：数字信号发生器
效果图：
回调函数列表：
控件列表：
代码如下：
properties (Access = private)%私有属性
Timer_id;
type;
end
methods (Access = private)%私有函数
function my_callback_fcn(app)
global Fs
global x
global N
global F
global A
global z
Fs = 44100;
dt=1.0/Fs;
T=1;
N=T/dt;
t=linspace(0,T,N);
F=app.FreSlider.Value;
A=app.AmpSlider.Value;
if app.type==1
x=0.5*A*randn(1,N);
end
if app.type==2
x=A*sin(2*pi*F*t+z);
end
if app.type==3
x=A*square(2*pi*F*t+z);
end
if app.type==4
x=A*sawtooth(2*pi*F*t+z);
end
z=z+0.25;
plot(app.UIAxes,t,x);
grid(app.UIAxes);
end
%定时器初始化，重复模式
function timer_init(app)
app.Timer_id=timer;
app.Timer_id.Period=1.0; %周期
app.Timer_id.ExecutionMode='fixedRate';
app.Timer_id.TasksToExecute=Inf;
app.Timer_id.TimerFcn=@(~,~) my_callback_fcn(app);
end
%定时器启动
function timer_start(app)
start(app.Timer_id);
end
%定时停止
function timer_stop(app)
stop(app.Timer_id);
end
%删除定时器
function timer_delete(app)
delete(app.Timer_id);
end
end
% Callbacks that handle component events
methods (Access = private)
% Value changed function: FreSlider
function FreSliderValueChanged(app, event)
app.EditFreEditField.Value=event.Value;
end
% Value changed function: AmpSlider
function AmpSliderValueChanged(app, event)
app.EditAmpEditField.Value=event.Value;
end
% Value changed function: EditFreEditField
function EditFreEditFieldValueChanged(app, event)
app.FreSlider.Value=event.Value;
end
% Value changed function: EditAmpEditField
function EditAmpEditFieldValueChanged(app, event)
app.AmpSlider.Value=event.Value;
end
% Button pushed function: PlayButton
function PlayButtonPushed(app, event)
global player
global Fs
global x
y=x*0.01;
player=audioplayer(y,Fs);
play(player);
end
% Button pushed function: NoiseButton
function NoiseButtonPushed(app, event)
app.type=1;
my_callback_fcn(app);
end
% Button pushed function: SineButton
function SineButtonPushed(app, event)
app.type=2;
my_callback_fcn(app);
end
% Button pushed function: SquareButton
function SquareButtonPushed(app, event)
app.type=3;
my_callback_fcn(app);
end
% Button pushed function: TrangleButton
function TrangleButtonPushed(app, event)
app.type=4;
my_callback_fcn(app);
end
% Button pushed function: RunButton
function RunButtonPushed(app, event)
timer_init(app);
timer_start(app);
end
% Button pushed function: StopButton
function StopButtonPushed(app, event)
ts=timerfind;
if ~isempty(ts)
timer_stop(app);
timer_delete(app);
end
end
end
P24- 作业----------------------------------------------
本人将两个作业合并为一个。
效果图：
所使用到的控件：
程序代码如下：
properties (Access = private)
x % Description
dt
N
t
Fs=5120;
type_num=1;
end
methods (Access = private)
function wave_sin(app)
pinlv=app.Slider.Value;
fuzhi=app.Slider_2.Value;
app.dt=1.0/app.Fs;
T=1;
app.N=T/app.dt;
app.t=(0:app.N-1)/app.N*T;
app.x=fuzhi*sin(2*pi*pinlv*app.t);
jisuan(app);
end
function daoshu(app)%对信号先微分再积分
y1=app.x;
y1(1)=app.x(1);
for k=2:1:app.N-1
y1(k)=(app.x(k+1) - app.x(k-1))/(2*app.dt);
end
y1(app.N)=app.x(app.N);
plot(app.UIAxes2,app.t,y1);
y2=y1;
y2(1)=0;
for k=2:1:app.N-1
y2(k)=y2(k-1) + (app.dt*(y1(k)+y1(k-1)))/2;
end
plot(app.UIAxes3,app.t,y2);
end
function jisuan(app)
%分析幅值
fz=max(app.x);
app.EditField.Value=fz;
%分析峰峰值
ffz=max(app.x)-min(app.x);
app.EditField_2.Value=ffz;
%分析平均值
pjz=mean(app.x);
app.EditField_3.Value=pjz;
%分析有效值
cz=app.x-pjz;%波形数组每一个元素与平均值的差值
czpf=cz.^2;%差值平方
czpfpj=sum(czpf)/app.N;%差值平方求和取平均
yxz=sqrt(czpfpj);%开方求得有效值
app.EditField_5.Value=yxz;
%标准差
app.EditField_6.Value=std(app.x);
%频率和初始相位
p=max(app.x);q=min(app.x);n=1;ti=(0);
at=0.8*(p-q)+q;
for k=2:1:app.N-1
if(app.x(k-1)<at && app.x(k)<=at && app.x(k+1)>at && app.x(k+2)>at)
ti(n) = k;
n = n+1;
end
end
T=((ti(2) - ti(1)))*app.dt;
F=1.0/T;
Q=(T-ti(1)*app.dt);%计算出初始相位
app.EditField_9.Value=Q;
app.EditField_4.Value=F;
plot(app.UIAxes,app.t,app.x);
end
function wave_square(app)
pinlv=app.Slider.Value;
fuzhi=app.Slider_2.Value;
app.dt=1.0/app.Fs;
T=1;
app.N=T/app.dt;
app.t=(0:app.N-1)/app.N*T;
app.x=fuzhi*square(2*pi*pinlv*app.t);
jisuan(app);
end
function wave_sawtooth(app)
pinlv=app.Slider.Value;
fuzhi=app.Slider_2.Value;
app.dt=1.0/app.Fs;
T=1;
app.N=T/app.dt;
app.t=(0:app.N-1)/app.N*T;
app.x=fuzhi*sawtooth(2*pi*pinlv*app.t);
jisuan(app);
end
function select(app)
switch (app.type_num)
case 1
try
wave_sin(app);
daoshu(app);
catch ErrorInfo
warndlg('数据输入有误！！');
end
case 2
try
wave_square(app);
daoshu(app);
catch ErrorInfo
warndlg('数据输入有误！！');
end
case 3
try
wave_sawtooth(app);
daoshu(app);
catch ErrorInfo
warndlg('数据输入有误！！');
end
end
end
end
% Callbacks that handle component events
methods (Access = private)
% Value changed function: Slider_2
function Slider_2ValueChanged(app, event)
app.EditField_7.Value = app.Slider_2.Value;
select(app);
end
% Value changed function: EditField_7
function EditField_7ValueChanged(app, event)
app.Slider_2.Value = app.EditField_7.Value;
select(app);
end
% Value changed function: Slider
function SliderValueChanged(app, event)
app.EditField_8.Value = app.Slider.Value;
select(app);
end
% Value changed function: EditField_8
function EditField_8ValueChanged(app, event)
app.Slider.Value = app.EditField_8.Value;
select(app);
end
% Button pushed function: Button
function ButtonPushed(app, event)
try
app.type_num=1;
wave_sin(app);
daoshu(app);
catch ErrorInfo
warndlg('数据输入有误！！');
end
end
% Button pushed function: Button_2
function Button_2Pushed(app, event)
try
app.type_num=2;
wave_square(app);
daoshu(app);
catch ErrorInfo
warndlg('数据输入有误！！');
end
end
% Button pushed function: Button_3
function Button_3Pushed(app, event)
try
app.type_num=3;
wave_sawtooth(app);
daoshu(app);
catch ErrorInfo
warndlg('数据输入有误！！');
end
end
end
P27---------------------------------------------------
Fs=5120;N=1024;
dt=1.0/Fs;T=dt*N;
t=linspace(0,T,N);
x=10*sin(2*pi*100*t)+10/3*sin(2*pi*3*100*t);
plot(t,x);
title('原始信号');
y=fft(x,N);
a=real(y);b=imag(y);
figure;%创建图窗窗口
subplot(2,1,1);
plot(a);
title('实频谱');
subplot(2,1,2);
plot(b);
title('虚频谱');
f=linspace(0,Fs/2,N/2);%x坐标换为频率
A1=abs(y)/(N/2);%幅值量纲还原
Q1=angle(y)*180/pi;
figure;%创建图窗窗口
subplot(2,1,1);
plot(f,A1(1:N/2));%不显示负频率部分
title('幅频谱');
subplot(2,1,2);
plot(f,Q1(1:N/2));%不显示负频率部分
title('相频谱');
A2=A1.^2;%功率谱
P2=20*log10(A2);%对数功率谱
figure;
subplot(2,1,1);
plot(f,A2(1:N/2));
title('功率谱');
subplot(2,1,2);
plot(f,P2(1:N/2));
title('对数功率谱');
P27 作业---------------------------------------------
Fs=11025;N=1024;
dt=1.0/Fs;T=dt*N;
t=linspace(0,T,N);
%正弦波
x=10*sin(2*pi*200*t);
y=fft(x,N);
f=linspace(0,Fs/2,N/2);%x坐标换为频率
A1=abs(y)/(N/2);%幅值量纲还原
subplot(2,1,1);
plot(t,x);
title('正弦波');
subplot(2,1,2);
plot(f,A1(1:N/2));%不显示负频率部分
title('幅频谱');
%方波
x1=10*square(2*pi*200*t);
y1=fft(x1,N);
A2=abs(y1)/(N/2);%幅值量纲还原
figure;
subplot(2,1,1);
plot(t,x1);
title('方波');
subplot(2,1,2);
plot(f,A2(1:N/2));%不显示负频率部分
title('幅频谱');
%三角波
x2=10*sawtooth(2*pi*200*t);
y2=fft(x2,N);
A3=abs(y2)/(N/2);%幅值量纲还原
figure;
subplot(2,1,1);
plot(t,x2);
title('三角波');
subplot(2,1,2);
plot(f,A3(1:N/2));%不显示负频率部分
title('幅频谱');
P34------------------------------------------------------
【备注】：只是完成了课程中老师所提及到的功能，更多细节可以再优化优化，就让各位读者自行探索。
效果图：
程序：
properties (Access = private)
Timer_id
wave
num=0;
N=24810; % 24810 = length(doubleArray)/2
Z1=0.1;
Z2=0.1;
end
methods (Access = private)
%音频处理
function Recording(app)
recObj = audiorecorder(11025,16,1);
recordblocking(recObj,5);
doubleArray = getaudiodata(recObj);
doubleArray = doubleArray(5506:end);%去除录音数据头
ylim(app.UIAxes,[-app.Z1 app.Z1]);
plot(app.UIAxes,doubleArray);
switch app.num
case 0
pinlv(app,doubleArray);
case 1
%注意两个数组的大小
H=2*hamming(app.N*2).*doubleArray;%2为修正系数
pinlv(app,H);
case 2
B=2.381*blackman(app.N*2).*doubleArray;
pinlv(app,B);
case 3
P=4.545*flattopwin(app.N*2).*doubleArray;
pinlv(app,P);
end
end
%fft
function pinlv(app,wave)
y_fft=fft(wave);
f=linspace(0,5505,app.N);
A=abs(y_fft)/app.N;
ylim(app.UIAxes2,[0 0.01*app.Z2]);
plot(app.UIAxes2,f,A(1:app.N));
end
%定时器初始化，重复模式
function timer_init(app)
app.Timer_id=timer;
app.Timer_id.Period=5.2; %周期
app.Timer_id.ExecutionMode='fixedRate';
app.Timer_id.TasksToExecute=Inf;
app.Timer_id.TimerFcn=@(~,~) Recording(app);
end
%定时器启动
function timer_start(app)
start(app.Timer_id);
end
%定时停止
function timer_stop(app)
stop(app.Timer_id);
end
%删除定时器
function timer_delete(app)
delete(app.Timer_id);
end
end
% Callbacks that handle component events
methods (Access = private)
% Button pushed function: Button
function ButtonPushed(app, event)
timer_init(app);
timer_start(app);
end
% Button pushed function: Button_2
function Button_2Pushed(app, event)
timer_delete(app);
end
% Button pushed function: Button_3
function Button_3Pushed(app, event)
app.num = 1;
end
% Button pushed function: Button_4
function Button_4Pushed(app, event)
app.num = 2;
end
% Button pushed function: Button_5
function Button_5Pushed(app, event)
app.num = 3;
end
% Value changed function: Slider
function SliderValueChanged(app, event)
app.Z1 = 0.1*2*app.Slider.Value;
end
% Value changed function: Slider_2
function Slider_2ValueChanged(app, event)
app.Z2 = 0.1*2*app.Slider_2.Value;
end
end
P36 第五章-----------------------------------------
信号的时差域相关分析
信号的相关分析是一种分析两个信号之间或一个信号自身的时间依存关系和相似程度的方法。
工程上，人们关心的是信号不同时刻的相似程度，但不太关心其具体值，这时相关函数可简化为：
x=y 自相关；x!=y 互相关
数字信号离散计算公式：
声音素材网址：https://sc.chinaz.com/yinxiao/220314192200.htm
程序：
[y,Fs] = audioread('beat.wav');%读取声音文件，需要在同一文件夹下
y1=y(:,1);%取单声道
subplot(2,1,1);
plot(y1);
s=xcorr(y1,'unbiased');%自相关函数
% ‘unbiased' 避免重叠失真
subplot(2,1,2);
plot(s);
结果：下图中相关函数图形，横坐标中点为零点。
相关函数的性质：
（1）自相关函数是偶函数
（2）当t=0时，自相关函数具有最大值
（3）周期信号的自相关函数仍然是同频率的周期信号，但不保留原信号的相位信息
（4）随机噪声信号的自相关函数将随t的增大快速衰减
（5）两周期信号的互相关函数仍然是同频率的周期信号，且保留原有信号的相位信息
（6）两个非同频率的周期信号互不相关
P40----------------------------------------------------
在数学中，概率密度函数(Probability Density Function)是一个描述信号的取值在某个确定的取值点附近的可能性的函数，定义为概率分布函数的倒数。
以幅值大小为横坐标，以每个幅值间隔内出现的概率为纵坐标进行统计分析的方法。它反映了信号落在不同幅值强度区域内的概率情况。
概率分布函数
概率分布函数(Cumulative distribution function)是信号幅值小于或等于某值R的概率，其定义为：
概率分布函数又称之为累积概率，表示了落在某一区间的概率。
直方图
以幅值为横坐标，以每个幅值间隔内出现的频次为纵坐标进行统计分析，称为直方图。
直方图归一化后可转换为概率密度函数。
直方图均衡，主要作用是提高图像的清晰度和纹理的识别度，可能会使得色彩失真，更多的是用在灰度图像的处理。
P42------------------------------------------------------
效果图：
代码：
function ButtonValueChanged(app, event)
[filename, filepath]=uigetfile({'*.bmp;*.jpg;*.png'},'选择图片文件'); %选择文件
if isequal(filename,0) || isequal(filepath,0)%判断是否选中图片
errordlg('没有选中文件，请重新选择','出错');
return;
else
set(app.EditField,'Value',[filepath, filename]);%获取文件的路径和名字
end
road=get(app.EditField,'Value');
set(gcbf,'UserData',road);% GUI中参数的传递，详情https://www.cnblogs.com/jmliao/p/5628521.html
end
% Value changed function: Button_2
function Button_2ValueChanged(app, event)
road=get(gcbf,'UserData');
app.Image.ImageSource=road;%显示原图像
O_image=imread(road);%读入原始图像
Image_gray=im2gray(O_image);%转换为灰度图像
imwrite(Image_gray,'myGray.png');%将灰度图像保存，图像控件只能显示真色彩图像，暂未找到直接显示灰度图像的方法
app.Image_3.ImageSource='myGray.png';
Image_A_gray=histeq(Image_gray);%直方图均衡强化对比度
imwrite(Image_A_gray,'myHisteq.png');
app.Image_4.ImageSource='myHisteq.png';
end
P43------------------------------------------------
白噪声信号的幅值域分析
效果图：
程序：
x=-4:0.02:4;
y=randn(10000,1);
subplot(4,1,1);
plot(y);
[N, edges] = histcounts(y,x);%N-bin中的计数，edges-bin的边界
centers = (edges(1:end-1) + edges(2:end))./2;%取边界中心值
subplot(4,1,2);stem(centers,N); %绘制离散序列数据
pdf=N/length(y);%归一化,也可以直接使用示例程序中的方式进行归一化
subplot(4,1,3);plot(centers,pdf);
subplot(4,1,4);cdfplot(y);%累计分布函数图像
补充：(可以参考一下，便于理解直方图的参数和属性）
示例程序（来源网络）：
x = randn(1,1000); %Generate some data to plot
figure(1); clf;
subplot(2,1,1); hold on;
h = histogram(x,'normalization','probability');
title('Plotted as a histogram');
ylabel('Frequency');
subplot(2,1,2); hold on;
[N, edges] = histcounts(x,'normalization','probability');
centers = (edges(1:end-1) + edges(2:end))./2; %histcounts gives the bin edges, but we want to plot the bin centers
plot(centers,N,'ko');
title('Plotted as points');
ylabel('Frequency');
自己画的，便于理解：
直方图，概率质量函数和概率密度函数
https://blog.csdn.net/lovehua365/article/details/79058815
频谱、能谱、功率谱、倍频程谱、1/3 倍频程谱
https://blog.csdn.net/liyuanbhu/article/details/42675765
2倍频程、1倍频程、1/3倍频程介绍
https://blog.csdn.net/chenhuanqiangnihao/article/details/122079324
P46 信号的数字滤波-------------------------------
测量中除传感器信号外，干扰也会出现，他们与信号叠加在一起，扭曲测量结果。滤波器是一种选频装置，可以使信号中特定频率成分通过，而衰减其他频率成分。
滤波有一个前提条件，信号的频率是不重叠的。
P47---------------------------------------------
数字滤波和模拟滤波的区别
滤波位置的不同；采样频率不同；有效A/D位数；
当干扰信号频率较高时，应优先模拟滤波将高频的干扰信号滤除，这样可以降低采样频率。
当干扰信号的量程较大时，直接进行数字滤波会使得A/D满量程精度较小。
滤波器评价指标：过渡区陡峭；脉冲响应紧支集。
P48---------------------------------------------
Fs=2048;dt=1.0/Fs;T=1;
N=T/dt;t=[0:N-1]/N;
x1=sin(2*pi*50*t)+sin(2*pi*300*t)+sin(2*pi*500*t);
subplot(4,1,1);
plot(t,x1);axis([0,0.1,-2,2]);
P=fft(x1,N);
Pyy=2*sqrt(P.*conj(P))/N;%幅值量纲还原
f=linspace(0,Fs/2,N/2);
subplot(4,1,2);plot(f,Pyy(1:N/2));
for k=1:N
P1(k)=real(P(k))+i*(imag(P(k)));
end
for k=400:N-400
P1(k)=0;
end
for k=1:200
P1(k)=0;
end
for k=N-200:N
P1(k)=0;
end
Pyy=2*sqrt(P1.*conj(P1))/N;
f=linspace(0,Fs/2,N/2);
subplot(4,1,3);plot(f,Pyy(1:N/2));
x2=ifft(P1);
subplot(4,1,4);plot(t,x2);
axis([0,0.1,-2,2]);
P49---------------------------------------------
%数字差分--简单的高通滤波器
Fs=5000;a=5;
f=2;T=1;
dt=1.0/Fs;
N=T/dt;
t=linspace(0,T,N);
y=a*sin(2*3.14*f*t)+0.3*sin(100*3.14*f*t);
plot(t,y);
x(1)=0;x(N)=0;
for i=2:N-1
x(i)=(y(i+1)-y(i-1))/2;
end
figure;
plot(t,x);
%数据平滑--简单的低通滤波器
Fs=5000;a=5;
f=2;T=1;
dt=1.0/Fs;
N=T/dt;
t=linspace(0,T,N);
y=a*sin(2*pi*f*t)+0.8*sin(1000*pi*f*t);
subplot(2,1,1);
plot(t,y);
x=y;
for i=1:N-5
x(i)=[y(i)+y(i+1)+y(i+2)+y(i+3)+y(i+4)]/5;
end
subplot(2,1,2);
plot(t,x);
作业
Fs=5120;
T=1;
dt=1.0/Fs;
N=T/dt;
t=linspace(0,T,N);
y=sin(2*pi*10*t)+sin(2*pi*500*t);
subplot(5,1,1);
plot(t,y);
%5点平滑滤波
x=y;
for a=1:N-5
x(a)=[y(a)+y(a+1)+y(a+2)+y(a+3)+y(a+4)]/5;
end
subplot(5,1,2);
plot(t,x);
%频域滤波
P=fft(y,N);
Pyy=sqrt(P.*conj(P))/N;%幅值量纲还原
d=linspace(0,Fs,N);
subplot(5,1,3);
plot(d,Pyy);
for k=1:N
P1(k)=real(P(k))+i*(imag(P(k)));
end
for k=200:N-200
P1(k)=0;
end
Pyy1=sqrt(P1.*conj(P1))/N;
d=linspace(0,Fs,N);
subplot(5,1,4);
plot(d,Pyy1);
x2=ifft(P1);
subplot(5,1,5);plot(t,x2);
P54--------------------------------------------------------------------------
Fs=2048;dt=1/Fs;
T=1;N=T/dt;t=[0:N-1]/N;
x1=sin(2*pi*50*t)+sin(2*pi*300*t)+sin(2*pi*500*t);
subplot(3,1,1);plot(t,x1);
axis([0,0.1,-2,2]);P=fft(x1,N);
Pyy = 2*sqrt(P.*conj(P))/N;
f=linspace(0,Fs/2,N/2);
subplot(3,1,2);plot(f,Pyy(1:N/2));
b=fir1(48,0.1);
x2=filter(b,1,x1);
subplot(3,1,3);plot(t,x2);axis([0,0.1,-2,2])
figure;
c=fir1(48,[0.2,0.4]);
x3=filter(c,1,x1);
subplot(2,1,1);plot(t,x3);axis([0,0.1,-2,2])
d=fir1(48,0.4,'high');
x4=filter(d,1,x1);
subplot(2,1,2);plot(t,x4);axis([0,0.1,-2,2])
%% 原始信号
Fs=5120;dt=1/Fs;
T=1;N=T/dt;t=[0:N-1]/N;
%含有低(100Hz)、中(600Hz)、高(1000Hz)三个不同频率成分的测试信号
x1=sin(2*pi*100*t)+sin(2*pi*600*t)+sin(2*pi*1000*t);
subplot(2,1,1);plot(t,x1);
axis([0,0.1,-3,3]);P=fft(x1,N);
Pyy = 2*sqrt(P.*conj(P))/N;
f=linspace(0,Fs/2,N/2);
subplot(2,1,2);plot(f,Pyy(1:N/2));
%% fir1
figure;
b=fir1(48,0.1);
%48-代表滤波器的阶数，指在滤波器的传递函数中有几个极点，同时决定了转折区的下降速度
x2=filter(b,1,x1);
subplot(3,1,1);plot(t,x2);axis([0,0.1,-2,2])
c=fir1(48,[0.1,0.3]);
x3=filter(c,1,x1);
subplot(3,1,2);plot(t,x3);axis([0,0.1,-2,2])
d=fir1(48,0.3,'high');
x4=filter(d,1,x1);
subplot(3,1,3);plot(t,x4);axis([0,0.1,-2,2])
%% fir2
figure
f1=[0 0.1 0.1 1];%低通
mhi1=[1 1 0 0];
b1=fir2(48,f1,mhi1);
x2_fir2=filter(b1,1,x1);
subplot(3,1,1);plot(t,x2_fir2);axis([0,0.1,-2,2])
f2=[0 0.1 0.1 0.3 0.3 1];
mhi2=[0 0 1 1 0 0];
c1=fir2(48,f2,mhi2);
x3_fir2=filter(c1,1,x1);
subplot(3,1,2);plot(t,x3_fir2);axis([0,0.1,-2,2])
f3=[0 0.3 0.3 1];
mhi3=[0 0 1 1];
d1=fir2(48,f3,mhi3);
x4_fir2=filter(d1,1,x1);
subplot(3,1,3);plot(t,x4_fir2);axis([0,0.1,-2,2])
%% firls
figure
% f1=[0 0.1 0.1 1];%低通
% mhi1=[1 1 0 0];
b2=firls(48,f1,mhi1);
x2_firls=filter(b2,1,x1);
subplot(3,1,1);plot(t,x2_firls);axis([0,0.1,-2,2])
% f2=[0 0.1 0.1 0.3 0.3 1];
% mhi2=[0 0 1 1 0 0];
c2=firls(48,f2,mhi2);
x3_firls=filter(c2,1,x1);
subplot(3,1,2);plot(t,x3_firls);axis([0,0.1,-2,2])
% f3=[0 0.3 0.3 1];
% mhi3=[0 0 1 1];
d2=firls(48,f3,mhi3);
x4_firls=filter(d2,1,x1);
subplot(3,1,3);plot(t,x4_firls);axis([0,0.1,-2,2])
