西安电子科技大学数字信号处理大作业
西安电子科技大学数字信号处理上机作业
数字信号处理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大作业任务报告
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的⾼斯噪声。
数字信号处理实验报告(西电)
数字信号处理实验报告班级：****姓名：郭**学号：*****联系方式：*****西安电子科技大学电子工程学院绪论数字信号处理起源于十八世纪的数学，随着信息科学和计算机技术的迅速发展，数字信号处理的理论与应用得到迅速发展，形成一门极其重要的学科。
当今数字信号处理的理论和方法已经得到长足的发展，成为数字化时代的重要支撑，其在各个学科和技术领域中的应用具有悠久的历史，已经渗透到我们生活和工作的各个方面。
数字信号处理相对于模拟信号处理具有许多优点，比如灵活性好，数字信号处理系统的性能取决于系统参数，这些参数很容易修改，并且数字系统可以分时复用，用一套数字系统可以分是处理多路信号；高精度和高稳定性，数字系统的运算字符有足够高的精度，同时数字系统不会随使用环境的变化而变化，尤其使用了超大规模集成的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');。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数字信号处理大作业班级：021231 学号：姓名：指导老师：吕雁写出奈奎斯特米样率和和信号稀疏米样的学习报告和体会1、采样定理___ T在进行A/D信号的转换过程中，当采样频率fs.max大于信号中最高频率fmax的2倍时(fs.max>2fmax)，采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的5〜10倍；采样定理又称奈奎斯特定理。
(1) 在时域频带为F的连续信号f(t)可用一系列离散的采样值f(t1),f(t1 ±A t), f(t1 ±2A t),...来表示,只要这些采样点的时间间隔A t < 1/2F，便可根据各采样值完全恢复原始信号。
(2) 在频域当时间信号函数f(t)的最高频率分量为fmax时,f(t)的值可由一系列采样间隔小于或等于1/2fo的采样值来确定，即采样点的重复频率fs >2fmax。
2、奈奎斯特采样频率(1)概述奈奎斯特采样定理：要使连续信号采样后能够不失真还原，采样频率必须大于信号最咼频率的两倍(即奈奎斯特频率)。
奈奎斯特频率(Nyquist frequency )是离散信号系统采样频率的一半，因哈里•奈奎斯特(Harry Nyquist )或奈奎斯特—香农采样定理得名。
采样定理指出，只要离散系统的奈奎斯特频率高于被采样信号的最高频率或带宽，就可以真实的还原被测信号。
反之，会因为频谱混叠而不能真实还原被测信号。
采样定理指出，只要离散系统的奈奎斯特频率高于采样信号的最高频率或带宽，就可以避免混叠现象。
从理论上说，即使奈奎斯特频率恰好大于信号带宽，也足以通过信号的采样重建原信号。
但是，重建信号的过程需要以一个低通滤波器或者带通滤波器将在奈奎斯特频率之上的高频分量全部滤除，同时还要保证原信号中频率在奈奎斯特频率以下的分量不发生畸变，而这是不可能实现的。
在实际应用中，为了保证抗混叠滤波器的性能，接近奈奎斯特频率的分量在采样和信号重建的过程中可能会发生畸变。
因此信号带宽通常会略小于奈奎斯特频率，具体的情况要看所使用的滤波器的性能。
需要注意的是，奈奎斯特频率必须严格大于信号包含的最高频率。
如果信号中包含的最高频率恰好为奈奎斯特频率，那么在这个频率分量上的米样会因为相位模糊而有无穷多种该频率的正弦波对应于离散采样，因此不足以重建为原来的连续时间信号。
(2) 奈奎斯特频率的应用除了奈奎斯特频率之外，还有一个指标非常重要，这个指标就是测量装置的带宽。
严格讲，带宽包含上限和下限两个数值，但是，由于许多宽频带的测量设备，比如说变频功率分析仪，其带宽的频率上限远远大于频率下限，或者频率下限为零，因此，一般以频率上限作为该仪器的带宽。
一般而言，带宽指3db带宽。
-3db带宽并不表明高于带宽上限频率的信号不能通过测量仪器。
举例而言，某功率分析仪的带宽上限为100kHz,那么，100kHz的正弦波通过测量仪器的AD转换器之前的电路时，幅值衰减为原信号幅值的70.7%,功率衰减为原信号的50%此外，对于非正弦波形，其含有的谐波频率高于信号频率(基波频率)。
因此，不能简单的认为，100kHz带宽的仪器可以用于测量100kHz的正弦波, 更不能认为100kHz带宽的仪器可以用于测量100kHz的方波或畸变波形。
要让采样过程符合奈奎斯特采样定理，测量仪器的带宽应该小于奈奎斯特频率。
若测量仪器的电路固有带宽高于奈奎斯特频率，应该在AD转换器之间加上截至频率小于奈奎斯特频率的防混叠滤波器。
对于后者，防混叠滤波器的截至频率就是仪器的带宽。
3、稀疏采样目前，Can des, Romberg Tao和Do noho等人提出了一种全新的理论一压缩感知理论(Compressed Sensing)。
该理论是一种崭新的信号采样、信号编码和信号解码理论。
采样速率不再像Nyquist速率一样，与信号的带宽密切相关,而是与信息在信号中的结构和位置息息相关。
编码过程是围绕观测器即观测矩阵展开的，而解码过程是一个优化计算过程。
该理论已经被证明能够用较低采样速率准确的进行信号采样，并且能够以很高的概率重构原始信号。
目前国内已经有科研单位的学者对其展开研究。
如我们学校课题组基于该理论提出采用超低速率采样检测超宽带回波信号。
其CS理论如图：稀疏采样，也被称为压缩感知、压缩传感或压缩采样，是一种利用稀疏的或可压缩的信号进行信号重构的技术。
或者可以说是信号在采样的同时被压缩，从而在很大程度上降低了采样率。
稀疏采样跳过了采集N个样本这一步骤，直接获得压缩的信号的表示。
其理论利用到了许多自然信号在特定的基上具有紧凑的表示。
即这些信号是“稀疏”的或“可压缩”的。
由于这一特性，稀疏采样理论的信号编解码框架和传统的压缩过程大不一样，主要包括信号的稀疏表示、编码测量和重构算法等三个方面。
简单地说，压缩感知理论指出：只要信号是可压缩的或在某个变换域是稀疏的，那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上，然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原信号，可以证明这样的投影包含了重构信号的足够信息。
在该理论框架下，采样速率不再取决于信号的带宽，而在很大程度上取决于两个基本准则：稀疏性和非相干性，或者稀疏性和等距约束性。
显然，在压缩感知理论中，图像/信号的采样和压缩同时以低速率进行，使传感器的采样和计算成本大大降低，而信号的恢复过程是一个优化计算的过程•因此，该理论指出了将模拟信号直接采样压缩为数字形式的有效途径。
从理论上讲任何信号都具有可压缩性，只要能找到其相应的稀疏表示空间，就可以有效地进行压缩采样。
当前，压缩感知理论主要涉及三个核心问题：(1) 具有稀疏表示能力的过完备字典设计；(2) 满足非相干性或等距约束性准则的测量矩阵设计；(3) 快速鲁棒的信号重建算法设计。
压缩感知理论必将给信号采样方法带来一次新的革命。
这一理论的引人之处还在于它对应用科学的许多领域具有重要的影响，如统计学、信息论、编码等。
目前，学者们已经在模拟-信息采样、合成孔径雷达成像、遥感成像、核磁共振成像、深空探测成像、无线传感器网络、信源编码、人脸识别、语音识别、探地雷达成像等诸多领域对压缩感知展开了广泛的应用研究。
Rice大学已经成功设计出了一种基于压缩感知的新型单像素相机，在实践中为取代传统相机迈出了实质性的一步。
(1) 压缩感知理论框架传统的信号采集、编解码过程如图所示：编码端先对信号进行采样，再对所有采样值进行变换，并将其中重要系数的幅度和位置进行编码，最后将编码值进行存储或传输：信号的解码过程仅仅是编码的逆过程，接收的信号经解压缩、反变换后得到恢复信号。
采用这种传统的编解码方法，由于信号的采样速率不得低于信号带宽的2倍，使得硬件系统面临着很大的采样速率的压力。
此外在压缩编码过程中，大量变换计算得到的小系数被丢弃，造成了数据计算和内存资源的浪费。
压缩感知理论对信号的采样、压缩编码发生在同一个步骤，利用信号的稀 疏性，以远低于Nyquist 采样率的速率对信号进行非自适应的测量编码。
测量 值并非信号本身，而是从高维到低维的投影值，从数学角度看，每个测量值是 传统理论下的每个样本信号的组合函数，即一个测量值已经包含了所有样本信 号的少量信息。
解码过程不是编码的简单逆过程，而是在盲源分离中的求逆思 想下。
利用信号稀疏分解中已有的重构方法在概率意义上实现信号的精确重构 或者一定误差下的近似重构。
解码所需测量值的数目远小于传统理论下的样本 数。
压缩感知理论的编解码框图(2) 压缩感知的基本理论及核心问题假设有一信号f(f R N )，长度为N ，基向量为i (i 12，…,N)，对信号进 行变换：N f a i i 或 fi 1显然f 是信号在时域的表示，是信号在 域的表示。
信号是否具有稀疏 性或者近似稀疏性是运用压缩感知理论的关键问题，若式中的 只有K 个是非 零值(N K)者仅经排序后按指数级衰减并趋近于零，可认为信号是稀疏的。
信号的可稀疏表示是压缩感知的先验条件。
在已知信号是可压缩的前提下，压 缩感知过程可分为两步：设计一个与变换基不相关的M N(MN)维测量矩阵对信号进行观测,得到M 维的测量向量。
由M维的测量向量重构信号(3) 信号的稀疏表示稀疏的数学定义：信号X在正交基下的变换系数向量为T x，假如对于° P 2和R 0，这些系数满足：II li p ( I i l P)1/P Ri则说明系数向量在某种意义下是稀疏的。
给出另一种定义：如果变换系数i X，i的支撑域｛i; i °｝的势小于等于K ,则可以说信号X是K项稀疏。
如何找到信号最佳的稀疏域？这是压缩感知理论应用的基础和前提，只有选择合适的基表示信号才能保证信号的稀疏度，从而保证信号的恢复精度。
在研究信号的稀疏表示时，可以通过变换系数衰减速度来衡量变换基的稀疏表示能力。
Can des和Tao研究表明，满足具有幕次(power-law)速度衰减的信号，可利用压缩感知理论得到恢复。
最近几年，对稀疏表示研究的另一个热点是信号在冗余字典下的稀疏分 解•这是一种全新的信号表示理论：用超完备的冗余函数库取代基函数，称之 为冗余字典，字典中的元素被称为原子•字典的选择应尽可能好地符合被逼近 信号的结构，其构成可以没有任何限制•从从冗余字典中找到具有最佳线性组 合的K 项原子来表示一个信号，称作信号的稀疏逼近或高度非线性逼近。
目前信号在冗余字典下的稀疏表示的研究集中在两个方面：(1) 如何构造一个适合某一类信号的冗余字典；(2) 如何设计快速有效的稀疏分解算法。
这两个问题也一直是该领域研究的热点，学者们对此已做了一些探索，其 中以非相干字典为基础的一系列理论证明得到了进一步改进•西安电子科技大 学的石光明教授也对稀疏表示问题进行了认真研究，并基于多组正交基级联而 成的冗余字典提出一种新的稀疏分解方法。
关于布莱克曼窗函数的研究1. 布莱克曼窗布莱克曼窗的时域形式可表示为：w( n) [0.42 0.5(cos ) 0.08cos( )]R N 1 N 1 N(n) (2-3)它的频域特性为：其中W R ()为矩形窗函数的幅度频率特性。
增加一个二次谐波余弦分量，可进一步降低旁瓣，但主瓣宽度进一步增加, 12为 N 。
加N 可减少过渡带。
布莱克曼窗函数的最大旁瓣之比主瓣值低 但是主瓣宽度是矩形窗函数的主瓣宽度的三倍。
布莱克曼窗主瓣宽，旁瓣小， 频率识别精度最低，但幅值识别精度最咼。
W ) 0.42W R ( ) 0.25[W R ( N 1) N 1]0.04[W RN 1) W R ( N 1)] (2-4)57db ,一个理想数字滤波器的频率响应为 Hd(ej co ),对应的时域序列为滤波器的 单位脉冲响应hd(n)，是无限长非因果的。