数字信号实习报告.doc

数字信号上机实习报告 学院：机电学院 班级：071103-27 姓名：高玉环 指导老师：王晓莉 专题一一、实验内容设线性时不变（LTI）系统的冲激响应为h(n)，输入序列为x(n)1、h(n)=(0.8)n，0n4； x(n)=u(n)-u(n-4)2、h(n)=(0.8)nu(n), x(n)=u(n)-u(n-4)3、h(n)=(0.8)nu(n), x(n)=u(n)求以上三种情况下系统的输出y(n)。二、实验目的1、掌握离散卷积计算机实现。2、进一步对离散信号卷积算法的理解。三、原理及算法概要算法：把冲激响应h(n)与输入序列x(n)分别输入到程序中，然后调用离散卷积函数y=conv(x.,h)即可得到所要求的结果。原理：离散卷积定义为 当序列为有限长时，则 四程序x1=1 1 1 1 ;nx1=0:3;h1=1 0.8 0.64 0.83 0.84;nh1=0:4;y1=conv(x1,h1);subplot(3,3,1);stem(nx1,x1);title(序列x1);xlabel(n);ylabel(x1(n);subplot(3,3,2);stem(nh1,h1);title(序列h1);xlabel(n);ylabel(h1(n);subplot(3,3,3);stem(y1);title(序列y1);xlabel(n);ylabel(y1(n);x2=1 1 1 1;nx2=0:3;nh2=0:1:20;h2=(0.8).nh2;y2=conv(x2,h2);subplot(3,3,4);stem(x2);title(序列x2);xlabel(n);ylabel(x2(n);subplot(3,3,5);stem(h2);title(序列h2);xlabel(n);ylabel(h2(n);subplot(3,3,6);stem(y2);title(序列y2);xlabel(n);ylabel(y2(n)nx3=0:1:20;x3=1.nx3;nh3=0:1:20;h3=(0.8).nh3;y3=conv(x3,h3);subplot(3,3,7);stem(nx3,x3);title(序列x3);xlabel(n);ylabel(x3(n);subplot(3,3,8);stem(nh3,h3);title(序列h3);xlabel(n);ylabel(h3(n);subplot(3,3,9);stem(y3);title(序列y3);xlabel(n);ylabel(y3(n)六程序运行结果六结果分析有限长序列的离散卷积计算结果与理论值一致，而存在无限长序列做卷积时，由于在程序处理时是用比较长有限长序列代替的，所以与理论值基本相同。专题二、 离散傅里叶变换及其应用一、 实验内容设有离散序列 x(n)=cos(0.48n)+cos(0.52n)分析下列三种情况下的幅频特性。（1） 采集数据长度N=16，分析16点的频谱，并画出幅频特性。（2） 采集数据长度N=16，并补零到64点，分析其频谱，并画出幅频特性。（3） 采集数据长度N=64，分析46点的频谱，并画出幅频特性。观察三幅不同的幅频特性图，分析和比较它们的特点及形成原因。二、实验目的1、了解DFT及FFT的性质和特点2、利用FFT算法计算信号的频谱。三、关键算法 读入离散序列x(n)=cos(0.48n)+cos(0.52n)，采集长度为N=16的数据，调用matlab中的函数fft（x,16）与fft（x,64）对其作离散傅里叶变换得到16点、64点的频谱。采集数据长度为N=64，调用matlab中的函数fft（x,46）对其作离散傅里叶变换得到46点的频谱。原理概要：当抽样数N=2M时，以下为算法蝶形图。一般规律如下：1、 当N=2M时，则要进行M次分解，即进行M级蝶形单元的计算2、按自然顺序输入，输出是码位倒置。3、每一级包含N/2个基本蝶形运算4、第L级有2L-1个蝶群，蝶群间隔为N/2L-1如果是Matlab实现的话，可用以下两种方法计算信号频谱1、调用库函数为：fft（），直接计算X（k）2、进行矩阵运算四程序n=0:1:15;x1=cos(0.48*3.14*n)+cos(0.52*3.14*n);g1=abs(fft(x1,16);subplot(3,2,1);stem(x1);title(x1);subplot(3,2,2);stem(g1);title(g1);n2=0:1:15;x2=cos(0.48*3.14*n2)+cos(0.52*3.14*n2);x2=x2 zeros(1,48);g2=abs(fft(x2,64);subplot(3,2,3);stem(x2);title(x2);subplot(3,2,4);stem(g2);title(g2);n3=0:1:64;x3=cos(0.48*3.14*n3)+cos(0.52*3.14*n3);g3=abs(fft(x3,46);subplot(3,2,5);stem(x3);title(x3);subplot(3,2,6);stem(g3);title(g3);五程序运行结果六结果分析N 点DFT的频谱分辨率是2 /N。一节指出可以通过补零观察到更多的频点，但是这并不意味着补零能够提高真正的频谱分辨率。这是因为xn 实际上是x(t) 采样的主值序列，而将xn补零得到的xn 周期延拓之后与原来的序列并不相同，也不是x(t) 的采样。因此是不同离散信号的频谱。对于补零至M点的x的DFT，只能说它的分辨率2 /M仅具有计算上的意义，并不是真正的、物理意义上的频谱。频谱分辨率的提高只能通过提高采样频率实现。专题三 IIR滤波器的设计一、实验内容1、设计一个Butterworth数字低通滤波器，设计指标如下： 通带截止频率：0.2，幅度衰减不大于1分贝 阻带截止频率：0.3，幅度衰减大于15分贝2、让不同频率的正弦波通过滤波器，验证滤波器性能。3、分析不同滤波器的特点和结果。4、编程设计实现IIR滤波器。二、实验目的掌握不同IIR滤波器的性质、特点。通过实验学习如何设计各种常用的IIR滤波器，以便在实际工作中能根据具体情况使用IIR滤波器。三、算法算法：输入通带截止频率Wp，阻带截止频率Ws，通带衰减Rp，阻带衰减Rs，通过这些数值调用N Wn=buttord(Wp,Ws,Rp,Rs) 函数计算巴特沃斯数字滤波器的阶数N和截止频率wn，再根据阶数N通过函数b,a=butter(N,Wn)，即可得到所要的巴特沃斯滤波器。设计一个正弦波信号，再调用函数A=filter(b,a,I)让正弦波信号通过滤波器，得到滤波信号。四程序Wp=0.2;Ws=0.3;Rp=1;Rs=15;N Wn=buttord(Wp,Ws,Rp,Rs) %用于计算巴特沃斯数字滤波器的阶数N和截止频率wnb,a=butter(N,Wn); %计算N阶巴特沃斯数字滤波器系统函数分子、分母多项式的系数向量b、a,设计所需的低通滤波器h,omega=freqz(b,a,512);%返回量h包含了离散系统频响 ，调用中若N默认，默认值为512。plot(omega/pi,20*log10(abs(h);grid;xlabel(omega/pi);ylabel(Gain,dB);title(IIR Butterworth Lowpass Filter);Wp=0.2;Ws=0.3;Rp=1;Rs=15;N1,Wn1=buttord(Wp,Ws,Rp,Rs);%用于确定阶次b,a=butter(N,Wn);%用于直接设计巴特沃兹数字滤波器，即为IIR滤波器%freqz(b,a);t=1:300I=sin(0.1*pi*t)+sin(0.4*pi*t);%设计正弦波plot(I);figure;A=filter(b,a,I);%正弦波通过滤波器plot(A);五程序运行结果N = 6Wn = 0.2329N1 = 6Wn1 = 0.2329六结果分析Butterworth滤波器在通带内的频率特性是平坦的,并且随着频率的增加而衰减。正弦信号在经过IIR滤波器滤波后，高频信号被滤除，低频信号被保留了下来。专题四 用窗函数设计FIR 滤波器一、实验内容选取合适窗函数设计一个线性相位FIR低通滤波器，使它满足如下性能指标：通带截止频率：p=0.5,通带截止频率处的衰减不大于3分贝；阻带截止频率：s=0.66，阻带衰减不小于40分贝。二、实验目的1、掌握用窗函数法设计FIR滤波器的原理和方法。2、熟悉线性相位滤波器特性。3、了解各种窗函数对滤波器特性的影响。三、算法算法:通过其通带截止频率p与阻带截止频率s算出其过渡带的宽度与滤波器的长度，从而得到理想滤波器的截止频率，根据所要求的理想滤波器，得到hd（n）。由于其通带截止频率处的衰减不大于3分贝与阻带衰减不小于40分贝，我选择最接近的汉宁窗，最后调用函数h=hd.*win 及freqz(h,1,512)得到实际汉宁窗的响应和实际滤波器的幅度响应。四程序wp =0.5*pi; ws=0.66*pi; wdelta =ws-wp; %过渡带宽度 N=ceil(8*pi/wdelta) %滤波器长度 if rem(N,2)=0 N=N+1; end Nw =N; wc =(wp+ws)/2; %理想低通滤波器的截止频率 n =0: N-1; alpha =(N-1)/2; m =n-alpha+0.00001; hd =sin(wc*m)./(pi*m); %一个响应 win =(hanning(Nw); %汉宁窗 h=hd.*win; %实际汉宁窗的响应 freqz(h,1,512); %实际滤波器的幅度响应五程序运行结果N = 50五、综合一、实验内容录制一段自己的语音信号，时间为10s左右，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法或双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的语音信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，用MATLAB设计一信号处理系统界面。二算法调用函数function pushbutton1_Callback(hObject, eventdata, handles)实现一个信号处理系统界面。选择左键时，用双线性变换法设计滤波器来对信号进行处理，选择右键时，用窗函数法设计滤波器来对信号进行处理。读取语音信号，对语音信号进行f=8000的频率进行采样，调用函数y1=fft(x1,2048)对所采集的点做2048点FFT变换。先设计butterworth模拟滤波器，再用双线性变换法实现模拟滤波器到数字滤波器的转换。最后调用函数f1=filter(bz,az,x2)对加了噪声的语音信号进行滤波，得到滤波后的频谱图。3 程序function pushbutton1_Callback(hObject, eventdata, handles)% hObject handle to pushbutton1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)global f1global fsglobal x1global x2global wp1global ws1%fs=40000; Fs=fs;%语音信号采样频率为8000%x1=wavread(E:LrC专业相关课程专业课数字信号处理2013实习新建文件夹刘若英 - 为爱痴狂.wav);t=(0:length(x1)-1)/8000;y1=fft(x1,2048); %对信号做2048点FFT变换f=fs*(0:1023)/2048;%figure(1);axes(handles.axes1);%subplot(2,2,1);plot(t,x1); %做原始信号的时域波形grid on;axis tight;title(原始语音信号);xlabel(time(s);ylabel(幅度);%subplot(2,2,2);axes(handles.axes2);plot(f,abs(y1(1:1024) %做原始信号的FFT频谱grid on;axis tight;title(原始语音信号的FFT频谱);xlabel(Hz);ylabel(幅度);%双线性变换法设计的巴特沃斯滤波器A1=0.05;A2=0.10;d=A1*cos(2*pi*6000*t)+A2*sin(2*pi*8000*t);%x2=x1;n=size(d);x2(1:n,1)=x1(1:n,1)+d(1:n,1);%x2(1:n,2)=x1(1:n,2);%wp=0.10*pi;%ws=0.20*pi;wp=wp1*pi;ws=ws1*pi;Rp=1;Rs=15;Ts=1/Fs;wp1=2/Ts*tan(wp/2); %将模拟指标转换为数字指标ws1=2/Ts*tan(ws/2);N,Wn=buttord(wp1,ws1,Rp,Rs,s); %选择滤波器最小阶数Z,P,K=buttap(N); %创建butterworth模拟滤波器Bap,Aap=zp2tf(Z,P,K);b,a=lp2lp(Bap,Aap,Wn); bz,az=bilinear(b,a,Fs); %用双线性法实现模拟到数字的转换 H,W=freqz(bz,az); %绘制频率响应曲线axes(handles.axes3);%subplot(2,2,3);plot(W*Fs/(2*pi),20*log10(abs(H);grid on;axis tight;title(Butterworth滤波);xlabel(频率（Hz）);ylabel(频率响应); f1=filter(bz,az,x2);%figure(2);axes(handles.axes4);%subplot(2,2,1);plot(t,x2); %画出滤波前的时域图grid on;axis tight;title(滤波前的时域波形);axes(handles.axes5);%subplot(2,2,2);plot(t,f1); %画出滤波后的时域图grid on;axis tight;title(滤波后的时域波形);y3=fft(f1,2048);y2=fft(x2,2048);axes(handles.axes6);%subplot(2,2,3);plot(f,abs(y2(1:1024); %画出滤波前的频谱图grid on;axis tight;title(滤波前的频谱);xlabel(Hz);ylabel(幅度);axes(handles.axes7);%subplot(2,2,4);plot(f,abs(y3(1:1024); %画出滤波后的频谱图grid on;axis tight;title(滤波后的频谱);xlabel(Hz);ylabel(幅度); % - Executes on button press in pushbutton2.function pushbutton2_Callback(hObject, eventdata, handles)% hObject handle to pushbutton2 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)global f2global fsglobal x2global x1global wp2global ws2t=(0:length(x1)-1)/8000;y1=fft(x1,2048); %对信号做2048点FFT变换f=fs*(0:1023)/2048;Fs=fs;%figure(1);axes(handles.axes1);plot(t,x1); %做原始信号的时域波形grid on;axis tight;title(原始语音信号);xlabel(time(s);ylabel(幅度);axes(handles.axes2);plot(f,abs(y1(1:1024); %做原始信号的FFT频谱grid on;axis tight;title(原始语音信号的FFT频谱);xlabel(Hz);ylabel(幅度);%窗函数设计滤波器t=(0:length(x1)-1)/40000; A1=0.05;A2=0.10;d=A1*cos(2*pi*8000*t)+A2*sin(2*pi*10000*t);%x2=x1;n=size(d);x2(1:n,1)=x1(1:n,1)+d(1:n,1);%wp=wp2;ws=ws2;wdelta=ws-wp;N=ceil(8*pi/wdelta); %取整wn=(wp2+ws2)*pi/2;bz,az=fir1(N,wn/pi,hamming(N+1); %选择窗函数，并归一化截止频率H,W=freqz(bz,az); %绘制频率响应曲线axes(handles.axes3);%subplot(2,2,3);plot(W*Fs/(2*pi),20*log10(abs(H);grid on;axis tight;title(Hanming滤波);xlabel(频率（Hz）);ylabel(频率响应);%f2=filter(bz,az,x2);f2=filter(bz,az,x2);%figure(3);axes(handles.axes4);plot(t,x2);grid on;axis tight;title(滤波前的时域波形);axes(handles.axes5);plot(t,f2);grid on;axis tight;title(滤波后的时域波形);f=fs*(0:204