西安电子科技大学数字信号处理大作业

数字信号处理MATLAB上机作业M 2.21.题目The square wave and the sawtooth wave are two periodic sequences as sketched in figure ing the function stem. The input data specified by the user are: desired length L of the sequence, peak value A, and the period N. For the square wave sequence an additional user-specified parameter is the duty cycle, which is the percent of the period for which the signal is positive. Using this program generate the first 100 samples of each of the above sequences with a sampling rate of 20 kHz ,a peak value of 7, a period of 13 ,and a duty cycle of 60% for the square wave.2.程序% 用户定义各项参数参数A = input('The peak value =');L = input('Length of sequence =');N = input('The period of sequence =');FT = input('The desired sampling frequency =');DC = input('The square wave duty cycle = ');% 产生所需要的信号t = 0:L-1;T = 1/FT;x = A*sawtooth(2*pi*t/N);y = A*square(2*pi*(t/N),DC);% Plotsubplot(2,1,1)stem(t,x);ylabel('幅度');xlabel('n');subplot(2,1,2)stem(t,y);ylabel('幅度');xlabel('n');3.结果4.结果分析M 2.41.题目(a)Write a matlab program to generate a sinusoidal sequence x[n]= Acos(ω0 n+Ф) and plot thesequence using the stem function. The input data specified by the user are the desired length L, amplitude A, the angular frequency ω0 , and the phase Фwhere 0<ω0 <pi and 0<=Ф<=2pi. Using this program generate the sinusoidal sequences shown in figure 2.15. (b)Generate sinusoidal sequences with the angular frequencies given in Problem 2.22.Determine the period of each sequence from the plot and verify the result theoretically. 2.程序%用户定义的参数L = input('Desired length = ');A = input('Amplitude = ');omega = input('Angular frequency = ');phi = input('Phase = ');%信号产生n = 0:L-1;x = A*cos(omega*n + phi);stem(n,x);xlabel('n');ylabel('幅度');title(['\omega_{o} = ',num2str(omega)]);3.结果(a)ω0=0ω0=0.1πω0=0.8πω0=1.2π(b)ω0=0.14πω0=0.24πω0=0.34πω0=0.68πω0=0.75π4.结果分析M 2.51.题目Generate the sequences of problem 2.21(b) to 2.21(e) using matlab.2.程序(b)n = 0 : 99;x=sin(0.6*pi*n+0.6*pi);stem(n,x);xlabel('n');ylabel('幅度');(c)n = 0 : 99;x=2*cos(1.1*pi*n-0.5*pi)+2*sin(0.7*pi*n);stem(n,x);xlabel('n');ylabel('幅度');(d)n = 0 : 99;x=3*sin(1.3*pi*n-4*cos(0.3*pi*n+0.45*pi));stem(n,x);xlabel('n');ylabel('幅度');(e)n = 0 : 99;x=5*sin(1.2*pi*n+0.65*pi)+4*sin(0.8*pi*n)-cos(0.8*pi*n);stem(n,x);xlabel('n');ylabel('幅度');(f)n = 0 : 99;x=mod(n,6);stem(n,x);xlabel('n');ylabel('幅度');3.结果(b)(c)(d)(e)(f)4.结果分析M 2.61.题目Write a matlab program to plot a continuous-time sinusoidal signal and its sampled version and verify figure 2.19. You need to use the hold function to keep both plots.2.程序%用户定义的参数fo = input('Frequency of sinusoid in Hz = ');FT = input('Samplig frequency in Hz = ');%产生信号t = 0:0.001:1;g1 = cos(2*pi*fo*t);plot(t,g1,'-')xlabel('时间t');ylabel('幅度')holdn = 0:1:FT;gs = cos(2*pi*fo*n/FT);plot(n/FT,gs,'o');hold off3.结果4.结果分析M 3.11.题目Using program 3_1 determine and plot the real and imaginary parts and the magnitude and phase spectra of the following DTFT for various values of r and θ:G(e jω)=1, 0<r<1.1−2r(cosθ)e−jω+r2e−2jω2.程序%program 3_1%discrete-time fourier transform computatition%k=input('Number of frequency points = ');num=input('Numerator coefficients= ');den=input('Denominator coefficients= ');%computer the frequency responsew=0:pi/k:pi;h=freqz(num,den,w);%plot the frequency responsesubplot(221)plot(w/pi,real(h));gridtitle('real part')xlabel('\omega/\pi');ylabel('Amplitude') subplot(222)plot(w/pi,imag(h));gridtitle('imaginary part')xlabel('\omega/\pi');ylabel('Amplitude') subplot(223)plot(w/pi,abs(h));gridtitle('magnitude spectrum')xlabel('\omega/\pi');ylabel('magnitude') subplot(224)plot(w/pi,angle(h));gridtitle('phase spectrum')xlabel('\omega/\pi');ylabel('phase,radians')3.结果(a)r=0.8 θ=π/6(b)r=0.6 θ=π/34.结果分析M 3.41.题目Using matlab verify the following general properties of the DTFT as listed in Table 3.2:(a) Linearity, (b) time-shifting, (c) frequency-shifting, (d) differentiation-in-frequency, (e) convolution, (f) modulation, and (g) Parseval’s relation. Since all data in matlab have to be finite-length vectors, the sequences to be used to verify the properties are thus restricted to be of finite length.2.程序%先定义两个信号N = input('The length of the sequence = ');k = 0:N-1;%g为正弦信号g = 2*sin(2*pi*k/(N/2));%h为余弦信号h = 3*cos(2*pi*k/(N/2));[G,w] = freqz(g,1);[H,w] = freqz(h,1);%*************************************************************************%% 线性性质alpha = 0.5;beta = 0.25;y = alpha*g+beta*h;[Y,w] = freqz(y,1);figure(1);subplot(211),plot(w/pi,abs(Y));xlabel('\omega/\pi');ylabel('|Y(e^j^\omega)|');title('线性叠加后的频率特性');grid;% 画出Y 的频率特性subplot(212),plot(w/pi,alpha*abs(G)+beta*abs(H));xlabel('\omega/\pi');ylabel('\alpha|G(e^j^\omega)|+\beta|H(e^j^\omega)|');title('线性叠加前的频率特性');grid;% 画出alpha*G+beta*H 的频率特性%*************************************************************************% % 时移性质n0 = 10;%时移10个的单位y2 = [zeros([1,n0]) g];[Y2,w] = freqz(y2,1);G0 = exp(-j*w*n0).*G;figure(2);subplot(211),plot(w/pi,abs(G0));xlabel('\omega/\pi');ylabel('|G0(e^j^\omega)|');title('G0的频率特性');grid;% 画出G0的频率特性subplot(212),plot(w/pi,abs(Y2));xlabel('\omega/\pi');ylabel('|Y2(e^j^\omega)|');title('Y2的频率特性');grid;% 画出Y2 的频率特性%*************************************************************************% % 频移特性w0 = pi/2; % 频移pi/2r=256; %the value of w0 in terms of number of samplesk = 0:N-1;y3 = g.*exp(j*w0*k);[Y3,w] = freqz(y3,1);% 对采样的512个点分别进行减少pi/2,从而生成G(exp(w-w0))k = 0:511;w = -w0+pi*k/512;G1 = freqz(g,1,w);figure(3);subplot(211),plot(w/pi,abs(Y3));xlabel('\omega/\pi');ylabel('|Y3(e^j^\omega)|');title('Y3的频率特性');grid;% 画出Y3的频率特性subplot(212),plot(w/pi,abs(G1));xlabel('\omega/\pi');ylabel('|G1(e^j^\omega)|');title('G1的频率特性');grid;% 画出G1 的频率特性%*************************************************************************% % 频域微分k = 0:N-1;y4 = k.*g;[Y4,w] = freqz(y4,1);%在频域进行微分y0 = ((-1).^k).*g;G2 = [G(2:512)' sum(y0)]';delG = (G2-G)*512/pi;figure(4);subplot(211),plot(w/pi,abs(Y4));xlabel('\omega/\pi');ylabel('|Y4(e^j^\omega)|');title('Y4的频率特性');grid;% 画出Y4的频率特性subplot(212),plot(w/pi,abs(delG));xlabel('\omega/\pi');ylabel('|delG(e^j^\omega)|');title('delG的频率特性');grid;% 画出delG的频率特性%*************************************************************************% % 相乘性质y5 = conv(g,h);%时域卷积[Y5,w] = freqz(y5,1);figure(5);subplot(211),plot(w/pi,abs(Y5));xlabel('\omega/\pi');ylabel('|Y5(e^j^\omega)|');title('Y5的频率特性');grid;% 画出Y5的频率特性subplot(212),plot(w/pi,abs(G.*H));%频域乘积xlabel('\omega/\pi');ylabel('|G.*H(e^j^\omega)|');title('G.*H的频率特性');grid;% 画出G.*H的频率特性%*************************************************************************% % 帕斯瓦尔定理y6 = g.*h;%对于freqz函数，在0到2pi直接取样[Y6,w] = freqz(y6,1,512,'whole');[G0,w] = freqz(g,1,512,'whole');[H0,w] = freqz(h,1,512,'whole');% Evaluate the sample value at w = pi/2% and verify with Y6 at pi/2H1 = [fliplr(H0(1:129)') fliplr(H0(130:512)')]';val = 1/(512)*sum(G0.*H1);% Compare val with Y6(129) i.e sample at pi/2 % Can extend this to other points similarly% Parsevals theoremval1 = sum(g.*conj(h));val2 = sum(G0.*conj(H0))/512;% Comapre val1 with val23.结果(a)(b)(c)(d)(e)4.结果分析M 3.81.题目Using matlab compute the N-point DFTs of the length-N sequences of Problem 3.12 for N=3, 5, 7, and 10. Compare your results with that obtained by evaluating the DTFTs computed in Problem 3.12 at ω= 2pik/N, k=0, 1,……N-1.2.程序%用户定义N的长度N = input('The value of N = ');k = -N:N;y1 = ones([1,2*N+1]);w = 0:2*pi/255:2*pi;Y1 = freqz(y1, 1, w);%对y1做傅里叶变换Y1dft = fft(y1);k = 0:1:2*N;plot(w/pi,abs(Y1),k*2/(2*N+1),abs(Y1dft),'o');grid;xlabel('归一化频率');ylabel('幅度');(a)clf;N = input('The value of N = ');k = -N:N;y1 = ones([1,2*N+1]);w = 0:2*pi/255:2*pi;Y1 = freqz(y1, 1, w);Y1dft = fft(y1);k = 0:1:2*N;plot(w/pi,abs(Y1),k*2/(2*N+1),abs(Y1dft),'o');xlabel('Normalized frequency');ylabel('Amplitude');(b)%用户定义N的长度N = input('The value of N = ');k = -N:N;y1 = ones([1,2*N+1]);y2 = y1 - abs(k)/N;w = 0:2*pi/255:2*pi;Y2 = freqz(y2, 1, w);%对y1做傅里叶变换Y2dft = fft(y2);k = 0:1:2*N;plot(w/pi,abs(Y2),k*2/(2*N+1),abs(Y2dft),'o');grid;xlabel('归一化频率');ylabel('幅度');(c)%用户定义N的长度N = input('The value of N = ');k = -N:N;y3 =cos(pi*k/(2*N));w = 0:2*pi/255:2*pi;Y3 = freqz(y3, 1, w);%对y1做傅里叶变换Y3dft = fft(y3);k = 0:1:2*N;plot(w/pi,abs(Y3),k*2/(2*N+1),abs(Y3dft),'o');grid;xlabel('归一化频率');ylabel('幅度');3.结果(a)N=3N=5 N=7N=10 (b)N=3N=5 N=7N=10 (c)N=3N=5 N=7N=104.结果分析M 3.191.题目Using Program 3_10 determine the z-transform as a ratio of two polynomials in z-1 from each of the partial-fraction expansions listed below:(a)X1(z)=−2+104+z−1−82+z−1,|z|>0.5,(b)X2(z)=3.5−21−0.5z−1−3+z−11−0.25z−2,|z|>0.5,(c)X3(z)=5(3+2z−1)2−43+2z−1+31+0.81z−2,|z|>0.9,(d)X4(z)=4+105+2z−1+z−16+5z−1+z−2,|z|>0.5.2.程序% Program 3_10% Partical-Fraction Expansion to rational z-Transform %r = input('Type in the residues = ');p = input('Type in the poles = ');k = input('Type in the constants = ');[num, den] = residuez(r,p,k);disp('Numberator polynominal coefficients');disp(num) disp('Denominator polynomial coefficients'); disp(den)4.结果分析M 4.61.题目Plot the magnitude and phase responses of the causal IIR digital transfer functionH(z)=0.0534(1+z−1)(1−1.0166z−1+z−2) (1−0.683z−1)(1−1.4461z−1+0.7957z−2).What type of filter does this transfer function represent? Determine the difference equation representation of the above transfer function.2.程序b=[0.0534 -0.00088644 -0.00088644 0.0534];a=[1 -2.1291 1.7833863 -0.5434631];figure(1)freqz(b,a);figure(2)[H,w]=freqz(b,a);plot(w/pi,abs(H)),grid;xlabel('Normalized Frequency (\times\pi rad/sample)'),ylabel('Magnitude');幅度化成真值之后：4.结果分析H(z)=0.0534−0.00088644z−1−0.00088644z−2+0.0534z−31−2.1291z−1+1.7833863z−2−0.5434631z−3M 4.71.题目Plot the magnitude and phase responses of the causal IIR digital transfer functionH(z)=(1−z−1)4(1−1.499z−1+0.8482z−2)(1−1.5548z−1+0.6493z−2).2.程序b=[1 -4 6 -4 1];a=[1 -3.0538 3.8227 -2.2837 0.5472]; figure(1)freqz(b,a);figure(2)[H,w]=freqz(b,a);plot(w/pi,abs(H)),grid;xlabel('Normalized Frequency (\times\pi rad/sample)'), ylabel('Magnitude');3.结果4.结果分析。

数字信号处理上机第一次实验实验一：设给定模拟信号()1000t a x t e -=，的单位是ms 。

(1) 利用MATLAB 绘制出其时域波形和频谱图（傅里叶变换），估计其等效带宽（忽略谱分量降低到峰值的3%以下的频谱）。

(2) 用两个不同的采样频率对给定的进行采样。

○1。

○2 。

比较两种采样率下的信号频谱，并解释。

实验一MATLAB 程序：（1）○1 clc; fs=5000;ts=1/fs;N=1000;t=(-N:N)*ts;s=exp(-abs(t));plot(t,s,'linewidth',1.5)xlabel('时间')ylabel('幅度')set(gca,'fontweight','b','fontsize',12)SPL=N*100;figuresp=fftshift(fft(s,SPL));sp=sp/max(sp)*100;freqb=-fs/2:fs/SPL:fs/2-fs/SPL;plot(freqb,abs(sp))xlabel('频率')ylabel('频谱幅度')set(gca,'fontweight','b','fontsize',12)yy=abs(abs(sp)-3);[aa,freqind]=min(yy);(freqind-SPL/2)*fs/SPLt ()a x t ()()15000s a f x t x n =以样本秒采样得到。

()()11j x n X e ω画出及其频谱()()11000s a f x t x n =以样本得到。

()()11j x n X e ω画出及其频谱○2 clc;fs=1000;ts=1/fs;N=1000;t=(-N:N)*ts;s=exp(-abs(t));plot(t,s,'linewidth',1.5)xlabel('时间')ylabel('幅度')set(gca,'fontweight','b','fontsize',12) SPL=N*100;figuresp=fftshift(fft(s,SPL));sp=sp/max(sp)*100;freqb=-fs/2:fs/SPL:fs/2-fs/SPL;plot(freqb,abs(sp))xlabel('频率')ylabel('频谱幅度')set(gca,'fontweight','b','fontsize',12)yy=abs(abs(sp)-3);[aa,freqind]=min(yy);(freqind-SPL/2)*fs/SPL实验三：设，，编写MATLAB 程序，计算：（１） ５点圆周卷积；（２） ６点圆周卷积；（３） 线性卷积；（４） 画出的，和时间轴对齐。

数字信号处理大作业班级：021231学号：姓名：指导老师：吕雁一写出奈奎斯特采样率和和信号稀疏采样的学习报告和体会1、采样定理在进行A/D信号的转换过程中，当采样频率fs.max大于信号中最高频率fmax的2倍时(fs.max>2fmax)，采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的5～10倍；采样定理又称奈奎斯特定理。

(1)在时域频带为F的连续信号 f(t)可用一系列离散的采样值f(t1),f(t1±Δt)，f(t1±2Δt)，...来表示,只要这些采样点的时间间隔Δt≤1/2F，便可根据各采样值完全恢复原始信号。

(2)在频域当时间信号函数f(t)的最高频率分量为fmax时,f(t)的值可由一系列采样间隔小于或等于1/2fo的采样值来确定,即采样点的重复频率fs ≥2fmax。

2、奈奎斯特采样频率(1)概述奈奎斯特采样定理：要使连续信号采样后能够不失真还原，采样频率必须大于信号最高频率的两倍（即奈奎斯特频率）。

奈奎斯特频率（Nyquist frequency）是离散信号系统采样频率的一半，因哈里·奈奎斯特（Harry Nyquist）或奈奎斯特－香农采样定理得名。

采样定理指出，只要离散系统的奈奎斯特频率高于被采样信号的最高频率或带宽，就可以真实的还原被测信号。

反之，会因为频谱混叠而不能真实还原被测信号。

采样定理指出，只要离散系统的奈奎斯特频率高于采样信号的最高频率或带宽，就可以避免混叠现象。

从理论上说，即使奈奎斯特频率恰好大于信号带宽，也足以通过信号的采样重建原信号。

但是，重建信号的过程需要以一个低通滤波器或者带通滤波器将在奈奎斯特频率之上的高频分量全部滤除，同时还要保证原信号中频率在奈奎斯特频率以下的分量不发生畸变，而这是不可能实现的。

在实际应用中，为了保证抗混叠滤波器的性能，接近奈奎斯特频率的分量在采样和信号重建的过程中可能会发生畸变。

第二次编程作业加载信号(EEGdata.txt)；分析信号的幅度谱,确定低通FIR数字滤波器的指标；分别利用各种窗函数（Rectangular, Hamming, Kaiser）设计此FIR滤波器；对信号加随机噪声，并用设计的FIR滤波器对含噪声信号进行滤波。

要求：1）给出程序，画出滤波器幅度谱及其增益图(dB)；分析设计的滤波器是否达到指标。

2）利用设计的三种滤波器对加载的信号分别进行滤波，对比分析滤波结果。

程序：data = load('EEGdata.txt');n=1:length(data);xn=0.1*rand(length(data),1)-0.05;x=xn+data;[h0,w0]=freqz(data);[h1,w1]=freqz(x);figure(1)plot(w0/pi,abs(h0),w1/pi,abs(h1),'r');legend('Ô-dataÐźÅƵÆ×ͼ','¼ÓÔëÐźÅƵÆ×ͼ');wp = 0.15*pi;ws = 0.2*pi;DB = ws-wp;N = ceil(2*3.11*pi/DB+1);wc = (wp+ws)/2;wc = wc/pi;win = rectwin(N+1);b = fir1(N, wc, 'low', win);figure(2)freqz(b,1);wp = 0.15*pi;ws = 0.2*pi;DB = ws-wp;N = ceil(2*3.11*pi/DB+1);wc = (wp+ws)/2;wc = wc/pi;win = hamming(N+1);b1 = fir1(N, wc, 'low', win);figure(3)freqz(b1,1);fpts = [0.15, 0.2];mag = [1, 0];dev = [0.01, 0.01];[N, Wn, beta, ftype] = kaiserord(fpts, mag, dev)Kwin = kaiser(N+1, beta);b2 = fir1(N, Wn, Kwin);[h,omega] = freqz(b2,1,512);figure(4)freqz(b2,1);y=filter(b, 1, x);figure(5)subplot(2,1,1);plot(n,data,n,x,'g',n,y,'r');title('¾ØÐδ°Â˲¨·ù¶ÈÆ×Ч¹û¶Ô±Èͼ');lege nd('Ô-ÐźÅdata','¼ÓÔëÐźÅx','Â˲¨ºóÐźÅy');[H,W]=freqz(y);subplot(2,1,2);plot(w0/pi,abs(h0),W/pi,abs(H),'r');title('¾ØÐδ°Â˲¨ÆµÆ×Ч¹û¶Ô±Èͼ') ;legend('Ô-ÐźÅdata','Â˲¨ºóÐźÅy');y1=filter(b1, 1, x);figure(6)subplot(2,1,1);plot(n,data,n,x,'g',n,y1,'r');title('Hamming´°Â˲¨·ù¶ÈÆ×Ч¹û¶Ô±Èͼ'); legend('Ô-ÐźÅdata','¼ÓÔëÐźÅx','Â˲¨ºóÐźÅy1');[H,W]=freqz(y1);subplot(2,1,2);plot(w0/pi,abs(h0),W/pi,abs(H),'r');title('Hamming´°Â˲¨ÆµÆ×Ч¹û¶Ô±Èͼ');legend('Ô-ÐźÅdata','Â˲¨ºóÐźÅy1');y2=filtfilt(b2, 1, x); figure(7) subplot(2,1,1);plot(n,data,n,x,'g',n,y2,'r');title('Kaiser´°Â˲¨·ù¶ÈÆ×Ч¹û¶Ô±Èͼ');l egend('Ô-ÐźÅdata','¼ÓÔëÐźÅx','Â˲¨ºóÐźÅy2'); [H,W]=freqz(y2); subplot(2,1,2);plot(w0/pi,abs(h0),W/pi,abs(H),'r');title('Kaiser´°Â˲¨ÆµÆ×Ч¹û¶Ô±Èͼ');legend('Ô-ÐźÅdata','Â˲¨ºóÐźÅy2');figure(8)plot(n,data,n,y,'c',n,y1,'r',n,y2,'g');title('ÈýÖÖ´°º¯ÊýÂ˲¨Ð§¹û¶Ô±Èͼ');legend('Ô-dataÐźÅ','¾ØÐδ°Â˲¨','Hamming´°Â˲¨','Kaiser´°Â˲¨');0.10.20.30.40.50.60.70.80.9102468101214161800.10.20.30.40.50.60.70.80.91-3000-2000-1000Normalized Frequency (⨯π rad/sample)P h a s e (d e g r e e s )00.10.20.30.40.50.60.70.80.91-100-5050Normalized Frequency (⨯π rad/sample)M a g n i t u d e (d B )00.10.20.30.40.50.60.70.80.91-3000-2000-1000Normalized Frequency (⨯π rad/sample)P h a s e (d e g r e e s )0.10.20.30.40.50.60.70.80.91-150-100-50050Normalized Frequency (⨯π rad/sample)M a g n i t u d e (d B )0.10.20.30.40.50.60.70.80.91-2000-1500-1000-500Normalized Frequency (⨯π rad/sample)P h a s e (d e g r e e s )0.10.20.30.40.50.60.70.80.91-150-100-50050Normalized Frequency (⨯π rad/sample)M a g n i t u d e (d B )-0.2-0.100.10.2矩形窗滤波幅度谱效果对比图0.10.20.30.40.50.60.70.80.91051015矩形窗滤波频谱效果对比图-0.2-0.100.10.2Hamming 窗滤波幅度谱效果对比图0.10.20.30.40.50.60.70.80.91051015Hamming 窗滤波频谱效果对比图-0.2-0.100.10.2Kaiser 窗滤波幅度谱效果对比图0.10.20.30.40.50.60.70.80.9105101520Kaiser 窗滤波频谱效果对比图结果分析：滤波效果皆不尽如人意，原因是加载的噪声信号是随机信号，各频率皆有，滤波器通带部分的噪声无法滤掉。

数字信号处理大作业院系：电子工程学院学号：02115043姓名：邱道森实验一：信号、系统及系统响应一、实验目的(1) 熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。

(2) 熟悉时域离散系统的时域特性。

(3) 利用卷积方法观察分析系统的时域特性。

(4) 掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。

二、实验原理采样是连续信号数字处理的第一个关键环节。

对连续信号()a x t 进行理想采样的过程可用(1.1)式表示：()()()ˆa a xt x t p t =⋅ 其中()t xa ˆ为()a x t 的理想采样，()p t 为周期冲激脉冲，即 ()()n p t t nT δ∞=-∞=-∑()t xa ˆ的傅里叶变换()j a X Ω为 ()()s 1ˆj j j a a m X ΩX ΩkΩT ∞=-∞=-∑进行傅里叶变换，()()()j ˆj e d Ωt a a n X Ωx t t nT t δ∞∞--∞=-∞⎡⎤=-⎢⎥⎣⎦∑⎰ ()()j e d Ωtan x t t nT t δ∞∞--∞=-∞=-∑⎰()j e ΩnTan x nT ∞-=-∞=∑式中的()a x nT 就是采样后得到的序列()x n ， 即()()a x n x nT =()x n 的傅里叶变换为()()j j e enn X x n ωω∞-=-∞=∑比较可知()()j ˆj e aΩTX ΩX ωω==为了在数字计算机上观察分析各种序列的频域特性，通常对()j e X ω在[]0,2π上进行M 点采样来观察分析。

对长度为N 的有限长序列()x n ，有()()1j j 0eekk N nn X x n ωω--==∑其中2π,0,1,,1k k k M Mω==⋅⋅⋅-一个时域离散线性时不变系统的输入/输出关系为()()()()()m y n x n h n x m h n m ∞=-∞=*=-∑上述卷积运算也可以转到频域实现()()()j j j e e e Y X H ωωω= (1.9)三、实验内容及步骤(1) 认真复习采样理论、 离散信号与系统、 线性卷积、 序列的傅里叶变换及性质等有关内容， 阅读本实验原理与方法。

DSP实验课程序设计报告学院：电子工程学院学号：1202121013姓名：赵海霞指导教师：苏涛DSP实验课大作业设计一实验目的在DSP上实现线性调频信号的脉冲压缩、动目标显示（MTI）和动目标检测(MTD)，并将结果与MATLAB上的结果进行误差仿真。

二实验内容2.1 MATLAB仿真设定带宽、脉宽、采样率、脉冲重复频率，用MATLAB产生16个脉冲的LFM，每个脉冲有4个目标（静止，低速，高速），依次做2.1.1 脉压2.1.2 相邻2脉冲做MTI，产生15个脉冲2.1.3 16个脉冲到齐后，做MTD，输出16个多普勒通道2.2 DSP实现将MATLAB产生的信号，在visual dsp中做脉压，MTI、MTD，并将结果与MATLAB作比较。

三实验原理3.1 线性调频线性调频脉冲压缩体制的发射信号其载频在脉冲宽度内按线性规律变化即用对载频进行调制(线性调频)的方法展宽发射信号的频谱，在大时宽的前提下扩展了信号的带宽。

若线性调频信号中心频率为f，脉宽为τ，带宽为B，幅度为A，μ为调频斜率，则其表达式如下：]212cos[)()(20t t f t rect A t x μπτ+••=；)(为矩形函数rect 在相参雷达中，线性调频信号可以用复数形式表示，即)]212(exp[)()(20t t f j t rect A t x μπτ+••= 在脉冲宽度内，信号的角频率由220μτπ-f 变化到220μτπ+f 。

3.2 脉冲压缩原理脉冲雷达信号发射时，脉冲宽度τ决定着雷达的发射能量，发射能量越大， 作用距离越远；在传统的脉冲雷达信号中，脉冲宽度同时还决定着信号的频率宽度B ，即带宽与时宽是一种近似倒数的关系。

脉冲越宽，频域带宽越窄，距离分辨率越低。

脉冲压缩的主要目的是为了解决信号的作用距离和信号的距离分辨率之间的矛盾。

为了提高信号的作用距离，我们就需要提高信号的发射功率，因此，必须提高发射信号的脉冲宽度，而为了提高信号的距离分辨率，又要求降低信号的脉冲宽度。

一系统响应及稳定性的实验报告一. 实验目的:（1）掌握 求系统响应的方法。

（2）掌握时域离散系统的时域特性。

（3）分析、观察及检验系统的稳定性。

二. 实验原理与方法:1.在时域中，描写系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。

已知输入信号可以由差分方程、单位脉冲响应或系统函数求出系统对于该输入信号的响应。

在计算机上可用filter 函数求差分方程的解， conv 函数计算输入信号和系统的单位脉冲响应的线性卷积，求出系统的响应。

2.系统的时域特性指的是系统的线性时不变性质、因果性和稳定性。

3.系统的稳定性是指对任意有界的输入信号，系统都能得到有界的系统响应。

或者系统的单位脉冲响应满足绝对可和的条件。

系统的稳定性由其差分方程的系数决定。

实际中检查系统是否稳定，不可能检查系统对所有有界的输入信号，输出是否都是有界输出，或者检查系统的单位脉冲响应满足绝对可和的条件。

可行的方法是在系统的输入端加入单位阶跃序列，如果系统的输出趋近一个常数（包括零），就可以断定系统是稳定的[19]。

系统的稳态输出是指当∞→n 时，系统的输出。

如果系统稳定，信号加入系统后，系统输出的开始一段称为暂态效应，随n 的加大，幅度趋于稳定，达到稳态输出。

三．实验内容及步骤:1.编制程序，包括产生输入信号、单位脉冲响应序列的子程序，用filter 函数或conv 函数求解系统输出响应的主程序。

程序中要有绘制信号波形的功能。

2.给定一个低通滤波器的差分方程为)1(9.0)1(05.0)(05.0)(-+-+=n y n x n x n y a) 分别求出系统对)()(81n R n x =和)()(2n u n x =的响应序列，并画出其波形。

b) 求出系统的单位冲响应，画出其波形。

3.给定系统的单位脉冲响应为)()(101n R n h =)3()2(5.2)1(5.2)()(2-+-+-+=n n n n n h δδδδ 用线性卷积法分别求系统h 1(n)和h 2(n)对)()(81n R n x =的输出响应，并画出波形。

西电数字信号处理⼤作业第⼆章2.25 已知线性时不变系统的差分⽅程为若系统的输⼊序列x(x)={1,2,3,4,2,1}编写利⽤递推法计算系统零状态响应的MATLAB程序，并计算出结果。

代码及运⾏结果：>> A=[1,-0.5];>> B=[1,0,2];>> n=0:5;>> xn=[1,2,3,4,2,1];>> zx=[0,0,0];zy=0;>> zi=filtic(B,A,zy,zx);>> yn=filter(B,A,xn,zi);>> figure(1)>> stem(n,yn,'.');>> grid on;2.28图所⽰系统是由四个⼦系统T1、T2、T3和T4组成的，分别⽤单位脉冲响应或差分⽅程描述为T1：其他T2:其他T3：T4:编写计算整个系统的单位脉冲响应h(n)，0≤n≤99的MATLAB程序，并计算结果。

代码及结果如下：>> a=0.25;b=0.5;c=0.25;>> ys=0;>> xn=[1,zeros(1,99)];>> B=[a,b,c];>> A=1;>> xi=filtic(B,A,ys);>> yn1=filter(B,A,xn,xi);>> h1=[1,1/2,1/4,1/8,1/16,1/32]; >> h2=[1,1,1,1,1,1];>> h3=conv(h1,h2);>> h31=[h3,zeros(1,89)]; >> yn2=yn1+h31;>> D=[1,1];C=[1,-0.9,0.81]; >> xi2=filtic(D,C,yn2,xi); >> xi2=filtic(D,C,ys);>> yn=filter(D,C,yn2,xi); >> n=0:99;>> figure(1)>> stem(n,yn,'.');>> title('单位脉冲响应'); >> xlabel('n');ylabel('yn');2.30 利⽤MATLAB画出受⾼斯噪声⼲扰的正弦信号的波形，表⽰为其中v(n)是均值为零、⽅差为1的⾼斯噪声。

Answer to “Digital Signal Processing of 2005”Problem 1(a) even part: };5.0,1,7,7,5,7,7,1,5.0{---=e X odd part: };5.0,1,3,1,0,1,3,1,5.0{----=o X(b) };20,16,11,94,36,40,31,16,12,0{-----=y (c) MATLAB Programn=-4:2;x=[1 -2 4 6 -5 8 10]; [x11,n11]=sigshift(x,n,2); [x12,n12]=sigshift(x,n,-1); [x13,n13]=sigfold(x,n); [x13,n13]=sigshift(x13,n13,-2); [x12,n12]=sigmult(x,n,x12,n12); [y,n]=sigadd(2*x11,n11,x12,n12); [y,n]=sigadd(y,n,-1*x13,n13)Problem 2(a)w j w j w j w j jw jw e e e e e e X 65424210124)(-----++++++=，()j X e ωis periodic in ω with period 2π(b) MATLAB Program :clear; close all;n = 0:6; x = [4,2,1,0,1,2,4]; w = [0:1:1000]*pi/1000;X = x*exp(-j*n'*w); magX = abs(X); phaX = angle(X); % Magnitude Response Plotsubplot(2,1,1); plot(w/pi,magX);grid;xlabel('frequency in pi units'); ylabel('|X|'); title('Magnitude Response'); % Phase response plotsubplot(2,1,2); plot(w/pi,phaX*180/pi);grid;xlabel('frequency in pi units'); ylabel('Degrees'); title('Phase Response'); axis([0,1,-180,180])(c) Because the given sequence x (n)={4,2,1,0,1,2,4} (n=0,1,2,3,4,5,6) is symmetric about 132N α-==,the phase response ()j H e ω< satisfied the condition ：()3j H e ωαωω<=-=- so the phase response is a linear function in ω.(d) 150,350Hz Hz Ω=-;(e) The difference of amplitude and magnitude response:Firstly, the amplitude response is a real function, and it may be both positive and negative. The magnitude response is always positive.Secondly, the phase response associated with the magnitude response is a discontinuous function. While the associated with the amplitude is a continuous linear function.Problem 3(a) )9.09.01/()1()(211------=z z z z HZero:0 and 1; Pole:-0.6 and 1.5; (b)1116151()212110.61 1.5H z z z--=⨯+⨯+-， 165()((0.6)(1.5))()2121n nh n u n =-+ (c) ROC : 0.6 1.5Z <<，()()()163531215212n nh n u n u n ⎛⎫⎛⎫=⨯--⨯--- ⎪ ⎪⎝⎭⎝⎭ Problem 4(a) y(n)={50,44,34,52};(b) y(n)={5,16,34,52,45,28,0}; (c) N=6;(d) MATLAB Program :Function y=circonv(x1,x2,N) If (length(x1)>N)error(“N must not be smaller than the length of sequence ”) elsex1=[x1,zeros(1,N-length(x1))]; endif(length(x2)>N)error(“N must not be smaller than the length of sequence ”)elsex2=[x2,zeros(1,N-length(x2))]; endy1=dft(x1,N).*dft(x2,N); y=idft(y,N);(e) DTFT is discrete in time domain, but continuous in frequency domain. The DFT is discrete both in time and frequency domain.The FFT is a very efficient method for calculating DFT.Problem 5(a) Direct form II uses the little delay and it can decrease the space of the compute. (b)The advantage of the linear-phase form:1. For frequency-selective filters, linear-phase structure is generally desirable to have a phase-responsethat is a linear function of frequency.2. This structure requires 50% fewer multiplications than the direct form. (c) Block diagrams are shown as under:1z -1z -1z -1z -1z -1z -()x n )n()x n1-1-Problem 6(a) we use Hamming or Blackman window to design the bandpass filter because it can provide us attenuationexceed 60dB .(b) According to Blackman window :first, Determine transition width =p s W W - ;second, Determine the type of the window according to s A ;third, Compute M using the formula MW W p s π2=- ; fourth, Compute ideal LPF2sp W W Wc +=;fifth, design the window needed, multiply point by point; sixth, determines p A R ,(c) MATLAB Program :%% Specifications about Blackman window: ws1 = 0.2*pi; % lower stopband edge wp1 = 0.3*pi; % lower passband edge wp2 = 0.6*pi; % upper passband edge ws2 = 0.7*pi; % upper stopband edge Rp = 0.5; % passband ripple As = 60; % stopband attenuation %tr_width = min((wp1-ws1),(ws2-wp2));M = ceil(6.6*pi/tr_width); M = 2*floor(M/2)+1, % choose odd M n = 0:M-1;w_ham = (hamming(M))';wc1 = (ws1+wp1)/2; wc2 = (ws2+wp2)/2; hd = ideal_lp(wc2,M)-ideal_lp(wc1,M); h = hd .* w_ham;[db,mag,pha,grd,w] = freqz_m(h,1); delta_w = pi/500;Asd = floor(-max(db([1:floor(ws1/delta_w)+1]))), % Actual AttnRpd = -min(db(ceil(wp1/delta_w)+1:floor(wp2/delta_w)+1)), % Actual passband ripple （5） %%% Filter Response Plotssubplot(2,2,1); stem(n,hd); title('Ideal Impulse Response: Bandpass'); axis([-1,M,min(hd)-0.1,max(hd)+0.1]); xlabel('n'); ylabel('hd(n)') set(gca,'XTickMode','manual','XTick',[0;M-1],'fontsize',10) subplot(2,2,2); stem(n,w_ham); title('Hamming Window'); axis([-1,M,-0.1,1.1]); xlabel('n'); ylabel('w_ham(n)')set(gca,'XTickMode','manual','XTick',[0;M-1],'fontsize',10) set(gca,'YTickMode','manual','YTick',[0;1],'fontsize',10)subplot(2,2,3); stem(n,h); title('Actual Impulse Response: Bandpass'); axis([-1,M,min(hd)-0.1,max(hd)+0.1]); xlabel('n'); ylabel('h(n)') set(gca,'XTickMode','manual','XTick',[0;M-1],'fontsize',10)subplot(2,2,4); plot(w/pi,db); title('Magnitude Response in dB');axis([0,1,-As-30,5]); xlabel('frequency in pi units'); ylabel('Decibels') set(gca,'XTickMode','manual','XTick',[0;0.3;0.4;0.5;0.6;1])set(gca,'XTickLabelMode','manual','XTickLabels',['0';'0.3';'0.4';'0.5';'0.6';'1'],... 'fontsize',10)set(gca,'TickMode','manual','YTick',[-50;0])set(gca,'YTickLabelMode','manual','YTickLabels',['-50';'0']);gridProblem 7Firstly, we use the given specifications ofs p s p A R ,,,ωω to design an analog lowpass IIR filter.Secondly, we change the analog lowpass IIR filter into the analog highpass IIR filter. Thirdly, we change the analog highpass IIR filter into the digital highpass IIR filter.。

1.题目及选项如下图所示A.AB.BC.CD.D【参考答案】: A2.下列关于FFT的说法中错误的是（）。

A.FFT是一种新的变换B.FFT是DFT的快速算法C.FFT基本上可以分成按时间抽取法和按频率抽取法两类D.基2 FFT要求序列的点数为（其中为整数）【参考答案】: A3.下列关于FIR滤波器的说法中正确的是（）。

A.FIR滤波器难以设计成线性相位特性B.FIR滤波器的单位冲激响应长度是无限的C.FIR滤波器的单位冲激响应长度是有限的D.对于相同的幅频特性要求，用FIR滤波器实现要比用IIR滤波器实现阶数低【参考答案】: C4.题目及选项如下图所示A.AB.BC.CD.D【参考答案】: A5.题目及选项如下图所示A.AB.BC.CD.D 【参考答案】: A6.题目及选项如下图所示A.AB.BC.CD.D 【参考答案】: A7.题目及选项如下图所示A.AB.BC.CD.D 【参考答案】: A8.题目及选项如下图所示A.AB.BC.CD.D【参考答案】: A9.在IIR数字滤波器的设计中，需将模拟参考滤波器从S平面到Z平面作单值映射，应采用（）。

A.脉冲响应不变法B.双线性变换法C.窗函数法D.频率采样法【参考答案】: B10.已知一连续时间信号的最高截止频率是4kHz，则要从该连续时间信号的均匀离散采样值无失真地恢复原信号，则采样频率为（）。

A.2kHzB.4kHzC.6kHzD.8kHz【参考答案】: D11.下列关于窗函数设计法的说法中错误的是（）。

A.窗函数的长度增加，则主瓣宽度减小，旁瓣宽度减小B.窗函数的旁瓣相对幅度取决于窗函数的形状，与窗函数的长度无关C.为减小旁瓣相对幅度而改变窗函数的形状，通常主瓣的宽度会增加D.对于长度固定的窗，只要选择合适的窗函数就可以使主瓣宽度足够窄，旁瓣幅度足够小【参考答案】: B12.已知某序列的Z变换收敛域为全Z平面，则该序列为（）。

A.有限长序列B.双边序列C.左边序列D.右边序列【参考答案】: A13.题目及选项如下图所示A.AB.BC.CD.D【参考答案】: A14.使用FFT算法计算两有限长序列的线性卷积，则需要调用（）次FFT算法。

数字信号处理实验报告班级：****姓名：郭**学号：*****联系方式：*****西安电子科技大学电子工程学院绪论数字信号处理起源于十八世纪的数学，随着信息科学和计算机技术的迅速发展，数字信号处理的理论与应用得到迅速发展，形成一门极其重要的学科。

当今数字信号处理的理论和方法已经得到长足的发展，成为数字化时代的重要支撑，其在各个学科和技术领域中的应用具有悠久的历史，已经渗透到我们生活和工作的各个方面。

数字信号处理相对于模拟信号处理具有许多优点，比如灵活性好，数字信号处理系统的性能取决于系统参数，这些参数很容易修改，并且数字系统可以分时复用，用一套数字系统可以分是处理多路信号；高精度和高稳定性，数字系统的运算字符有足够高的精度，同时数字系统不会随使用环境的变化而变化，尤其使用了超大规模集成的DSP 芯片，简化了设备，更提高了系统稳定性和可靠性；便于开发和升级，由于软件可以方便传送，复制和升级，系统的性能可以得到不断地改善；功能强，数字信号处理不仅能够完成一维信号的处理，还可以试下安多维信号的处理；便于大规模集成，数字部件具有高度的规范性，对电路参数要求不严格，容易大规模集成和生产。

数字信号处理用途广泛，对其进行一系列学习与研究也是非常必要的。

本次通过对几个典型的数字信号实例分析来进一步学习和验证数字信号理论基础。

实验一主要是产生常见的信号序列和对数字信号进行简单处理，如三点滑动平均算法、调幅广播（AM ）调制高频正弦信号和线性卷积。

实验二则是通过编程算法来了解DFT 的运算原理以及了解快速傅里叶变换FFT 的方法。

实验三是应用IRR 和FIR 滤波器对实际音频信号进行处理。

实验一●实验目的加深对序列基本知识的掌握理解●实验原理与方法1.几种常见的典型序列：0()1,00,0(){()()(),()sin()j n n n n u n x n Aex n a u n a x n A n σωωϕ+≥<====+单位阶跃序列：复指数序列：实指数序列：为实数 正弦序列：2.序列运算的应用：数字信号处理中经常需要将被加性噪声污染的信号中移除噪声，假定信号 s(n)被噪声d(n)所污染，得到了一个含噪声的信号()()()x n s n d n =+。

程序清单及波形显示： clc;close all;clear all;%======内容1：调用filter 解差分方程，由系统对u(n)的响应判断稳定性======A=[1,-0.9];B=[0.05,0.05]; %系统差分方程系数向量B 和A x1n=[1 1 1 1 1 1 1 1 zeros(1,50)]; %产生信号x1(n)=R8(n) x2n=ones(1,128); %产生信号x2(n)=u(n)y1n=filter(B,A,x1n); %求系统对x1(n)的响应y1(n) n=0:length(y1n)-1;subplot(2,2,1);stem(n,y1n,'.'); title('(a) 系统对R8(n)的响应y1(n)');xlabel('n');ylabel('y1(n)');y2n=filter(B,A,x2n); %求系统对x2(n)的响应y2(n) n=0:length(y2n)-1;subplot(2,2,2);stem(n,y2n,'.'); title('(b) 系统对u(n)的响应y2(n)');xlabel('n');ylabel('y2(n)');hn=impz(B,A,58); %求系统单位脉冲响应h(n) n=0:length(hn)-1;subplot(2,2,3);y=hn;stem(n,hn,'.'); title('(c) 系统单位脉冲响应h(n)');xlabel('n');ylabel('h(n)');(a) 系统对R8(n)的响应y1(n)ny 1(n )(b) 系统对u(n)的响应y2(n)ny 2(n )(c) 系统单位脉冲响应h(n)nh (n )%===内容2：调用conv 函数计算卷积============================x1n=[1 1 1 1 1 1 1 1 ]; %产生信号x1(n)=R8(n) h1n=ones(1,10); h2n=[1 2.5 2.5 1 ];y21n=conv(h1n,x1n); y22n=conv(h2n,x1n); figure(2)n=0:length(h1n)-1;subplot(2,2,1);stem(n,h1n); title('(d) 系统单位脉冲响应h1n');xlabel('n');ylabel('h1(n)'); n=0:length(y21n)-1;subplot(2,2,2);stem(n,y21n); title('(e) h1(n)与R8(n)的卷积y21n');xlabel('n');ylabel('y21(n)'); n=0:length(h2n)-1;subplot(2,2,3);stem(n,h2n); title('(f) 系统单位脉冲响应h2n');xlabel('n');ylabel('h2(n)'); n=0:length(y22n)-1;subplot(2,2,4);stem(n,y22n); title('(g) h2(n)与R8(n)的卷积y22n');xlabel('n');ylabel('y22(n)');nh 1(n )ny 21(n )(f) 系统单位脉冲响应h2nnh 2(n)(g) h2(n)与R8(n)的卷积y22nny 22(n )%=========内容3：谐振器分析======================== un=ones(1,256); %产生信号u(n) n=0:255;xsin=sin(0.014*n)+sin(0.4*n); %产生正弦信号A=[1,-1.8237,0.9801];B=[1/100.49,0,-1/100.49]; %系统差分方程系数向量B 和Ay31n=filter(B,A,un); %谐振器对u(n)的响应y31(n) y32n=filter(B,A,xsin); %谐振器对u(n)的响应y31(n) figure(3)n=0:length(y31n)-1;subplot(2,1,1);stem(n,y31n,'.'); title('(h) 谐振器对u(n)的响应y31n');xlabel('n');ylabel('y31(n)'); n=0:length(y32n)-1;subplot(2,1,2);stem(n,y32n,'.'); title('(i) 谐振器对正弦信号的响应y32n');xlabel('n');ylabel('y32(n)');050100150200250300(h) 谐振器对u(n)的响应y31nny 31(n )(i) 谐振器对正弦信号的响应y32nny 32(n )程序清单及波形显示：% DTMF 双频拨号信号产生6位电话号码 %clear all;clc;tm=[1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68]; % DTMF 信号代表的16个数N=205;K=[18,20,22,24,31,34,38,42];f1=[697,770,852,941]; % 行频率向量 f2=[1209,1336,1477,1633]; % 列频率向量 TN=input('键入6位电话号码= '); % 输入6位数字TNr=0; %接收端电话号码初值为零for l=1:6;d=fix(TN/10^(6-l)) TN=TN-d*10^(6-l); for p=1:4; for q=1:4;if tm(p,q)==abs(d); break,end % 检测码相符的列号q endif tm(p,q)==abs(d); break,end % 检测码相符的行号p endn=0:1023; % 为了发声，加长序列x = sin(2*pi*n*f1(p)/8000) + sin(2*pi*n*f2(q)/8000);% 构成双频信号sound(x,8000); % 发出声音 pause(0.1)% 接收检测端的程序X=goertzel(x(1:205),K+1); % 用Goertzel 算法计算八点DFT样本val = abs(X); % 列出八点DFT向量subplot(3,2,l);stem(K,val,'.');grid;xlabel('k');ylabel('|X(k)|') % 画出DFT(k)幅度axis([10 50 0 120])limit = 80; %for s=5:8;if val(s) > limit, break, end % 查找列号endfor r=1:4;if val(r) > limit, break, end % 查找行号endTNr=TNr+tm(r,s-4)*10^(6-l);enddisp('接收端检测到的号码为：') % 显示接收到的字符disp(TNr)显示结果：键入6位电话号码= 123456d = 1d = 2d = 3d = 4d = 5d = 6接收端检测到的号码为：123456050100|X (k )||X (k )|050100k|X (k )|050100k|X (k )|050100k|X (k )|050100k|X (k )|% DTMF 双频拨号信号产生8位电话号码 %clear all;clc;tm=[1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68]; % DTMF 信号代表的16个数 N=205;K=[18,20,22,24,31,34,38,42];f1=[697,770,852,941]; % 行频率向量 f2=[1209,1336,1477,1633]; % 列频率向量 TN=input('键入8位电话号码= '); % 输入8位数字TNr=0; %接收端电话号码初值为零for l=1:8;d=fix(TN/10^(8-l)) TN=TN-d*10^(8-l); for p=1:4; for q=1:4;if tm(p,q)==abs(d); break,end % 检测码相符的列号q endif tm(p,q)==abs(d); break,end % 检测码相符的行号p endn=0:1023; % 为了发声，加长序列 x = sin(2*pi*n*f1(p)/8000) + sin(2*pi*n*f2(q)/8000);% 构成双频信号 sound(x,8000); % 发出声音 pause(0.1)% 接收检测端的程序X=goertzel(x(1:205),K+1); % 用Goertzel 算法计算八点DFT 样本val = abs(X); % 列出八点DFT 向量 subplot(4,2,l);stem(K,val,'.');grid;xlabel('k');ylabel('|X(k)|') % 画出DFT(k)幅度axis([10 50 0 120])limit = 80; % for s=5:8;if val(s) > limit, break, end % 查找列号 endfor r=1:4;if val(r) > limit, break, end % 查找行号 endTNr=TNr+tm(r,s-4)*10^(8-l); enddisp('接收端检测到的号码为：') % 显示接收到的字符 disp(TNr) 显示结果：键入8位电话号码= 12345678 d = 1 d = 2 d = 3 d = 4 d = 5 d = 6 d = 7 d = 8接收端检测到的号码为：12345678|X (k )||X (k )||X (k )||X (k )||X (k )||X (k )|k|X (k )|k|X (k )|程序清单及波形显示： % 时域采样理论验证程序Tp=64/1000; %观察时间Tp=64微秒Fs=1000;T=1/Fs; M=Tp*Fs;n=0:M-1;t=n*T;A=444.128;alph=pi*50*2^0.5;omega=pi*50*2^0.5;xat=A*exp(-alph*t).*sin(omega*t);Xk=T*fft(xat,M); %M点FFT[xat)]subplot(3,2,1); stem(n,xat,'.'); xlabel('n');ylabel('x1(n)'); title('(a) Fs=1000Hz');k=0:M-1;fk=k/Tp;subplot(3,2,2);plot(fk,abs(Xk));title('(a) T*FT[xa(nT)],Fs=1000Hz');xlabel('\omega/hz');ylabel('(H1(ejw))');axis([0,Fs,0,1.2*max(abs(Xk))]);Fs=300;T=1/Fs; M=Tp*Fs;n=0:M-1;t=n*T;A=444.128;alph=pi*50*2^0.5;omega=pi*50*2^0.5;xat=A*exp(-alph*t).*sin(omega*t);Xk=T*fft(xat,M); %M点FFT[xat)]subplot(3,2,3); stem(n,xat,'.'); xlabel('n');ylabel('x2(n)'); title('(b)Fs=300Hz');k=0:M-1;fk=k/Tp;subplot(3,2,4);plot(fk,abs(Xk));title('(a) T*FT[xa(nT)],Fs=300Hz');xlabel('\omega/hz');ylabel('(H2(ejw))');axis([0,Fs,0,1.2*max(abs(Xk))]);Fs=200;T=1/Fs; M=Tp*Fs;n=0:M-1;t=n*T;A=444.128;alph=pi*50*2^0.5;omega=pi*50*2^0.5;xat=A*exp(-alph*t).*sin(omega*t);Xk=T*fft(xat,M); %M点FFT[xat)]subplot(3,2,5); stem(n,xat,'.'); xlabel('n');ylabel('x3(n)'); title('(c) Fs=200Hz');k=0:M-1;fk=k/Tp;subplot(3,2,6);plot(fk,abs(Xk));title('(a) T*FT[xa(nT)],Fs=200Hz');xlabel('\omega/hz');ylabel('(H3(ejw))');axis([0,Fs,0,1.2*max(abs(Xk))])nx 1(n )(a) Fs=1000Hz(a) T*FT[xa(nT)],Fs=1000Hzω/hz(H 1(e j w ))nx 2(n )(b) Fs=300Hz(a) T*FT[xa(nT)],Fs=300Hzω/hz(H 2(e j w ))nx 3(n )(c) Fs=200Hz(a) T*FT[xa(nT)],Fs=200Hzω/hz(H 3(e j w ))%频域采样理论验证程序 clc;clear;close all; M=27;N=32;n=0:M; xa=0:(M/2); xb= ceil(M/2)-1:-1:0; xn=[xa,xb]; %产生M 长三角波序列x(n) Xk=fft(xn,1024); %1024点FFT[x(n)], 用于近似序列x(n)的TF X32k=fft(xn,32) ;%32点FFT[x(n)] x32n=ifft(X32k); %32点IFFT[X32(k)]得到x32(n) X16k=X32k(1:2:N); %隔点抽取X32k 得到X16(K) x16n=ifft(X16k,N/2); %16点IFFT[X16(k)]得到x16(n) subplot(3,2,2);stem(n,xn,'.'); title('(b) 三角波序列x(n)');xlabel('n');ylabel('x(n)');axis([0,32,0,20]) k=0:1023;wk=2*k/1024; subplot(3,2,1);plot(wk,abs(Xk));title('(a)FT[x(n)]'); xlabel('\omega/\pi');ylabel('|X(e^j^\omega)|');axis([0,1,0,200]) k=0:N/2-1; subplot(3,2,3);stem(k,abs(X16k),'.'); title('(c) 16点频域采样');xlabel('k');ylabel('|X_1_6(k)|');axis([0,8,0,200]) n1=0:N/2-1; subplot(3,2,4);stem(n1,x16n,'.'); title('(d) 16点IDFT[X_1_6(k)]');xlabel('n');ylabel('x_1_6(n)');axis([0,32,0,20]); k=0:N-1; subplot(3,2,5);stem(k,abs(X32k),'.'); title('(e) 32点频域采样');xlabel('k');ylabel('|X_3_2(k)|');axis([0,16,0,200]);n1=0:N-1; subplot(3,2,6);stem(n1,x32n,'.');box on title('(f) 32点IDFT[X_3_2(k)]');xlabel('n');ylabel('x_3_2(n)');axis([0,32,0,20])(b) 三角波序列x(n)nx (n )0100200(a)FT[x(n)]ω/π|X (e j ω)|(c) 16点频域采样k|X 16(k)|102030(d) 16点IDFT[X 16(k)]nx 16(n)(e) 32点频域采样k|X 32(k )|(f) 32点IDFT[X 32(k)]nx 32(n )程序清单及波形显示：% 用FFT 对信号作频谱分析 clear all;close all %实验内容(1)=================================================== x1n=[ones(1,4)]; %产生序列向量x1(n)=R4(n) M=8;xa=1:(M/2); xb=(M/2):-1:1; x2n=[xa,xb]; %产生长度为8的三角波序列x2(n) x3n=[xb,xa]; X1k8=fft(x1n,8); %计算x1n 的8点DFT X1k16=fft(x1n,16); %计算x1n 的16点DFT X2k8=fft(x2n,8); %计算x1n 的8点DFT X2k16=fft(x2n,16); %计算x1n 的16点DFT X3k8=fft(x3n,8); %计算x1n 的8点DFT X3k16=fft(x3n,16); %计算x1n 的16点DFT %以下绘制幅频特性曲线 subplot(1,2,1);stem(X1k8,'.'); %绘制8点DFT 的幅频特性图 title('(1a) 8点DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度'); subplot(1,2,2);stem(X1k16,'.'); %绘制16点DFT 的幅频特性图 title('(1b)16点DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度');figure(2)subplot(2,2,1);stem(X2k8,'.'); %绘制8点DFT 的幅频特性图 title('(2a) 8点DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度');subplot(2,2,2);stem(X2k16,'.'); %绘制16点DFT 的幅频特性图 title('(2b)16点DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度'); subplot(2,2,3);stem(X3k8,'.'); %绘制8点DFT 的幅频特性图 title('(3a) 8点DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度');subplot(2,2,4);stem(X3k16,'.'); %绘制16点DFT 的幅频特性图 title('(3b)16点DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度');-2.5-2-1.5-1-0.500.511.522.5(1a) 8点DFT[x 1(n)]ω/π幅度-2.5-2-1.5-1-0.500.511.522.5(1b)16点DFT[x 1(n)]ω/π幅度2468-4-2024(2a) 8点DFT[x 2(n)]ω/π幅度5101520-20-1001020(2b)16点DFT[x 2(n)]ω/π幅度2468-4-2024(3a) 8点DFT[x 3(n)]ω/π幅度-10-50510(3b)16点DFT[x 3(n)]ω/π幅度%实验内容(2) 周期序列谱分析==================================N=8;n=0:N-1; %FFT 的变换区间N=8 x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k8=fft(x4n); %计算x4n 的8点DFT X5k8=fft(x5n); %计算x5n 的8点DFT N=16;n=0:N-1; %FFT 的变换区间N=16 x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k16=fft(x4n); %计算x4n 的16点DFT X5k16=fft(x5n); %计算x5n 的16点DFT figure(3)subplot(2,2,1);stem(X4k8,'.'); %绘制8点DFT 的幅频特性图 title('(4a) 8点DFT[x_4(n)]');xlabel('ω/π');ylabel('幅度');subplot(2,2,3);stem(X4k16,'.'); %绘制16点DFT 的幅频特性图 title('(4b)16点DFT[x_4(n)]');xlabel('ω/π');ylabel('幅度');subplot(2,2,2);stem(X5k8,'.'); %绘制8点DFT 的幅频特性图 title('(5a) 8点DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度');subplot(2,2,4);stem(X5k16,'.'); %绘制16点DFT 的幅频特性图title('(5b)16点DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度');2468-1012x 10-15(4a) 8点DFT[x 4(n)]ω/π幅度5101520-4-2024x 10-15(4b)16点DFT[x 4(n)]ω/π幅度2468-4-2024(5a) 8点DFT[x 5(n)]ω/π幅度5101520-4-2024x 10-15(5b)16点DFT[x 5(n)]ω/π幅度%实验内容(3) 模拟周期信号谱分析=============================== figure(4) Fs=64;T=1/Fs; N=16;n=0:N-1; %FFT 的变换区间N=16 x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)16点采样 X6k16=fft(x6nT); %计算x6nT 的16点DFT X6k16=fftshift(X6k16); %将零频率移到频谱中心 Tp=N*T;F=1/Tp; %频率分辨率F k=-N/2:N/2-1;fk=k*F; %产生16点DFT 对应的采样点频率（以零频率为中心） subplot(3,1,1);stem(fk,abs(X6k16),'.');box on %绘制8点DFT 的幅频特性图 title('(6a) 16点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度'); axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16))])N=32;n=0:N-1; %FFT的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)32点采样X6k32=fft(x6nT); %计算x6nT的32点DFTX6k32=fftshift(X6k32); %将零频率移到频谱中心Tp=N*T;F=1/Tp; %频率分辨率Fk=-N/2:N/2-1;fk=k*F; %产生16点DFT对应的采样点频率（以零频率为中心）subplot(3,1,2);stem(fk,abs(X6k32),'.');box on %绘制8点DFT的幅频特性图title('(6b) 32点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度');axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32))])N=64;n=0:N-1; %FFT的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)64点采样X6k64=fft(x6nT); %计算x6nT的64点DFTX6k64=fftshift(X6k64); %将零频率移到频谱中心Tp=N*T;F=1/Tp; %频率分辨率Fk=-N/2:N/2-1;fk=k*F; %产生16点DFT对应的采样点频率（以零频率为中心）subplot(3,1,3);stem(fk,abs(X6k64),'.'); box on%绘制8点DFT的幅频特性图title('(6a) 64点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度');axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k64))])-30-20-1001020300510(6a) 16点|DFT[x 6(nT)]|f(Hz)幅度010(6b) 32点|DFT[x 6(nT)]|f(Hz)幅度020(6a) 64点|DFT[x 6(nT)]|f(Hz)幅度程序清单及波形显示： clc;clear all;close all;fc1=250; fm1=15; fc2=500; fm2=50; fc3=1000; fm3=100;N=800; Fs=10000;Ts=1/Fs; n=[0:N-1];t=n*Ts;x11=cos(2*pi*fc1*t); x12=cos(2*pi*fm1*t); x1=x11.*x12; subplot(3,1,1);plot(t,x11,'g');plot(t,x12,'r');plot(t,x1,'b'); legend('载波','调制波 ','已调 ');xlabel('t/s');ylabel('波形')x=cos(2*pi*fc1*t).*cos(2*pi*fm1*t)+cos(2*pi*fc2*t).*cos(2*pi*fm2*t)+cos(2*pi*fc3*t).*cos(2*pi*fm3*t); subplot(3,1,2);plot(t,x);X=fft(x)subplot(3,1,3)k=[0:(N-1)/2]stem(k*2/N,abs(X(k+1))/max(abs(X(k+1))),'.');axis([0,0.3,0,1]);xlabel('\omeg a/\pi');ylabel('幅度');wp=[0.04,0.06];ws=[0.03,0.07];rp=0.1;rs=60;[N1,wp]=ellipord(wp,ws,rp,rs);[B,A]=ellip(N1,rp,rs,wp);y1=filter(B,A,x);figuresubplot(3,1,1);plot(t,x);xlabel('t/s');title('3路混合信号波形')subplot(3,1,2);[H1,w1]=freqz(B,A,N);plot(w1/pi,20*log10(abs(H1)));axis([0,0.5,-80,1]);xlabel('\omega/\pi');ylabel('|H(e^j\omega)|');title('中心频率为250Hz的频率响应');subplot(3,1,3);plot(t,y1);xlabel('t/s');ylabel('y1(t)');title('中心频率为250H的滤波信号')wp=[0.08,0.12];ws=[0.07,0.13];rp=0.1;rs=60;[N1,wp]=ellipord(wp,ws,rp,rs);[B,A]=ellip(N1,rp,rs,wp);y1=filter(B,A,x);figuresubplot(3,1,1);plot(t,x);xlabel('t/s');title('3路混合信号波形');subplot(3,1,2)[H1,w1]=freqz(B,A,N);plot(w1/pi,20*log10(abs(H1)));axis([0,0.5,-90,2]);xlabel('\omega/\pi');ylabel('|H(e^j\omega)|');title('中心频率为500Hz的频率响应')subplot(3,1,3);plot(t,y1);xlabel('t/s');ylabel('y1(t)');title('中心频率为500H的滤波信号')wp=[0.17,0.23];ws=[0.16,0.24];rp=0.1;rs=60;[N1,wp]=ellipord(wp,ws,rp,rs);[B,A]=ellip(N1,rp,rs,wp);y1=filter(B,A,x);figuresubplot(3,1,1);plot(t,x);xlabel('t/s');title('3路混合信号波形');subplot(3,1,2);[H1,w1]=freqz(B,A,N);plot(w1/pi,20*log10(abs(H1)));axis([0,0.5,-100,10]);xlabel('\omega/\pi');ylabel('|H(e^j\omega)|');title('中心频率为1000Hz的频率响应')subplot(3,1,3);plot(t,y1);xlabel('t/s');ylabel('y1(t)');title('中心频率为1000H的滤波信号')00.010.020.030.040.050.060.070.08-101t/s波形00.010.020.030.040.050.060.070.08-505ω/π幅度-505t/s3路混合信号波形-80-60-40-200ω/π|H (e j ω)|中心频率为250Hz 的频率响应-202t/sy 1(t )中心频率为250H 的滤波信号0.010.020.030.040.050.060.070.08-505t/s3路混合信号波形0.050.10.150.20.250.30.350.40.450.5-80-60-40-200ω/π|H (e j ω)|中心频率为500Hz 的频率响应0.010.020.030.040.050.060.070.08-202t/sy 1(t )中心频率为500H 的滤波信号0.010.020.030.040.050.060.070.08-505t/s3路混合信号波形0.050.10.150.20.250.30.350.40.450.5-100-500ω/π|H (e j ω)|中心频率为1000Hz 的频率响应0.010.020.030.040.050.060.070.08-202t/sy 1(t )中心频率为1000H 的滤波信号程序清单及波形显示： clc;clear;clear allN=1000; Fs=1000;T=1/Fs;Tp=N*T; t=0:T:(N-1)*T;fc=Fs/10;f0=fc/10; %载波频率fc=Fs/10，单频调制信号频率为f0=Fc/10; mt=cos(2*pi*f0*t); %产生单频正弦波调制信号mt ，频率为f0 ct=cos(2*pi*fc*t); %产生载波正弦波信号ct ，频率为fc xt=mt.*ct; %相乘产生单频调制信号xt nt=2*rand(1,N)-1; %产生随机噪声ntfp=150; fs=200;Rp=0.1;As=70; % 滤波器指标fb=[fp,fs];m=[0,1]; % 计算remezord 函数所需参数f,m,devdev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];[n,fo,mo,W]=remezord(fb,m,dev,Fs); % 确定remez 函数所需参数hn=remez(n,fo,mo,W); % 调用remez 函数进行设计,用于滤除噪声nt 中的低频成分yt=filter(hn,1,10*nt); %滤除随机噪声中低频成分，生成高通噪声yt xt=xt+yt; %噪声加信号 fst=fft(xt,N);k=0:N-1;f=k/Tp;subplot(2,1,1);plot(t,xt);grid;xlabel('t/s');ylabel('x(t)');axis([0,Tp/5,min(xt),max(xt)]);title('(a) 信号加噪声波形')subplot(2,1,2);plot(f,abs(fst)/max(abs(fst)));grid;title('(b) 信号加噪声的频谱')axis([0,Fs/2,0,1.2]);xlabel('f/Hz');ylabel('幅度')-10-50510t/sx (t )(a) 信号加噪声波形0501001502002503003504004505000.51(b) 信号加噪声的频谱f/Hz幅度%==调用xtg 产生信号xt, xt 长度N=1000,并显示xt 及其频谱,========= N=1000;xt=xtg(N);fp=120; fs=150;Rp=0.2;As=60;Fs=1000; % 输入给定指标% (1) 用窗函数法设计滤波器wc=(fp+fs)/Fs; %理想低通滤波器截止频率(关于pi归一化）B=2*pi*(fs-fp)/Fs; %过渡带宽度指标Nb=ceil(11*pi/B); %blackman窗的长度Nhn=fir1(Nb-1,wc,blackman(Nb));Hw=abs(fft(hn,1024)); % 求设计的滤波器频率特性ywt=fftfilt(hn,xt,N); %调用函数fftfilt对xt滤波subplot(4,1,1);plot(f,Hw);xlabel('f/Hz');ylabel('幅度'); title('(a)低通滤波器幅频特性');subplot(4,1,2);plot(t,ywt); title('(b)滤除噪声后的信号波形');xlabel('t/s');ylabel('ywt');% (2) 用等波纹最佳逼近法设计滤波器fb=[fp,fs];m=[1,0]; % 确定remezord函数所需参数f,m,dev dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)];[Ne,fo,mo,W]=remezord(fb,m,dev,Fs); % 确定remez函数所需参数hn=remez(Ne,fo,mo,W); % 调用remez函数进行设计Hw=abs(fft(hn,1024)); % 求设计的滤波器频率特性yet=fftfilt(hn,xt,N); % 调用函数fftfilt对xt滤波subplot(3,1,1);plot(f,Hw);xlabel('f/Hz');ylabel('幅度'); title('(c)低通滤波器幅频特性');subplot(3,1,2);plot(t,yet);title('(d)滤除噪声后的信号波形');xlabel('t/s');ylabel('yet');。

浅谈奈奎斯特频率采样和压缩感知信息技术的飞速发展使得人们对信息的需求量剧增。

现实世界的模拟化和信号处理工具的数字化决定了信号采样是从模拟信源获取数字信息的必经之路。

在信号和图像处理领域，凡是涉及到计算机作为处理工具的场合，所面临的首要问题就是模拟信号的数字化问题，然后再对得到的离散的样本进行各种处理。

连续信号转化为离散的数字化信号的过程称为采样。

对模拟信号采样所得的离散数字信号能否代表并恢复成原来的连续模拟信号呢？如能恢复应具备什么样的条件呢？这个问题直接关系到是否可以用数字处理工具和数字化的方法处理模拟信号。

一奈奎斯特频率采样奈奎斯特采样定理给我们提供了如何采样的重要理论基础。

它指出，如果信号是带限的，采样速率必须达到信号带宽的两倍以上才能精确重构信号。

事实上，在音频和可视电子设备、医学图像设备、无线接收设备等设备中的所有信号采样协议都隐含了这样的限制。

奈奎斯特采样定理至出现以来一直是数字信号和图像处理领域的重要理论基础，它支撑着几乎所有的信号和图像处理过程，包括信号和图像的获取、存储、处理、传输等。

采样定理，又称香农采样定理，奈奎斯特采样定理，是信息论，特别是通讯与信号处理学科中的一个重要基本结论.E.T.Whittaker (1915年发表的统计理论),克劳德·香农与Harry Nyquist都对它作出了重要贡献。

另外,V. A. Kotelnikov 也对这个定理做了重要贡献。

采样是将一个信号（即时间或空间上的连续函数）转换成一个数值序列（即时间或空间上的离散函数）。

采样定理指出，如果信号是带限的，并且采样频率高于信号带宽的一倍，那么，原来的连续信号可以从采样样本中完全重建出来。

带限信号变换的快慢受到它的最高频率分量的限制，也就是说它的离散时刻采样表现信号细节的能力是有限的。

采样定理是指，如果信号带宽小于奈奎斯特频率（即采样频率的二分之一），那么此时这些离散的采样点能够完全表示原信号。

数字信号处理上机实验报告14020710021 张吉凯第一次上机实验一：设给定模拟信号()1000t a x t e -=，t 的单位是ms 。

(1) 利用MATLAB 绘制出其时域波形和频谱图（傅里叶变换），估计其等效带宽（忽略谱分量降低到峰值的3%以下的频谱）。

(2) 用两个不同的采样频率对给定的()a x t 进行采样。

○1()()15000s a f x t x n =以样本秒采样得到。

()()11j x n X e ω画出及其频谱。

○2()()11000s a f x t x n =以样本秒采样得到。

()()11j x n X e ω画出及其频谱。

比较两种采样率下的信号频谱，并解释。

（1）MATLAB 程序：N=10; Fs=5; T s=1/Fs;n=[-N:T s:N];xn=exp(-abs(n)); w=-4*pi:0.01:4*pi; X=xn*exp(-j*(n'*w)); subplot(211) plot(n,xn);title('x_a(t)时域波形');xlabel('t/ms');ylabel('x_a(t)'); axis([-10, 10, 0, 1]); subplot(212);plot(w/pi,abs(X)); title('x_a(t)频谱图');xlabel('\omega/\pi');ylabel('X_a(e^(j\omega))');ind = find(X >=0.03*max(X))*0.01; eband = (max(ind) -min(ind)); fprintf('等效带宽为%fKHZ\n',eband); 运行结果：等效带宽为12.110000KHZ（2）MATLAB程序：N=10;omega=-3*pi:0.01:3*pi;%Fs=5000Fs=5;T s=1/Fs;n=-N:T s:N;xn=exp(-abs(n));X=xn*exp(-j*(n'*omega));subplot(2,2,1);stem(n,xn);grid on;axis([-10, 10, 0, 1.25]); title('时域波形(f_s=5000)');xlabel('n');ylabel('x_1(n)');subplot(2,2,2);plot(omega/pi,abs(X));title('频谱图(f_s=5000)');xlabel('\omega/\pi');ylabel('X_1(f)');grid on;%Fs=1000Fs=1;T s=1/Fs;n=-N:T s:N;xn=exp(-abs(n));X=xn*exp(-j*(n'*omega));subplot(2,2,3);stem(n,xn);grid on;axis([-10, 10, 0, 1.25]); title('时域波形(f_s=1000)');xlabel('n');ylabel('x_2(n)');grid on;subplot(2,2,4);plot(omega/pi,abs(X));title('频谱图(f_s=1000)');xlabel('\omega/\pi');ylabel('X_2(f)');grid on;运行结果：实验二：给定一指数型衰减信号()()0cos 2at x t e f t π-=，采样率1s f T=，T 为采样周期。

西安电子科技大学数字信号处理（DSP）课程实验报告实验名称 VISUAL DSP++的使用入门电子工程学院 1402071 班同作者崔健孟智超杨伟祺实验日期 2017 年 5 月 16 日实验一：VISUAL DSP++的使用入门一、实验目的：熟悉VISUAL DSP++的开发环境二、实验内容:练习一：启动Visual DSP++，建立一个用C源代码的工程（Project），同时用调试器来评估用C语言所编写代码的性能；练习二：创立一个新的工程，修改源码来调用一个汇编（asm）程序，重新编译工程，用调试器来评估用汇编语言所写程序的性能；练习三：利用调试器的绘图（plot）功能来图形显示一个卷积算法中的多个数据的波形；练习四：利用调试器的性能统计功能（Statistical profile）来检查练习三中卷积算法的效率。

利用所收集到的性能统计数据就能看出算法中最耗时的地方。

三、实验步骤及实验结果：练习一：1)新建工程进入 Visual DSP＋＋，显示Visual DSP++的集成开发和调试环境窗口，选择菜单File 中Open 打开文件：…\unit_1\dot_product_c \dotprodc．dpj。

2)编译 dotprodc工程在菜单 Project中选择 Build Project来对工程进行编译。

在本例子中，编译器会检测到一个未定义的错误，显示为：“．＼dotprod_main.c”，line 115：error ＃20：identifier“itn”is undefined itn i；将该错误改正后，保存并重新编译，没有错误出现，编译成功。

3)运行VsualDSP++调试器在编译完成后，环境将自动进入调试状态，对于初次进入debugger，将显示对象选择对话框，在其中指定对象和处理器信息。

4)运行dotprod.c从 Debug菜单中选择 Run项，程序将被执行，其输出结果在 Output window中显示。

% 数字信号处理大作业语音信号提取程序加噪声程序
% 功能:从文件名为''语音文件中提取出数据，并加入声音最大音量五分之一的噪声% 参数：name-原始信号；Y-加噪声后的信号；
% 大作业内容：请自己录音，编写FIR滤波程序，消减噪声的影响。

% 作者：西安电子科技大学雷达信号处理国家重点实验室李军
% 时间：2015年
close all;clear all;clc;
namex=1e5;
% 读取原始声音文件
[name,fs0,bit]=wavread('',[100,namex]); %读取源文件
sound(name,fs0,bit);
namex=length(name);
noise=(max(name(:,1))/5)*randn(namex,2);
Y=name+noise; % 加入噪声
wavwrite(Y,fs0,bit,'mingyun_'); % 写入噪声文件
sound(Y,fs0,bit);
% 作图分析声音信号的时域特性和频域特性
figure(1)
subplot(2,1,1);
plot(abs(name));
title('（a）原声音信号的时域特性');
subplot(2,1,2);
plot(abs(Y));
title('(b)含有噪声声音信号的时域特性');。

数字信号处理上机大作业实验一：信号、系统及系统响应(1) 简述实验目的及实验原理。

1.实验目的●熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。

●熟悉时域离散系统的时域特性。

●利用卷积方法观察分析系统的时域特性。

●掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。

2.实验原理与方法●时域采样。

● LTI系统的输入输出关系。

(2)按实验步骤附上实验过程中的信号序列、系统单位脉冲响应及系统响应序列的时域和幅频特性曲线，并对所得结果进行分析和解释。

Matlab源程序如下：A=1;T1=1/1000;T2=1/300;T3=1/200;a=25*pi;w0=30*pi;n=0:99;x1=A*exp(-a*n*T1).*sin(w0*n*T1);x2=A*exp(-a*n*T2).*sin(w0*n*T2);x3=A*exp(-a*n*T3).*sin(w0*n*T3);m=linspace(-pi,pi,10000);X1=x1*exp(-j*n'*m);%n'与m构造矩阵，xi向量与矩阵每一列相乘对应元素相加，构成DTFT后的矩阵X2=x2*exp(-j*n'*m);X3=x3*exp(-j*n'*m);figure(1);subplot(3,2,1)plot(m/pi,abs(X1));xlabel('\omega/π');ylabel('|H(e^j^\omega)|');title('采样频率为1000Hz时的幅度谱');subplot(3,2,3)plot(m/pi,abs(X2));xlabel('\omega/π');ylabel('|H(e^j^\omega)|');title('采样频率为300Hz时的幅度谱');subplot(3,2,5)plot(m/pi,abs(X3));xlabel('\omega/π');ylabel('|H(e^j^\omega)|');title('采样频率为200Hz时的幅度谱');subplot(3,2,2)plot(n,abs(x1));xlabel('n');ylabel('x1(t)');title('采样频率为1000Hz时的时域波形');subplot(3,2,4)plot(n,abs(x2));xlabel('n');ylabel('x2(t)');title('采样频率为300Hz时的时域波形');subplot(3,2,6)plot(n,abs(x3));xlabel('n');ylabel('x3(t)');title('采样频率为200Hz时的时域波形');波形图如下：-1-0.8-0.6-0.4-0.200.20.40.60.81ω/π|H (e j ω)|采样频率为1000Hz 时的幅度谱ω/π|H (e j ω)|采样频率为300Hz 时的幅度谱ω/π|H (e j ω)|采样频率为200Hz 时的幅度谱102030405060708090100nx 1(t )采样频率为1000Hz 时的时域波形nx 2(t )采样频率为300Hz 时的时域波形nx 3(t )采样频率为200Hz 时的时域波形② 时域离散信号、 系统和系统响应分析。

实验一、信号的采样clc,clear;dt=0.001;tf=6;t=0:dt:tf;xa=sqrt(t)+cos(t);T=0.5;n=0:tf/T;x=sqrt(n*T)+cos(n*T);figure(1)subplot(2,1,1)plot(t,xa),grid on ;title('original image')subplot(2,1,2)stem(n*T,x),grid on ,title('digital image')实验二、信号与系统的时域分析差分方程为)()2()1()(21n bx n y a n y a n y +----=，其中8.01-=a ，64.02=a ，866.0=b 。

系统单位脉冲响应)(n ha1=-0.8;a2=0.64;b=0.866;ys=0;xn=[1,zeros(1,49)];B=1;A=[1,a1,a2];xi=filtic(B,A,ys);yn=filter(B,A,xn,xi);n=0:length(yn)-1;subplot(1,1,1);stem(n,yn,'.')title('(a)');xlabel('n');ylabel('y(n)')输入x(n)=cos(n)T=0.1;z=cos(n*T);zn=conv(yn,z); figure(2);n1=1:99;stem(n1,zn,'.')实验三、系统的频域和Z域分析程序代码（画出dtft的幅度和频率谱）clc,clear;n=0:1:7;x=(0.9*exp(j*pi/3)).^n;w=0:pi/200:pi;X=x*exp(-j).^(n'*w);realX=real(X);imagX=imag(X);angX=angle(X);magX=abs(X);subplot(2,2,1);plot(w/pi,magX);grid xlabel('frequency in pi unit');title('magnitude part');subplot(2,2,2);plot(w/pi,realX);grid xlabel('frequency in pi unit');title('real part');subplot(2,2,3);plot(w/pi,imagX);grid xlabel('frequency in pi unit');title('imaginary part');subplot(2,2,4);plot(w/pi,angX);grid xlabel('frequency in pi unit');title('angel part');clc,clear;a=[1,-0.5,0.06];b=[1,1,0];m=0:length(b)-1;l=0:length(a)-1;w=0:pi/500:pi;num=b*exp(-j*m'*w);den=a*exp(-j*l'*w);H=num./den;magH=abs(H);angH=angle(H);H1=freqz(b,a,w);magH1=abs(H1);angH1=angle(H1);subplot(2,2,2);plot(w/pi,angH);grid;xlabel('w(frequency in pi units)');ylabel('Ïàλrad/w');subplot(2,2,1);plot(w/pi,magH);grid;xlabel('w(frequency in pi units)');ylabel('·ù¶È|H|');subplot(2,2,3);plot(w/pi,magH1);grid;xlabel('w(frequency in pi units)');ylabel('·ù¶È|H1|');subplot(2,2,4);plot(w/pi,angH);grid;xlabel('w(frequency in pi units)');ylabel('Ïàλrad/w');axis([0,1,-0.8,0]); figure(2);zplane(b,a);实验四、信号的频谱分析程序代码clc,clear;n=0:7;k=0:7;N=8;w=n*(2*pi)/8;x=(0.9*exp(j*pi/3)).^n;X1=[x zeros(1,8)];X2=[X1 zeros(1,16)];XK=x*exp(-j*k'*w);k1=0:15;n1=0:15;w1=n1*(2*pi)/16;XK1=X1*exp(-j*k1'*w1);k2=0:31;n2=0:31;w2=n2*(2*pi)/16;XK2=X2*exp(-j*k2'*w2);w3=0:pi/200:2*pi;X=x*exp(-j*n'*w3);magX=abs(X);angX=angle(X);magXK=abs(XK);angXK=angle(XK);magXK1=abs(XK1);angXK1=angle(XK1);magXK2=abs(XK2);angXK2=angle(XK2);subplot(4,2,1);plot(w3/pi,magX);xlabel('w/pi');ylabel('·ù¶È|X|');grid on;subplot(4,2,2);plot(w3/pi,angX);xlabel('w/pi');ylabel('Ïàλrad/pi'); subplot(4,2,3);stem(n,magXK);xlabel('K');ylabel('·ù¶È|XK|');subplot(4,2,4);stem(n,magXK);xlabel('K');ylabel('Ïàλrad/pi'); subplot(4,2,5);stem(n1,magXK1);xlabel('K1');ylabel('·ù¶È|XK1|'); subplot(4,2,6);stem(n1,magXK1);xlabel('K1');ylabel('Ïàλrad/pi'); subplot(4,2,7);stem(n2,magXK2);xlabel('K2');ylabel('·ù¶È|XK2|'); subplot(4,2,8);stem(n2,magXK2);xlabel('K2');ylabel('Ïàλrad/pi');实验五、IIR数字滤波器设计IIR汉宁窗低通高通低通巴特沃斯通带截止频率wp=0.2pi 通带最大衰减R=1dB阻带截止频率wp=0.35pi 阻带最大衰减R=10dBclc,clear;Wp=0.2;Ws=0.35;Rp=1;Rs=100;[N,Wc]=buttord(Wp,Ws,Rp,Rs);[Bz,Az]=butter(N,Wc)w=0:0.1:pi;[H,w1]=freqz(Bz,Az,w);;ang=angle(H);H=20*log10(abs(H))subplot(4,2,1); plot(w/pi,H) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('µÍͨÂ˲¨Æ÷')subplot(4,2,2);plot(w/pi,ang);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')[Bz1,Az1]=butter(N,Wc,'high')w=0:0.1:pi;[H1,w2]=freqz(Bz1,Az1,w);ang1=angle(H1);H1=20*log10(abs(H1))subplot(4,2,3); plot(w/pi,H1) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('¸ßͨÂ˲¨Æ÷')subplot(4,2,4);plot(w/pi,ang1);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')Wp1=[0.2 0.8];Ws1=[0.35 0.65];[N2,Wc1]=buttord(Wp1,Ws1,Rp,Rs);[Bz2,Az2]=butter(N2,Wc1,'stop')w=0:0.1:pi;[H2,w3]=freqz(Bz2,Az2,w);ang2=angle(H2);H2=20*log10(abs(H2))subplot(4,2,5); plot(w/pi,H2) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('´ø×èÂ˲¨Æ÷')subplot(4,2,6);plot(w/pi,ang2);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')Wp1=[0.2 0.8];Ws1=[0.35 0.65];[N2,Wc1]=buttord(Wp1,Ws1,Rp,Rs);[Bz3,Az3]=butter(N2,Wc1)w=0:0.1:pi;[H3,w4]=freqz(Bz3,Az3,w);ang3=angle(H3);H3=20*log10(abs(H3))subplot(4,2,7); plot(w/pi,H3) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('´øͨÂ˲¨Æ÷')subplot(4,2,8);plot(w/pi,ang3);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')切比雪夫1型通带截止频率wp=0.7pi 通带最大衰减R=1dB阻带截止频率wp=0.5pi 阻带最大衰减R=40dBclc,clear;Wp=0.7;Ws=0.5;Rp=1;Rs=40;[N,Wpo]=cheb1ord(Wp,Ws,Rp,Rs);[Bz,Az]=cheby1(N,Rp,Wpo)w=0:0.1:pi;[H,w1]=freqz(Bz,Az,w);ang=angle(H);H=20*log10(abs(H))subplot(4,2,1); plot(w/pi,H) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('µÍͨÂ˲¨Æ÷')subplot(4,2,2);plot(w/pi,ang);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')[Bz1,Az1]=cheby1(N,Rp,Wpo,'high');w=0:0.1:pi;[H1,w2]=freqz(Bz1,Az1,w);ang1=angle(H1);H1=20*log10(abs(H1))subplot(4,2,3); plot(w/pi,H1) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('¸ßͨÂ˲¨Æ÷')subplot(4,2,4);plot(w/pi,ang1);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')Wp1=[0.2 0.8];Ws1=[0.35 0.65];[N2,Wpo1]=cheb1ord(Wp1,Ws1,Rp,Rs);[Bz2,Az2]=cheby1(N2,Rp,Wpo1,'stop')w=0:0.1:pi;[H2,w3]=freqz(Bz2,Az2,w);ang2=angle(H2);H2=20*log10(abs(H2))subplot(4,2,5); plot(w/pi,H2) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('´ø×èÂ˲¨Æ÷')subplot(4,2,6);plot(w/pi,ang2);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')Wp1=[0.2 0.8];Ws1=[0.35 0.65];[N2,Wpo1]=cheb1ord(Wp1,Ws1,Rp,Rs);[Bz3,Az3]=cheby1(N2,Rp,Wpo1)w=0:0.1:pi;[H3,w4]=freqz(Bz3,Az3,w);ang3=angle(H3);H3=20*log10(abs(H3))subplot(4,2,7); plot(w/pi,H3) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('´øͨÂ˲¨Æ÷')subplot(4,2,8);plot(w/pi,ang3);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')实验六、FIR数字滤波器设计FIR汉宁窗低通高通低通% 采用Hamming窗设计一个带阻FIR滤波器阻带：0~0.5pi，阻带最小衰减Rs=40dB；通带：0.5~pi，通带最大衰减：Rp=1dB。