数字信号处理课程设计-地震信号的干扰去除.doc

DONGFANG COLLEGE，FUJIAN AGRICULTURE AND FORESTRY UNIVERSITY课程名称：数字信号处理课设题目：地震信号的干扰去除系 别：计算机系年级专业：10电子信息工程2班学 号：10503020XX姓名：凯爷任课教师：罗志聪成绩：2013年月日一、作品简介现有 1 年的某地震信号，一共有 365 个文件，每个文件长度为 10M个采样点，采样率为 200Hz，由于受到 50hz 的信号干扰，请设计一 FIR 滤波器，对其进行滤波，为了加快分析的速度，请采用学过的分段方法进行分段分析。 如果信号截至带宽为 40hz，请设计一低通滤波器 FIR 滤波器对以上信号进行处理；可自行产生 10M 采样点数的复合要求的信号，放在某文档中即可，为了证明可以对 365 个文档中的任意一个文档可以处理，请至少产生 3 个文档。二、设计目的数字信号处理是一门理论和技术发展十分迅速、应用非常广泛的前沿性学科，在雷达、通信、医学、地震等众多领域都有广泛的应用。 数字信号处理课程是电子、通信、信号处理等专业的专业基础课程，数字信号处理课程设计是完成数字信号处理相关理论的学习后进行的综合性训练课程，使学生能够对信号的采集、处理、传输、显示和存储等理论及技术实现有一个系统的掌握和理解。 本课程设计的目的为： 1、巩固和加深学生对数字信号处理的基本概念、基本原理、基本分析方法的理解； 2、培养学生根据设计题目的要求，综合利用数字信号处理技术的能力，提高学生综合运用所学知识，解决实际工程技术问题的能力； 3、通过课程设计，培养学生查阅文献、资料，独立解决实际问题的能力。 通过课程设计，要求学生综合运用数字信号处理的理论知识进行频谱分析和滤波器设计，通过理论推导得出相应结论，再利用 MATLAB 或LabView 作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念，使学生基本掌握以下内容： 1、学会 MATLAB 或 LabView 的使用，掌握其程序设计方法。 2、掌握在 Windows 环境下语音信号采集的方法。 3、掌握数字信号处理的基本概念、基本理论和基本方法。 4、掌握设计 FIR 和 IIR 数字滤波器的方法。 5、学会用 MATLAB 或 LabView 对信号进行分析和处理。三、设计技术和方法截至频率为40hz的地震信号可以由正弦信号的线性相加得到，然后利用matlab中计算出的和得到文本文档保存，即可得到“保存的地震信号”。对工频信号的滤除可以用陷波滤波器。而陷波滤波器又分自适应、小波变换、模板匹配滤波器.同时实现FIR设计陷波滤波器又有窗函数法、切比雪夫逼近法、频率抽样法。本设计采用窗函数法，窗函数简单实用，且在设计中要先构造出低通滤波器在后续中也可能用到。低通滤波器可以将40hz以上的频率滤除。为了比较各个滤波器的优劣，本设计采用了三个典型的滤波器。其中观察傅里叶变换所得频域的结果更能明显看出实验结果。傅利叶变换原理 卷积 由于信号的时宽带宽积为一常数，所以主瓣宽度与窗口宽度成反比，N越大，的主瓣越窄。同时窗的形状也对傅里叶频谱有影响，比如吉布斯效应。虽然频率分辨率很高但由于第一旁瓣的衰减很小但是求得的与的偏差较大。所以対窗的选择要根据实际需要选择合适的窗及其大小形状。而在设计中对大量信号的处理为了加快处理熟读我们可以用分段卷积：分段卷积可以分为重叠相加法和重叠保留法。重叠相加法：设h（n）的长度为M，X（n）的长度为N，将长序列X按如下方法分： 根据卷积分配率： 因此分段后，h与x 的卷积可以认为是x的各个子段与h的卷积之和。 function y=ovrlpadd(x,h,Nfft) lenx=length(x); M=length(h); M1=M-1; N=Nfft-M+1; x=x,zeros(1,N-1); K=floor(Lenx+Nfft-1); for i=0:1:K ix=i*N; xseg=x(ix+1:ix+N); yseg=circonv(xseg,h,Nfft); y(ix+1:ix+Nfft)=y(ix+1:ix+Nfft)+yseg(1:Nfft); y=y(1:Lenx+M-1);重叠保留法：同样将序列X（n）按如下方式分成长度为L的字段： 即每个子段之间有M-1哥点重叠。将分段后的x与h卷积之和即为结果。 四、课程设计的步骤和结果1.信号的产生根据傅里叶变换，任何函数都可以写作正弦函数的叠加，所以地震信号可以用频率不同的正弦信号叠加得到。在matlab的workplace区有使用的所有变量的值，将所需要的值复制到文档里，比如txt文档，然后保存。由于变量在matlab中以矩阵的形式存在，所以保存在文档里的数据也是矩阵。Fs=str2double(get(handles.fs,string); 时域信号 频域信号陷波滤波器的设计由于要求的地震信号的截至频率为40HZ，而在实际地震的信号的获取时，工频信号的存在在所难免因此要设计一个陷波滤波器滤除50HZ的信号。利用matlab设计FIR滤波器的方法有窗函数法、频率抽样法和切比雪夫逼近法等，但窗函数法我们比较熟悉且简单实用。窗函数思想：根据要求选择适当的理想低通滤波器，由于理想低通滤波器得脉冲响应是非因果且无限长，用最优化窗函数来截取就可以达到目标。 function hd=ideal(Wc,M) a=(M-1)/2; n=0:(M-1); m=n-a+0.000001; hd=sin(Wc*m)./(pi*m); 以上为一个低通滤波器，Wc为归一化截止角频率，M为阶数。陷波滤波器可以看做是一个低通滤波器加上一个高通滤波器，而高通滤波器有可以由全通滤波器减去一个低通滤波器得到。 去除工频信号后的波形低通滤波器的设计FIR滤波器的设计问题实际上就是用要求设计的FIR滤波器的频率响应H（w）去逼近所要求的理想滤波器的响应Hd（w）。从单位取样序列来看，就是使所涉及的滤波器的h（n）逼近单位取样响应序列hd（n）。既 由于理想低通滤波器 ，其单位取样响应hd（n）是无限长的，是非因果序列。因此为了构造有限长度的线性滤波器需要将理想滤波器截取一段，这可以用理想低通滤波器与一个窗函数相乘得到。窗函数可以有矩形窗、汉宁窗、海明窗、布莱克曼窗等。窗函数的要求为： A.窗谱主瓣尽可能窄，可以获取较陡的过渡带 B.尽量减小窗谱的最大旁瓣的相对幅度；即使能量集中于主瓣，使肩峰和波纹减小，增大阻带的衰减。 本设计采用以下三种： 凯赛窗： rp=1;rs=-100;fcuts=39.8 40.2;rp为通带衰减，rs为阻带衰减，fcuts为过渡带宽度。而对于凯赛窗可以直接调用matlab中的函数。n,wn,beta,ftype=kaiserord(fcuts,mags,devs,Fs);hh=fir1(n,wn,ftype,kaiser(n+1,beta);Y=fft(p1);x=fftfilt(hh,p1); 凯赛窗频率响应 经过低通滤波器后的频谱布莱克曼：fp1=39; N=2*ceil(12*pi/(ws1-wp1)/2); window=blackman(N+1); wn=(wp1+ws1)/2; bf1=fir1(N,wn/pi,low,window); output1=fftfilt(bf1,p1); X=fft(output1); 布莱克曼窗函数频率响应 汉宁窗 bf2=fir1(n-1,wc,low,hanning(n); output2=fftfilt(bf2,p1);