一、           课程设计目的

 “信号与系统”是一门重要的专业基础课，MATLAB作为信号处理强有力的计算和分析工具是电子信息工程技术人员常用的重要工具之一。本课程设计基于MATLAB完成信号与系统综合设计实验，以提高学生的综合应用知识能力为目标，是“信号与系统”课程在实践教学环节上的必要补充。通过课设综合设计实验，激发学生理论课程学习兴趣，提高分析问题和解决问题的能力。

二、           课程设计时间

第十五、十六周。上机时间安排见附件一。第十六周周五提交课程设计报告并答辩。

三、           参考书目

1、  谷源涛、应启珩、郑君里著，信号与系统——MATLAB综合实验，北京：高等教育出版社，20##年1月。

2、  郑君里、应启珩 、杨为理，信号与系统引论，北京：高等教育出版社，20##年3月。

3、  梁虹等，信号与系统分析及Matlab实现，北京：电子工业出版社，20##年2月。

四、           注意事项

1、  基本部分，共三道题，每人都需要全部完成，要求十五周周五做完。

2、  提高部分，共八道题，每人按照学号分配（见附件二）只做其中的一题。

3、  第十六周周五所提交的课程设计报告如有雷同，一律退回重写。

一、基础部分

（一）傅里叶变换分析

1、周期信号的谱分析，要求任意给定单频周期信号，能够准确计算出其幅度谱和相位谱，并画出图形，要求正确显示频率。

%周期信号频谱分析

t=0:0.2:8;

x=cos(t);   %产生单频周期信号

subplot(311);

plot(x);

xlabel('t');

ylabel('cos(t)');

title('波形');

X=fft(x);   %计算DFS系数

a=abs(X); 

b=angle(X);  

subplot(312);

stem(a);   %画出幅度频谱

title('幅度谱');

xlabel('w');

ylabel('a');

subplot(313);

stem(b);   %画出相位频谱

title('相位谱');

xlabel('w');

ylabel('b');

图1

图1中是cos函数在一个周期内的波形、幅度谱和相位谱，由图可知余弦函数傅氏变换的幅度谱是两个冲击，相位谱为一条直线。

2、  非周期信号的频谱分析，要求分析语音信号的幅度谱和相位谱，并画出图形。

%非周期信号频谱分析

y=wavread('u');%读取音乐文件

Fn=fft(y);%作傅里叶变换

f=44000;

Fnf=abs(Fn); %幅度

Fnx=angle(Fn);%相位

N=length(y);

M=0:N-1;

figure;

subplot(3,1,1);

plot(M/7005000,y);%语音信号波形

title('语音信号的波形');

xlabel('t');

ylabel('y');

subplot(3,1,2);

stem(M/f,Fnf);%幅度谱

axis([0 0.0005 0 1500]);

title('幅度谱');

xlabel('Hz');

ylabel('Fnf');

subplot(3,1,3);

stem(M/f,Fnx);%相位谱

axis([0 0.0005 0 10]);

title('相位谱');

xlabel('Hz');

ylabel('Fnx');

图2

由图2可以看出，语音信号的幅度谱变化不大，相位谱是断续的。

3、对于方波，设计程序计算其傅里叶级数系数，仿真吉伯斯现象。

%仿真吉伯斯现象

t=0:0.01:20*pi;%定义时间长度

y=0;

num=input('叠加的正弦波个数=')

for k=1:num       %用for循环实现波形的叠加

    y=y+(4/pi)*(1/((k*2)-1))*sin(((k*2)-1)*t);

    r(k,:)=y;

end

figure;

plot(t,r);

title('吉伯斯现象');

xlabel('t');

ylabel('f(t)');

grid on;

图3

由图3可以发现，N个正弦波的叠加后的波形近似于方波。

（二）、连续时间系统分析：

1、自行设计一个有初始条件的微分方程，至少二阶。

2、解出其零输入解，并画出图形，与手工计算相比较。

3、解出其单位冲激响应，并画出图形，与手工计算相比较。

4、设定某一激励信号，用卷积方法解出其零状态响应，并画出图形，与手工计算相比较。

5、计算系统的幅频响应和相频响应，并画出图形。

a=[1,3,2];        %微分方程式左边的系数

b=[1,3];          %微分方程式右边系数

sys=tf(b,a);       %系统的传递函数

rsys=ss(sys);     %讲传递函数转换为状态空间

t=0:0.01:50;      %定义时间常数

x=0*t;   

e=[1,2];

rzi=lsim(rsys,x,t,e);   %零输入解

subplot(211);

plot(t,rzi);

xlabel('t');

ylabel('rzi(t)');

title('零输入解');

rzi1=4*exp(-t)-3*exp(-2*t);    %零输入解的手工计算

subplot(212);

plot(t,rzi1);

xlabel('t');

ylabel('rzi1(t)');

title('零输入解的手工计算');

 %单位冲激响应

rzs=impulse(sys,t);       %求冲激响应

figure;

subplot(211);

plot(rzs,'k:');

xlabel('t');

ylabel('rzs(t)');

title('单位冲激响应');

rzs0=2*exp(-t)-exp(-2*t);  %冲激响应的手工计算

subplot(212);

plot(rzs0,'k:');

xlabel('t');

ylabel('rzs1(t)');

title('单位冲激响应的手工计算');

x=exp(t);              %卷积求零状态响应 

r1=conv(x,rzs);

u=r1(1:5001);

figure;

subplot(211);

plot(r1);

xlabel('t');

ylabel('r1');

title('e^t的零状态响应');

r2=-exp(-t)+1/3*exp(-2*t)+2/3*exp(t); %零状态响应的手工计算

r3=u'-r2;

subplot(212);

plot(r3);

xlabel('t');

ylabel('r2');

title('e^t的零状态响应的手工计算');

F=fft(rzs);         %求系统函数

Ff=abs(F);          %系统函数的幅频响应

Fx=angle(F);        %系统函数的相频响应

figure;

subplot(211);

plot(Ff,'r-');

xlabel('w');

ylabel('Ff');

title('幅频响应');

subplot(212);

plot(Fx,'g-');

xlabel('w');

ylabel('Fx');

title('相频响应');

图4

图4是MATLAB计算的零输入解与手工计算的零输入解，比较可知，两者基本吻合。

图5

由图5可以看出微分方程的冲激响应为一条下降的曲线，即2*exp(-t)-exp(-2*t)。

图6

由图6，我们可以发现输入为exp(t)时的零状态响应由一段递增的曲线和一段递减的曲线构成，即-exp(-t)+1/3*exp(-2*t)+2/3*exp(t)。

图7

图7是系统函数的幅频响应和相频响应。

（三）、离散时间系统分析：

1、自行设计某离散时间系统函数，至少是二阶的系统，画出零极点图，判断系统的稳定性。

2、求出单位样值响应，并画出图形。

3、求出系统的幅频响应和相频响应，并画出图形。

b=[1,1];   %系统函数的分子

a=[1,0.2,-0.24];   %系统函数分母

[z,p,k]=tf2zp(b,a);

zplane(b,a);  %画零极点图

title('零极点图(稳定系统)');

xlabel('Rez');

ylabel('jImz');

figure;

h=impz(a,b,0:10);   %求单位样值响应

stem(h);         %画出单位样值响应

title('单位样值响应');

xlabel('n');

ylabel('h(n)');

figure;

y=freqz(a,b);    %求频谱

hf=abs(y);      %幅度谱

hx=angle(y);        %相位谱

subplot(211);

plot(hf,'k-');    %输出幅度谱

title('幅频响应');

xlabel('n');

ylabel('hf(n)');

subplot(212);

plot(hx,'g:');  %输出相位谱

title('相频响应');

xlabel('n');

ylabel('hx(n)');

图8

图8所示为零极点图，由图可知，系统函数的极点为z=0.4和z=-0.6，零点为z=-1 z=0。

图9

图9为系统的单位样值响应

图10

图10为系统函数的幅频响应和相频响应。

二、提高部分（频分复用）

1、自行给出二路语音信号，分别显示其频谱，并播放语音。

2、对二路语音信号进行频分复用，显示复用后的频谱，播放语音。

3、设计程序对频分复用的信号进行解调，显示解调结果，并回放语音。

N=10000;

n=[0:N-1];

bit=10000;

[e1,fs1,bits]=wavread('li.wav');

y1=e1(1:N);

[e2,fs2,bits]=wavread('xc.wav');%读取两路语音信号

y2=e2(1:N);

figure;

subplot(2,1,1);

plot(y1);   %输出a的时域波形

xlabel('t');

ylabel('y1');

title('a的波形');

sp1=max(abs(e1));

sound(e1,fs1,bits);

Fy1=fft(y1);

f=n/N*fs1;

subplot(2,1,2);

plot(f,abs(Fy1));%a的幅度谱

xlabel('t');

ylabel('Fy1');

title('a的幅度谱');    

bw1=fs1/2;

pause(10)

figure;

subplot(2,1,1);

plot(y2); %b的时域波形

xlabel('t');

ylabel('y2');

title('b的波形');

sp2=max(abs(e2));

sound(e2,fs2,bits);pause(2)

Fy2=fft(y2);

subplot(2,1,2);

plot(f,abs(Fy2));  %b的幅度谱

xlabel('f');

ylabel('Fy2');

title('b的幅度谱');

bw2=fs2/2;

%调制

fca=12000;

fcb=45000;

fs=120000;

Ts=1/fs;

fn=n/N*fs;

y1cos=cos(2*pi*fca.*n*Ts);%第一个载波信号

y2cos=cos(2*pi*fcb.*n*Ts);%第二个载波信号

y=y1'.*y1cos+y2'.*y2cos;%调制相加得到复合信号

F=fft(y);

figure;

subplot(211);

plot(y);

xlabel('t');

ylabel('y');

title('复合信号的波形');

subplot(212);

plot(fn,abs(F));

xlabel('fn');

ylabel('|F|');

title('复合信号的幅度谱');xlabel('HZ');

%设计滤波器

dat=zeros(1,N);

dat(round((fca-bw1)/fs*N:(fca+bw1)/fs*N))=1;%

dat(round((fs-fca-bw1)/fs*N:(fs-(fca-bw1))/fs*N))=1;

fi1=dat;     %a滤波器

dat=zeros(1,N);

dat(round((fcb-bw2)/fs*N:(fcb+bw2)/fs*N))=1;

dat(round((fs-(fcb-bw2)-2*bw2)/fs*N:(fs-(fcb-bw2))/fs*N))=1;

fi2=dat;    %b滤波器

figure;

subplot(211);

plot(fn,fi1);

axis([0 120000 0 1.5]);

title('滤波器fi1');

xlabel('HZ');

ylabel('fi1');

subplot(212);

plot(fn,fi2);

axis([0 120000 0 1.5]);

xlabel('HZ');

ylabel('fi2');

title('滤波器fi2');

xlabel('Hz');

fily1=F.*fi1; %滤波

fily2=F.*fi2;

%混频

hy1=ifft(fily1).*y1cos;

hy2=ifft(fily2).*y2cos;

figure;

subplot(211);

plot(abs(fft(hy1)));

subplot(212);

plot(abs(fft(hy2)));

%低通滤波器

filter=zeros(1,N);

filter(1:round(bw1/fs*N))=1;

filter(round((fs-bw1)/fs*N):N)=1;

filter1=filter;

filter(1:round(bw2/fs*N))=1;

filter(round((fs-bw2)/fs*N):N)=1;

filter2=filter;

%解调

Fya=fft(hy1).*filter1;

ya=ifft(Fya);

figure;

subplot(211);

plot(n,ya);

xlabel('t');

ylabel('ya');

title('a的解调');

sp3=max(abs(ya));

sound(ya/sp3,fs1,bits)

Fyb=fft(hy2).*filter2;

yb=ifft(Fyb);

subplot(212);

plot(n,yb);

title('b的解调');

xlabel('t');

ylabel('yb');

sp4=max(abs(yb));

sound(yb/sp4,fs1,bits)

图11

图11为语音信号a的波形和幅度谱，可知a的带宽为20K左右。

图12

图12为语音信号b的波形和幅度谱，可知其带宽跟a相差不大。但是幅度比a小。

图13

图示为a和b调制并复合后的波形和幅度谱，此时两路信号已经被调制到各自的载波频率左右。其中两边是a的幅度谱，中间为b的幅度谱。

图14

图14为两个带通滤波器，fi1用来滤出a信号，fi2用来滤出b信号。

图15

图15为解调后的混频信号，混频信号a是拿复合信号乘a的载波信号得到的频谱，混频信号b是拿复合信号b乘b的载波信号得到的频谱。

图16

图16为还原的a和b语音信号，是用低通滤波器对混频信号滤波后再进行傅里叶逆变换还原到时域，与图11和图12比较可知，还原后的信号跟原信号基本一致。

五、实验问题及体会：

   这次课程设计，尝到了各种滋味，有在困惑中的苦恼，在迷茫中的挣扎，也有在学习中的快乐，在成功中的兴奋，这是一个独立思考和挑战自己恒心的过程。实验中学到的不仅仅是MATLAB的应用和一些课题的解决方法，更重要的是锻炼了自己的意志，在做基础部分的时候，我在对MATLAB一无所知中苦苦摸索，一次一次地编写代码，试验函数的用法，慢慢地学会了怎么写一些简单的程序，实验过程中出现最多的问题是这个：

??? Error using ==> times

Matrix dimensions must agree.

这个问题让我特别困惑，后来发现时矩阵维度不一致，通过上网查找资料，才把问题解决。我发现课本知识在实践面前很脆弱，不是说课本知识不重要，我们在学好课本知识的同时更要注重联系实际，要能解决实际问题，把课本上学到的东西应用到课程设计里面来，比如说频分复用，频分复用就是课本上讲过的一个应用，但是具体到自己设计，就要考虑各种问题，比如说载波的选择、滤波器的设计，这些课本上只是提到但是怎么解决得靠自己想办法。这次课程设计对我的启发很大，我懂得了遇到困难首先要思考，查找解决办法，耐心分析错误原因，做事要有耐心，我会在以后的学习中注重实践。