现代信号处理大作业

[Y,fs]=audioread('D:\r.mp3'); L=length(Y); y=awgn(Y,20);%加高斯白噪声 t=0:1:L/2-1; %定义时间的范围和步长 f=7000*t/L; figure(2) subplot(211); plot(Y);xlabel('时间（t）');ylabel('幅度（Y）'); title('加噪前音频语音波形的时域图'); Y=Y(:,1); y=y(:,1); Y1=y+Y; subplot(212); plot(Y1); xlabel('时间（t）');ylabel('幅度（Y）'); title('加噪后音频语音波形的时域图');
现代信号处理新方法
时间：2016.5.3（第9周）
实验一
（1）蒲丰投针实验计算PI值
³ ³ ³ ³ ³
在平面上画有一组间距为a的平行线，将一根长 度为l（l<a）的针任意掷在这个平面上，求此 针与平行线中任一条相交的概率 投掷n次 相交次数为m 相交概率p=m/n=2l/(a* π) π=2l/(a*p)
MATLAB程序
运行结果
0.15 0.1 理想冲击响应 1 0.8 哈明窗
hd[n]
0.05 0 -0.05 0 20 40 60 80
w[n]
0.6 0.4 0.2 0 0 20 40 衰减幅度 60 80
0.15 0.1 0.05 0 -0.05 0
实际冲击响应
100 0
H[w]
20 40 60 80
n 100 500
实验结果（a=5,l=a/2）
pi 3.846153846 3.048780488 3.134796238 3.13283208 3.122073057 3.140309006 3.147722623 3.140703518 3.144673866 3.141929107
1000 5000
MATLAB程序
运行结果
0.5 0.4 0.3 hd[n] w[n] 0.2 0.1 0 -0.1 -0.2 0 5 10 15 20 25 30 35 40 45 50 0.2 0 0.6 0.4 理想冲击响应 1 0.8 哈明窗
0
5
10
15
20
25
30
35
40
45
50
0.5 0.4 0.3
实际冲击响应
fp=100;fz=200;fc=150;fs=2000; wc=2*pi*fc/fs;%将模拟滤波器技术指标转换为数字滤波器技术指标 M=ceil(3.44*fs/100);%求窗函数项数M alpha=(M-1)/2; n=[0:1:(M-1)]; m=n-alpha+eps; hd=sin(wc*m)./(pi*m);%理想冲击响应 w_ham=(hamming(M))';%M项哈明窗 h=hd.*w_ham;%加哈明窗，计算脉冲响应 [H,w]=freqz(h,1);%求频谱 mag=abs(H);%求实部 db=20*log(mag); subplot(221);stem(hd);title('理想冲击响应'); ylabel('hd[n]'); subplot(222);stem(w_ham);title('哈明窗'); ylabel('w[n]'); subplot(223);stem(h);title('实际冲击响应'); ylabel('h[n]'); subplot(224); plot(w/(2*pi),db);title('衰减幅度');ylabel('H[w]');
MATLAB程序
运行结果
4 2
幅度
信号滤波前时域图像
0 -2 -4 0 0.05 0.1 0.15 0.2 0.25
300 200
幅度
时 间 /秒 信号滤波前频域图
100 0
0
100
200
300
400 500 600 频 率 /赫 兹
700
800
900
1000
³ ³ ³ ³
³ ³ ³ ³ ³ ³

³
（2）蒙特卡洛方法计算PI值
³
如图正方形的边长为1,圆的半径为1 可知PI=4*圆面积/正方形面积













n=input('请输入n：');%输入实验次数 x=unifrnd(0,1,1,n); %产生n个在（0,1）之间均匀分布的随机数 y=unifrnd(0,1,1,n);%产生n个在（0,1）之间均匀分布的随机数 j=0; %判断（x,y）和圆的关系，若不在圆内，j=j+1 for i=1:n if x(i)^2+y(i)^2>1; x(i)=0;y(i)=0; j=j+1; end end scatter(x,y,'.','r') %描绘（x,y）散点图 axis equal % 使在每个方向的数据单位都相同 axis([0,1,0,1])%axis([xmin xmax ymin ymax])用来标注输出的图线的 最大值最小值 PI=4*(1-j/n)
MATLAB程序

运行结果
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 1
n
1000 5000 10000 50000 100000
PI
3.1560 3.1368 3.1600
3.1392 3.1475
n=10000运行结果图
实验结论
MATLAB程序
运行结果
2 1
幅度
信号滤波后的时域图
0 -1 -2 0 0.05 0.1 0.15 0.2 0.25
300 200
幅度
时 间 /s 信号滤波后的频域图
100 0
0
100
200
300
400 500 600 频 率 /赫 兹
700
800
900
1000
语音加噪滤波
数字低通滤波器 通带边缘频率fp=1500HZ; 阻带边缘频率fz=2000HZ; 采样频率fs=7000HZ 加入噪声为高斯白噪声
10000 50000
100000 500000
1000000 5000000
实验结论
从上述数据可知，随着模拟次数的 越来越多，PI的值逐渐稳定在π值附近， 即越来越趋近于π
实验一
给定一个变量p,用于描述事件在一定条件下发生 的概率，程序重复地生成这些条件。事件发生 的次数除以总产生次数应该近似地等于概率p
MATLAB程序
运行结果
1 0.5 幅 度 （ Y） 0 -0.5 -1 加噪前音频语音波形的时域图
0
0.5
1
1.5
2
时 间 （ t）
2.5
3
3.5
4
4.5 x 10
5
3 2 幅 度 （ Y） 1 0 -1 -2 -3 0 0.5 1 1.5
加噪后音频语音波形的时域图
2
时 间 （ t）
2.5
3
3.5
³
³ ³
³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³
fp=1500;fz=2000;fc=1750;fs=7000; wc=2*pi*fc/fs;%将模拟滤波器技术指标转换为数字滤波器技术指标 M=ceil(3.44*fs/500);%求窗函数项数M alpha=(M-1)/2; n=[0:1:(M-1)]; m=n-alpha+eps; hd=sin(wc*m)./(pi*m);%理想冲击响应 w_ham=(hamming(M))';%M项哈明窗 h=hd.*w_ham;%加哈明窗，计算脉冲响应 [H,w]=freqz(h,1);%求频谱 mag=abs(H);%求实部 db=20*log(mag); figure(1) subplot(221);stem(hd);title('理想冲击响应'); ylabel('hd[n]'); subplot(222);stem(w_ham);title('哈明窗'); ylabel('w[n]'); subplot(223);stem(h);title('实际冲击响应'); ylabel('h[n]'); subplot(224); plot(w/(2*pi),db);title('衰减幅度');ylabel('H[w]');
sf=filter(h,1,s); %使用filter函数对信号进行滤波 subplot(211); plot(t,sf); %滤波后信号的图像 xlabel('时间/s');ylabel('幅度');title('信号滤波后的时域图 '); subplot(212); Fsf=fft(sf,512); %滤波后信号频域图 AFst=abs(Fsf); %信号频域图的幅值 f=(0:255)*fs/512; % 频率采样 plot(f,AFst(1:256)) %滤波后的信号频域图 xlabel('频率/赫兹'); ylabel('幅度'); title('信号滤波后的频 域图');
4
4.5 x 10
5
³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³
Yn=fft(Y,L); FY=abs(Yn(1:1:L/2)); figure(3) subplot(211) plot(f,FY); grid on; xlabel('频率（f）');ylabel('幅度（FY）'); title('加噪前音频语音波形的频域图'); Y1n=fft(Y1,L); FY1=abs(Y1n(1:1:L/2)); subplot(212); plot(f,FY1); grid on; xlabel('频率（f）');ylabel('幅度（FY1）'); title('加噪后音频语音波形的频域图');
h[n]
-100 -200 -300 0 0.2 0.4 0.6 0.8
叠加的正弦信号
频率分别为f1=50HZ;f2=300HZ;f3=700HZ 采样频率fp=2000HZ
³
³
³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³
f1=50;f2=300;f3=700;%待滤波正弦信号频率 fs=2000;%采样频率 subplot(211); t=0:1/fs:0.25;%定义时间范围和步长 s=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t); plot(t,s);%滤波前的信号图像 xlabel('时间/秒'); ylabel('幅度') title('信号滤波前时域图像'); subplot(212); Fs=fft(s,512);%将信号变换到频域 AFs=abs(Fs);%信号频域图的幅值 f=(0:255)*fs/512;%频率采样 plot(f,AFs(1:256));%滤波前的信号频域图 xlabel('频率/赫兹'); ylabel('幅度') title('信号滤波前频域图')
50 0 -50 H[w] -100 -150 -200
衰减幅度
h[n]
0.2 0.1 0 -0.1 -0.2 0 5 10 15 20 25 30 35 40 45 50
-250
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³
随着实验次数n的增加，实验所得的数值 越来越接近PI,即n越大，结果越精确
实验二
数字滤波器设计
Fra Baidu bibliotek³ ³ ³ ³ ³
设计加窗设计低通FIR滤波器 通带边缘频率fp=100HZ 阻带边缘频率fz=200HZ 采样频率fp=2000HZ 选择哈明窗
³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³ ³
MATLAB程序
³ ³ ³ ³
³
³
a=input('请输入a：');%输入两条平行线间距 n=input('请输入n：');%输入投掷次数 x=unifrnd(0,a/2,[n,1]);%产生n个(0,a/2)之间均匀分 布的随机数 f=unifrnd(0,pi,[n,1]);% 产生n个(0,pi)之间均匀分布的 随机数，pi为投针到最近平行线的角度 y=x<0.25*a*sin(f);%只要满足x<l*sin(f)/2则相交， 这里l=a/2m=sum(y);%相交次数 PI=vpa(a*n/(a*m),10)%求PI，保留10位有效数字
MATLAB程序
运行结果
1500 幅 度 （ FY） 1000 500 0 加噪前音频语音波形的频域图
0
500
1000
3000 幅 度 （ FY1） 2000 1000 0
1500 2000 2500 频 率 （ f） 加噪后音频语音波形的频域图