Matlab课程实验报告.docx
《Matlab课程实验报告.docx》由会员分享，可在线阅读，更多相关《Matlab课程实验报告.docx（57页珍藏版）》请在冰豆网上搜索。
Matlab课程实验报告
西北工业大学
《基于MATLAB的数字信号处理》实验报告
学院：
学 号：
姓 名：
专业：
实验时间：
2014年7月
指导教师：
西北工业大学
2014年07月
实验一MATLAB基本编程实验
一、实验目的及要求
1．熟悉MATLAB运行环境；
2.掌握MATLAB的基本语法和函数；
3.掌握MATLAB的基本绘图功能
二、实验设备（环境）及要求
1．计算机
2．Matlab软件编程实验平台
三、实验内容与步骤
1．求下列线性方程组的解
2.分别用for或while循环结构编写程序，求出
。
并考虑一种避免循环语句的程序设计算法实现同样的运算。
3.在同一坐标系下绘制以下3条曲线,并作标记。
四、设计思想
设计思想1.1：
对于AX=B类型的方程组求解，需要先对非奇异矩阵A求逆，若存在逆矩阵，然后与B相乘即为所求
设计思想1.2：
用for或while循环结构编写程序时，首先对K初始化为0,确定循环的起始条件、步长和终值条件；用非循环语句时，先确定首项、公比，然后带入到等比数列求和公式中
设计思想1.3：
首先确定绘制曲线的每两点之间的间隔，然后在一幅图上绘制出所有曲线，并进行标注
五、程序代码及注释
程序代码1：
>>A=[634;
-257;
8-4-3];%矩阵A
>>X=inv（A）
X=
0.1300-0.07000.0100
0.5000-0.5000-0.5000
-0.32000.48000.3600
>>B=[3-4-7]';
X=X*B
程序代码2：
方法1：
for循环结构实现
%结果初始化0
K=0.0;
fori=1:
106%默认步长为1
K=K+sqrt（3）/（2^i）;
end;
fprintf（'K=%f\n',K）;
方法2：
while循环结构实现
K=0.0;
i=1;
whilei<107
K=K+sqrt（3）/（2^i）;
i=i+1;
end;
fprintf（'K=%f\n',K）;
方法3：
不使用循环，利用等比数列的前106项和公式实现
K=sqrt（3）/2.0;%数列首项
K=K*（1-0.5^106）/（1-0.5）;%公比为0.5
fprintf（'K=%f\n',K）;
程序代码3：
x=0:
0.01:
pi;%步长0.01
y1=sin（x）;
y2=sin（x）.*sin（10*x）;
y3=-cos（x）;
gridon%显示网格
holdon%在一幅图中绘制多条曲线
plot（x,y1,'r','LineWidth',2）;%y1红线表示，线宽为2
plot（x,y2,'g','LineWidth',2）;%y2绿线表示，线宽为2
plot（x,y3,'b','LineWidth',2）;%y3蓝线表示，线宽为2
%x轴范围0~3.5，y轴范围-1.2~1.2
axis（[0,3.5,-1.2,1.2]）;
%图例注释自动放在最佳位置
legend（'sin（x）','sin（x）sin（10x）','-cos（x）'）;
xlabel（'x坐标轴'）;%标注横坐标
ylabel（'y坐标轴'）;%标注纵坐标
title（'三曲线对比图'）;%标注标题
六、实验结果
实验结果1：
X=
0.6000
7.0000
-5.4000
所以，X=[
]=[0.6000,7.0000,-5.4000]’
实验结果2：
K=1.732051
实验结果3：
三曲线的对比图如下所示：
图1.1三曲线的对比
实验二基于MATLAB信号处理实验
一、实验目的及要求
1.回顾数字信号处理的主要内容；
2.掌握利用MATLAB进行信号处理的方法；
3.了解信号处理工具箱中一些函数的功能；
二、实验设备（环境）及要求
1．计算机
2．Matlab软件编程实验平台
三、实验内容
1．设序列
其中，
2.一数字滤波器由频率响应函数描述
1）写出其差分方程表示；
2）画出上面滤波器的幅频和相频图；
3）产生信号的200个样本，通过该滤波器得到输出y（n）,试将输出y（n）的稳态部分与x（n）作比较，说明这两个正弦信号的幅度和相位是如何受该滤波器影响的。
3.设计FIR滤波器（带通滤波器，四种窗函数，fs＝2000Hz,带通频率300Hz—600Hz，阶数自选,画频率特性并分析比较）。
四、设计思想
设计思想2.1：
将求互相关问题转化为求卷积的问题，互相关的区间大小为两信号的区间端点分别相加，区间长度为两信号区间长度之和
设计思想2.2:
由系统的差分方程
，得到方程的各个项的系数，进而可求出系统的幅频和相频响应
设计思想2.3：
首先选定带通滤波器的阶数，根据采样频率、上、下限截止频率求得滤波器的各个参数，然后得出窗函数的频率特性
五、程序代码及注释
程序代码2.1：
clc;
clearall;
x=[3,1,7,0,-1,4,2];%冲击信号x的幅度
nx=[-3:
3];%冲击信号x的取值范围
%对x进行右移2个单位，得到x（n-2）
ny=nx+2;
y=x;
w=normrnd（0,1,[1,7]）;%w为均值为0，方差为1的高斯随机序列
nw=ny;%得到y（n）=x（n-2）+w（n）
%对x进行翻折
x=fliplr（x）;
nx=-fliplr（nx）;
%卷积
nyb=ny
（1）+nx
（1）;
nye=ny（length（y））+nx（length（x））;
nrxy=[nyb:
nye];
rxy=conv（y,x）;
subplot（1,1,1）;
stem（nrxy,rxy）;
axis（[-5,9,0,90]）;
xlabel（'x'）;
title（'x（n）与y（n）的互相关'）;
程序代码2.2.2：
clc
clearall
fs=1000;
b=[10001];
a=[1000-.8145];
[h,f]=freqz（b,a,512,fs）;
mag=abs（h）;%幅度
ph=angle（h）;%相位
subplot（2,1,1）;
ph=ph*180/pi;%由弧度转换为角度
plot（f,mag）;
grid;
xlabel（'Frequency/Hz'）;
ylabel（'Magnitude'）;
title（'幅频响应'）;
subplot（2,1,2）;
plot（f,ph）;
grid;
xlabel（'Frequency/Hz'）;
ylabel（'Phase'）;
title（'相频响应'）;
程序代码2.2.3：
clc
clearall
N=200;
n=linspace（-100,100,N）;
x=sin（pi*n/2）+5*cos（pi*n）;
N_fft=2^nextpow2（2*N）;
w=linspace（0,2*pi,N_fft）;
h_fft=（1+exp（-1j*4*w））./（1-0.8145*exp（-1j*4*w））;
x_fft=fft（x,N_fft）;
y_fft=x_fft.*h_fft;
y_temp=fftshift（ifft（y_fft））;
y=y_temp（N_fft/2:
N_fft/2+N-1）;
figure;
plot（w,abs（h_fft）,'b','LineWidth',2）;
holdon;
plot（w,angle（h_fft）,'g','LineWidth',2）;
legend（'幅度','相位'）
figure;
plot（n,x,'b'）;
holdon;
plot（n,real（y）,'g'）;
legend（'x（n）','y（n）稳态部分'）
程序代码2.3：
fs=2000;%采样频率
fc1=300;%下限截止频率
fc2=600;%上限截止频率
N=31;%滤波器的阶数
wlp=fc1/（fs/2）;
whp=fc2/（fs/2）;
wn=[wlp,whp];;%滤波器归一化后的上下限截止频率
w1=boxcar（N）;%矩形窗的时域响应
w2=hanning（N）;%汉宁窗的时域响应
w3=hamming（N）;%海明窗的时域响应
w4=blackman（N）;%布莱克窗的时域响应
%用不同的窗函数设计N阶的滤波器
b1=fir1（N-1,wn,w1）;
b2=fir1（N-1,wn,w2）;
b3=fir1（N-1,wn,w3）;
b4=fir1（N-1,wn,w4）;
%求出滤波器的频率响应
[H1f1]=freqz（b1,1,512,fs）;
[H2f2]=freqz（b2,1,512,fs）;
[H3f3]=freqz（b3,1,512,fs）;
[H4f4]=freqz（b4,1,512,fs）;
figure;
subplot（2,1,1）;
plot（f1,20*log10（abs（H1）））;
xlabel（'频率/Hz'）;
ylabel（'振幅/dB'）;
title（'矩形窗的幅频特性'）;
gridon;
subplot（2,1,2）;
plot（f1,180/pi*unwrap（angle（H1）））;
xlabel（'频率/Hz'）;
ylabel（'相位'）;
title（'矩形窗的相频特性'）;
gridon;
figure;
subplot（2,1,1）;
plot（f2,20*log10（abs（H2）））;
xlabel（'频率/Hz'）;
ylabel（'振幅/dB'）;
title（'汉宁窗的幅频特性'）;
gridon;
subplot（2,1,2）;
plot（f2,180/pi*unwrap（angle（H2）））;
xlabel（'频率/Hz'）;
ylabel（'相位'）;
title（'汉宁窗的相频特性'）;
gridon;
figure;
subplot（2,1,1）;
plot（f3,20*log10（abs（H3）））;
xlabel（'频率/Hz'）;
ylabel（'振幅/dB'）;
title（'海明窗的幅频特性'）;
gridon;
subplot（2,1,2）;
plot（f3,180/pi*unwrap（angle（H3）））;
xlabel（'频率/Hz'）;
ylabel（'相位'）;
title（'海明窗的相频特性'）;
gridon;
figure;
subplot（2,1,1）;
plot（f4,20*log10（abs（H4）））;
xlabel（'频率/Hz'）;
ylabel（'振幅/dB'）;
title（'布莱克窗的幅频特性'）;
gridon;
subplot（2,1,2）;
plot（f4,180/pi*unwrap（angle（H4）））;
xlabel（'频率/Hz'）;
ylabel（'相位'）;
title（'布莱克窗的相频特性'）;
gridon;
六、实验结果
实验结果2.1：
图2-1x（n）与y（n）的互相关序列图
由实验结果可知，x（n）与y（n）的互相关只在区间[-4,8]上有能力，刚好是区间[-3,3]与右移后的区间[-1,5]两端点之和，与结论一致。
且互相关在2处达到最大。
实验结果2.2.1：
其表示的差分方程为：
y（n）-0.8145y（n-4）=x（n）+x（n-4）
实验结果2.2.2：
滤波器的幅频和相频图如下所示：
图2-2滤波器的幅频与相频图
实验结果2.2.3：
由下图实验结果可知，输出信号相对于输入信号有一小小的延迟，基本上x（n）的频点都通过了，滤波器是个梳状filter，正好在想通过的点附近相位为0，也就是附加延迟为0
图2-3滤波器的幅度和相位变化
图2-4两信号波形
实验结果2.3：
四种带通滤波器的窗函数的频率响应如下所示：
图2-5矩形窗的频率特性
图2-6汉宁窗的频率特性
图2-7海明窗的频率特性
图2-8布莱克曼窗的频率特性
实验三基于MATLAB的图像处理实验
一、实验目的及要求
1．了解图像处理的基本概念和功能；
2.掌握利用MATLAB进行图像处理的方法；
3.了解图像处理工具箱中一些函数的功能；
二、实验设备（环境）及要求
1．计算机
2．Matlab软件编程实验平台
三、实验内容
1．对一幅分别添加了椒盐噪声和高斯白噪声的图象进行均值和中值滤波（取三种不同大小的模板），并分析实验结果。
2．对一幅图像进行DCT变换，分别将小于5、10、20的变换系数置为0，进行反DCT得到重构图像，显示原图像和各重构图像并分别计算重构图像的峰值信噪比和置为0的变换系数的个数。
四、设计思想
设计思想3.1：
任意读取一幅图像，利用matlab自带函数imnoise（）向图像中加入椒盐噪声和高斯噪声，并分别用3*3、5*5、7*7模板进行均值、中值滤波处理。
均值滤波是一种线性滤波，也是低通滤波。
中值滤波是一种统计滤波器，是非线性的。
分别向图像中加入高斯噪声和椒盐噪声，利用不同的滤波方法，以及不同大小的窗口模板进行滤波，比较结果，得出结论。
由实验可以看出，中值滤波对椒盐噪声的处理效果好，均值滤波对高斯噪声的处理效果好。
无论是中值滤波还是均值滤波，都是其模板越大，图像越模糊。
所以选择模板大小时，以处理图片效果最佳为宜。
设计思想3.2：
由于DCT变换有使图像能量几种在左上方的特性，因此重构图像保留了原始图像大部分的图像特征，其视觉效果与原始图像相差不大。
对比重构前后的图像易知，重构后的图像稍显模糊，这是因为该压缩算法为有损压缩，压缩后的图像丢失了原始图像部分数据信息。
五、程序代码及注释
程序代码3.1：
%读取源图像
I=imread（'cameraman.tif'）;
%向图像中加入椒盐噪声
J=imnoise（I,'salt&pepper',0.1）;
%向图像中加入高斯噪声
K=imnoise（I,'gaussian',0.02）;
figure;
subplot（2,2,1）;
imshow（I）;
title（'原图像'）;
subplot（2,2,3）;
imshow（J）;
title（'加入椒盐噪声'）;
subplot（2,2,4）;
imshow（K）;
title（'加入高斯噪声'）;
%进行3*3均值滤波
J_average1=filter2（fspecial（'average',3）,J）/255;
%进行5*5均值滤波
J_average2=filter2（fspecial（'average',5）,J）/255;
%进行7*7均值滤波
J_average3=filter2（fspecial（'average',7）,J）/255;
figure;
subplot（2,2,1）;
imshow（J）;
title（'加入椒盐噪声'）;
subplot（2,2,2）;
imshow（J_average1）;
title（'3*3模板均值滤波'）;
subplot（2,2,3）;
imshow（J_average2）;
title（'5*5模板均值滤波'）;
subplot（2,2,4）;
imshow（J_average3）;
title（'7*7模板均值滤波'）;
%进行3*3中值滤波
J_Med1=medfilt2（J,[33]）;
%进行5*5中值滤波
J_Med2=medfilt2（J,[55]）;
%进行5*5中值滤波
J_Med3=medfilt2（J,[77]）;
figure;
subplot（2,2,1）;
imshow（J）;
title（'加入椒盐噪声'）;
subplot（2,2,2）;
imshow（J_Med1）;
title（'3*3模板中值滤波'）;
subplot（2,2,3）;
imshow（J_Med2）;
title（'5*5模板中值滤波'）;
subplot（2,2,4）;
imshow（J_Med3）;
title（'7*7模板中值滤波'）;
%进行3*3均值滤波
K_average1=filter2（fspecial（'average',3）,K）/255;
%进行5*5均值滤波
K_average2=filter2（fspecial（'average',5）,K）/255;
%进行7*7均值滤波
K_average3=filter2（fspecial（'average',7）,K）/255;
figure;
subplot（2,2,1）;
imshow（K）;
title（'加入高斯噪声'）;
subplot（2,2,2）;
imshow（K_average1）;
title（'3*3模板均值滤波'）;
subplot（2,2,3）;
imshow（K_average2）;
title（'5*5模板均值滤波'）;
subplot（2,2,4）;
imshow（K_average3）;
title（'7*7模板均值滤波'）;
%进行3*3中值滤波
K_Med1=medfilt2（K,[33]）;
%进行5*5中值滤波
K_Med2=medfilt2（K,[55]）;
%进行5*5中值滤波
K_Med3=medfilt2（K,[77]）;
figure;
subplot（2,2,1）;
imshow（K）;
title（'加入高斯噪声'）;
subplot（2,2,2）;
imshow（K_Med1）;
title（'3*3模板中值滤波'）;
subplot（2,2,3）;
imshow（K_Med2）;
title（'5*5模板中值滤波'）;
subplot（2,2,4）;
imshow（K_Med3）;
title（'7*7模板中值滤波'）;
程序代码3.2：
%读取源图像
I=imread（'cameraman.tif'）;
%%第2题离散余弦变换及反变换
%对图像进行离散余弦变换
J=dct2（I）;
%小于5的变换系数的索引值
[colrow]=size（find（abs（J）<5））;
A=col*row;%置为0的变换系数的个数
%小于10的变换系数的索引值
[colrow]=size（find（abs（J）<10））;
B=col*row;%置为0的变换系数的个数
%小于20的变换系数的索引值
[colrow]=size（find（abs（J）<20））;
C=col*row;%置为0的变换系数的个数
%将小于5的变换系数置为0后做离散余弦反变换
J（abs（J）<5）=0;I1=idct2（J）;
%将小于10的变换系数置为0后做离散余弦反变换
J（abs（J）<10）=0;I2=idct2（J）;
%将小于20的变换系数置为0后做离散余弦反变换
J（abs（J）<20）=0;I3=idct2（J）;
%显示原图及反变换结果
figure
（2）;
subplot（2,2,1）;
imshow（I）;
title（'原图'）;
subplot（2,2,2）;
imshow（I1,[0,255]）;
title（'系数置0的个数小于5'）;
subplot（2,2,3）;
imshow（I2,[0,255]）;
title（'系数置0的个数小于10'）;
subplot（2,2,4）;
imshow（I3,[0,255]）;
title（'系数置0的个数小于20'）;
%计算反重构时，DCT的变换系数的置0个数小于5时的峰值信噪比及置为0的变换系数的个数
I=double（I）;
I1=double（I1）;
[Row,Col]=size（I）;
[Row,Col]=size（I1）;
MSE1=sum（sum（（I-I1）.^2））/（Row*Col）;
PSNR1=10*log10（255^2/MSE1）;
fprintf（'图像的峰值信噪比：
MSE1=%f\n',MSE1）;
fprintf（'置为0的变换系数的个数为：
PSNR1=%f\n',PSNR1）;
%计算反重构时，DCT的变换系数的置0个数小于10时的峰值信噪比及置为0的变换系数的个数
I=double（I）;
I2=double（I2）;
[Row,Col]=size（I）;
[Row,Col]=size（I2）;
MSE2=sum（sum（（I-I2）.^2））/（Row*Col）;
PSNR2=10*log10（255^2/MSE2）;
fprintf（'图像的峰值信噪比：
MSE2=%f\n',MSE2）;
fprintf（'置为0的变换系数的个数为：
PSNR2=%f\n',PSNR2）;
%计算反重构时，DCT的变换系数的置0个数小于20时的峰值信噪比及置为0的变换系数的个数
I=double（I）;
I3=double（I3）;
[Row,Col]=size（I）;
[Row,Col]=size（I3）;
MSE3=sum（sum（（I-I3）.^2））/（Row*Col）;
PSNR3=10*log10（255^2/MSE3）;
fprintf（'图像的峰值信噪比：
MSE3=%f\n',MSE1）;
fprintf（'置为0的变换系数的个数为：
PSNR3=%f\n',PSNR3）;
六、实验结果
实验结果3.1：
由实验结果可以看出，对于椒盐噪声，中值滤波效果明显优于均值滤波，对于高斯噪声，上图中中值滤波比均值滤波效果略佳。
随着模板从3*3增加到7*7，对于均值滤波，滤噪效果变好，但图像变模糊；对于中值滤波，图像细节丢失变多，滤除高斯噪声的效果变好。
图3-1加噪前、后图像对比
图3-2加椒盐噪声的图像均值滤波前、后的图像对比
图3-3加椒盐噪声的图像中值滤波前、后的图像对比
图3-4加高斯噪声的图像均值滤波前、后的图像对比
图3-5加高斯噪声的图像中值滤波前、后的图像对比
实验结果3.2：
图3-6原图及重构图像
图3-7程序运行结果
由实验结果可知，当DCT变换的系数置0个数小于5时，重构图像与原图像的峰值信噪比为2.768259，重构图像置为0的变换系数的个数个数为：
43.708737；
当DCT变换的系数置0个数小于10时，重构图像与原图
升级会员 