东北大学现代信号处理新方法郎俊大作业

>> l=0:0.001:1; >> for i=1:length(l) f(i)=vpi(l(i),10^6); end >> scatter(l,f,'x');hold on;plot(get(gca, 'XLim'),[pi pi], '-r')
3.3
3.25
3.2
3.15
3.1
3.05
0
0.1
滤波前时域图 0.25 0.2 0.15 0.1 0.05 0.002 0.004 0.006 0.008 0.01 0 0 5000 10000 15000 滤波前频域图
2 0 -2 0
2 1 0 -1 -2
滤波后时域图
0.2 0.15 0.1 0.05 0
滤波后频域图
0
0.002
0.004
0.006
>> fp=2500;fst=3000;fs=24000; >> N=ceil(3.44*fs/(fst-fp)); >> if mod(N,2)==0 N=N+1; end >> n=-(N-1)/2:1:(N-1)/2; >> w=0.54+0.46.*cos(2*pi*n/(N-1)); >> f0=fp+(fst-fp)/2; >> O=2*pi*f0/fs; >> h1=sin(O*n)./(n*pi); >> h1((N+1)/2)=O/pi; >> h=h1.*w; >> figure,subplot(121);plot(h);title('滤波器 h(n)'); >> y=20*log(abs(fft(h,1024))); >> subplot(122);plot(y(1:1024/2));title('频率响应');
2
3.148 3.146 3.144 3.142 3.14 3.138 3.136 3.134
1
2
3
4
5
6
7
8
9
10 x 10
6
据图可得：在一定范围内，随着 n 的增大，估计的 pi 值越接近实际值；当 n 足够大时（大 概大于 8*10^6） ，|vpi-pi|将不随着 n 的增大而明显变化。
l/a 与 pi 估值关系
1
题目一
如何用计算机模拟一个随机事件， 比如蒲丰投针的实验， 并估计 随机事件发生的概率以计算圆周率 pi， 并分析所设计的随机实验中各 因素的影响。
1.模拟蒲丰掷针实验 在平面上画有等距离的一些平行线，平行线间的距离为 a(a>0)，向平面上随机投 一长为 l(l<a)的针。针与平行线相交的概率为 p=2l/(πa) ,则π=2l/(a*p).
function re=vpi(l,n) %投针法估算 pi a=1; x=unifrnd(0,a/2,1,n); %投针的中点到最近的平行线的距离 y=unifrnd(0,pi,1,n); %投针到最近的平行线的角度 m=0; for i=1:n if x(i)<l*sin(y(i))/2 m=m+1; end end p=m/n; re=2*l/(a*p); end
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
据图可得：当 l/a 较小时，估计值与实际值偏差较大；随着 l/a 的增大，估计值与实际值越 来越接近；当 l/a 的值稍大之后（大概大于 0.8） ，估计值不随 l/a 的增大明显变化。
2.蒙特卡洛方法计算π值
3
取一正方形以其一个顶点为圆心，边长为半径画圆，取四分之一圆。1/4 圆 的面积与正方形的面积之比为 pi/4。 用随机试验模拟在正方形内投点， 根据在圆内的点的个数与在正方形内的点 的个数之比得到两者面积比的近似值。再乘以 4 即可得到近似的 pi 值。
3.146 3.145 3.144 3.143 3.142 3.141 3.14 3.139 3.138
1
2
3
4
5
6
7
8
9
10 x 10
6
据图可得：同投针法类似，在一定范围内，随着 n 的增大，估计值越接近实际值；当 n 足够 大时（大概大于 6*10^6） ，随着 n 的增大，估计值将不再随 n 明显变化。
试验次数与 pi 估值关系
>> t=10^6:10000:10^7; >> L=length(t); >> for i=1:L m(i)=vpi(0.8,t(i)); end >> scatter(t,m,'x');hold on;plot(get(gca, 'XLim'),[pi pi], '-r')
东
北
大
学
考 试 试 卷
评 分
考试科目： 课程编号： 阅 卷 人： 考试日期： 姓 学 名： 号：
现代信号处理新方法
郎俊 2016 年 05 月 03 日 哈哈哈 哈哈哈
注 意 事 项
1．考 前 研 究 生 将 上 述 项 目 填 写 清 楚 2．字 迹 要 清 楚，保 持 卷 面 清 洁 3．交 卷 时 请 将 本 试 卷 和 题 签 一 起 上 交
滤 波 器 h(n) 频率响应
0.3 0.25 0.2 0.15
50 0 -50 -100 -150
0.1 -200 0.05 0 -0.05 -250 -300 -350
0
50
100
150
200
0
100
200
300
400
500
600
5
2.叠加正弦波通过滤波器
>> f1=1000;f2=3000;f3=4000; >> subplot(221); >> t=0:1/fs:0.01; >> fx=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t); >> plot(t,fx); >> axis([0 0.01 -3.5 3.5]); >> title('滤波前时域图'); >> subplot(222); >> x=(0:511)*fs/1024;y=abs(fft(fx,1024)); >> plot(x,y(1:512)*2./1024); >> title('滤波前频域图'); >> fh=filter(h,1,fx); >> subplot(223);plot(t,fh); >> title('滤波后时域图'); >> subplot(224); >> y2=abs(fft(fh,1024)); >> plot(x,y2(1:512)*2./1024); >> title('滤波后频域图');
function re=vpi2(n) %蒙特卡洛法估计 pi 值 x=rand(n,1); y=rand(n,1); m=0; for i=1:n if(x(i)^2+y(i)^2<=1) m=m+1; end end re=4*m/n; End
试验次数与 pi 估值关系
>> t=10^6:10000:10^7; >> for i=1:length(t) m(i)=vpi2(t(i)); end >> scatter(t,m,'x');hold on;plot(get(gca, 'XLim'),[pi pi], '-r')
7
0.5 0 -0.5 0
原信号时域图
1000 500
原信号频域图
2
4
6
0
0
0.5
1
1.5
2
2.5 x 10
4wk.baidu.com
0.5 0 -0.5 0
加噪后时域图
1000 500
加噪后频域图
2
4
4
题目二：数字滤波器设计
利用窗函数法、频率取样法以及优化设计方法设计 FIR 滤波器，绘制出滤波 器的特性图。利用所设计的滤波器对多个频带叠加的正弦信号进行处理，对比 滤波前后的信号时域和频域图，验证滤波器的效果。最后录制或下载一段语音 信号，对其进行采样和加噪，然后用设计滤波器对加噪后信号进行滤波，绘制 出信号的时域波形和频谱，分别计算均值、方差、自相关函数、功率谱，并对 滤波前后的信号进行对比。 1.低通滤波器的设计
东北大学研究生院
目录
题目一.................................................................................................................................. 2 1.模拟蒲丰掷针实验...................................................................................................2 2.蒙特卡洛方法计算π值.......................................................................................... 3 题目二：数字滤波器设计.................................................................................................5 1.低通滤波器的设计...................................................................................................5 2.叠加正弦波通过滤波器.......................................................................................... 6 3.音频信号处理........................................................................................................... 7 题目三.................................................................................................................................. 9 题目四................................................................................................................................ 10
0.008
0.01
0
5000
10000
15000
6
3.音频信号处理
>> clear >> [data,fs]=audioread('F:\Media\Music\audiotest.wav'); >> data=data(:,1);%取单声道 >> L=length(data); >> tn=0:1/fs:(L-1)/fs; >> f=0:fs/L:fs*(L-1)/L; >> yn=awgn(data,30);%加高斯白噪声 SNR=30dB >> F1=abs(fft(data,L)); >> F2=abs(fft(yn,L));%频域 %滤波器 >> fp=2500;fst=3000; >> N=ceil(3.44*fs/(fst-fp)); >> if mod(N,2)==0 N=N+1; end >> n=-(N-1)/2:1:(N-1)/2; >> w=0.54+0.46.*cos(2*pi*n/(N-1)); >> f0=fp+(fst-fp)/2; >> O=2*pi*f0/fs; >> h1=sin(O*n)./(n*pi); >> h1((N+1)/2)=O/pi; >> h=h1.*w; >> yh=filter(h,1,yn); >> figure,subplot(321); >> plot(tn,data);title('原信号时域图'); >> subplot(322);plot(f,F1);title('原信号频域图'); >> subplot(323);plot(tn,yn);title('加噪后时域图'); >> subplot(324);plot(f,F2);title('加噪后频域图'); >> subplot(325);plot(tn,yh);title('滤波后时域图'); >> F3=abs(fft(yh,L)); >> subplot(326);plot(f,F3);title('滤波后频域图');