用MATLAB实现数字信号处理

摘要

MATLAB(Matrix Laboratory-----矩阵实验室)是由美国Mathworks公司于1984年正式推出的一种以矩阵运算为基础的交互式程序语言。与其他的语言相比她有着“简洁和智能化，适应科技专业人员的思维方式和书写习惯”的优点，所以为科技人员所乐于接受。更因为它可适用于各种平台，并且随计算机硬、软件的更新而及时升级。所以MATLAB逐渐成为当前我国大学教学和科学研究中一种非常重要的工具软件。

首先通过的MATLAB和DSP(Digital Signal Processing)------数字信号处理的简单概括，在第二章里利用MATLAB中几个简单而又强大的命令（比如：plot,函数的赋值，M函数等）来对数字信号处理中信号的取样和还原进行模拟。本文先从理论上说明：原函数与冲激函数相乘得到的函数fs(t)其频谱是原函数频谱的无限个频移项组成，其幅值是原函数的1/Ts.而且，只有当取样频率fs大于原函数fm的两倍时，那么相后的频谱就不会混叠了，也就是说，取样函数（fs=f(t)*s(t)）包含了原函数的全部信息。那么我们就可以在接收端用一个低通滤波器（其频率响应的幅度是Ts）,就可以把原信号还原出来了。反之，如果采样的时候，采样频率小于2*fm,那么采样后得到的函数的频谱会出现混叠。这样，采样函数就不能包括原函数的全部信息，使得接收断的低通滤波器所截取的信号不能真实的反映原函数。

关键词：MATLAB,DSP,取样频率，信号最高频率

Abstract

MATLAB(Matrix Laboratory)was a scientist computer language which was invented by the American company Mathworks in1984.This language is based on Matrix https://www.wendangku.net/doc/f715475372.html,pared with other computer science languages,it has the advantages in its clearance in the expression of its language and fits the thinking method of many scientists.So,as the days goes by,it becomes an acceptable and enjoyable language for the scientists.For the most important aspects,it could be installed in many different operate systems and be updated its version according to the development of computer software and hardware,thus,MATLAB becomes an important teaching tools in many China university.

In the first chapter,this article provides us an overall description of the MATLAB and DSP (Digital Signal Processing).Then,in the second chapter,this article use some simple but powerful commands to simulate the experiment of the sampling and reconstruction.The reader will find in the second chapter,the article first provides the theory why the signal could be reconstructed: when the sampling frequency(fs)is larger or equal than the two times of the maximum frequency of the signal(fm).Then,the sampling signal could include the whole information of the original signal and when we use a low pass filter in the end,we could reconstruction the original signal.

Otherwise,if the fs<=2fm,the sampling signal could not provide us the whole information of the original signal,thus based on this we could not reconstruct the original signal.

Key word:MATLAB,DSP,sampling frequency，the maximum frequency of the signal

第一章绪论

第一节MATLAB语言简介

1.1.1MATLAB语言概述

MATLAB(Matrix Laboratory-----矩阵实验室)是一种科学计算软件，适用于工程应用各领域的分析设计与复杂计算，它使用方便，输入简捷运算高效且内容丰富，很容易由用户自行扩展。因此，当前以成为美国和其他发达国家大学教学和科学研究中最常用且必不可少的工具。

MATLAB是由美国Mathworks公司于1984年正式推出的，到1988年有了3.x（DOS）版本；1992年推出4.x（Windows）版本；1997年推出5.1（Windows）版本。以后又升级到5.3（也称R11）版本。2000年下半年，Mathworks公司推出了他们的最新产品MATLAB6.0（R12）试用版，并于2001年初推出了正式版。随着版本的升级，内容不断扩充，功能更加强大。另一方面对使用环境也提出了更高的要求。近几年来，Mathworks公司在将MATLAB 语言运用于系统仿真和实时运行方面，取得了很多成绩，更扩大了它的应用前景。

MATLAB是“矩阵实验室”（MATrix LABoratoy）的缩写，它是一种以矩阵运算为基础的交互式程序语言，着重针对科学运算，工程计算和绘图的需求。与其他计算机语言相比，其特点是简洁和智能化，适应科技专业人员的思维方式和书写习惯，使得编程和调试效率大大提高。它用解释方式工作，键入程序立即得出结果，人机交互性能好，为科技人员所乐于接受。特别是它可适用于各种平台，并且随计算机硬、软件的更新而及时升级。因此，MATLAB 语言在国外的大学工学院中，特别是数值计算用的最频繁的电子信息类学科中，已成为每个学生都应掌握的工具了。它大大提高了课程教学、解题作业、分析研究的效率。

1.1.2MATLAB语言的特点

MATLAB语言有以下特点：

1．起点高

(1)每个变量代表一个矩阵。当前的科学计算中，几乎无处不用矩阵运算，这使它的优势得到充分的体现。在MATLAB中，每个变量代表一个矩阵，它可以有n*m个元素。

(2)每个元素都看作复数，这个特点在其他语言中是不多见的。

(3)所有的运算都对矩阵和复数有效，包括加、减、乘、除、函数运算等。

2．人机界面适合科技人员。

(1)MATLAB的语言规则和笔算式相似。

(2)矩阵的行列无需定义。要输入一个矩阵，用其他语言时必须先定义矩阵的阶数，而

MATLAB则不必有阶数定义语句。

(3)键入算式立即得结果，无需编译。MATLAB是以解释方式工作的，即它对每条语句解释后立即执行，若有错误也立即做出反应，便于编程者立即改正。

3.强大而简易的做图功能。

(1)能根据输入数据自动确定坐标绘图。

(2)能规定多种坐标系。

(3)能绘制三维坐标中的曲线和曲面。

(4)可设置不同颜色、线型、视角等。

4.智能化程度高。

(1)绘图时自动选择最佳坐标，大大方便了用户。

(2)做数值积分时自动按精度选择步长。

(3)自动检测和显示程序错误的能力强，易于调试。

5.功能丰富，可扩展性能强。

扩展部分称为工具箱。它实际上是用MATLAB的基本语句编成的各种子程序集，用于解决某一方面的专门问题，或实现某一类的新算法。现在已经有控制系统、信号处理、图象处理、系统辩识、模糊集合、神经元网络及小波分析等工具箱。

第二节数字信号处理简介

1.2.1数字信号处理概述

数字信号处理，或者说对信号的数字处理，是60年代前后发展起来的一门新兴学科，进入70年代以来，随着电子计算机、大规模集成电路（LSI）和超大规模集成电路（VLSI），以及微处理器技术的迅速发展、数字信号处理无论在理论上还是在工程应用中，都是目前发展最快的学科之一，而且日趋完善和成熟。

目前对数字信号处理没有一个严格的定义，我们不妨把用包括软件在内的数字技术来处理模拟信号叫做数字信号处理。或者简单的说，数字信号处理就是用数字的方法对信号进行变换和处理。

从学科的内容来看，数字滤波器和快速傅立叶变换是数字信号处理的两个最重要的学科分支，也是数字信号处理的核心内容。事实上，对信号进行调制、滤波、变换、增强、压缩、检测、估计、识别等，都是数字信号处理的研究范围。用数字的方法实现具有各种功能的信号处理处理系统，可以使系统做到更加经济和小型化，进而有可能实现过去用模拟技术难以实现的各种功能。

图1给出了了模拟信号处理数字化处理系统的简化框图。此系统先将模拟信号变换为数

字信号，经数字信号处理后，在变换成模拟信号输出，其中抗混叠滤波器的作用是将输入信号x(t)中高于折叠频率（其值等于采样频率的一半）的分量滤除，以防信号频谱的混叠。随后，信号经采样和A/D 变换后，变成数字信号x(n),数字信号处理器对x(n)进行处理，得到输出数字信号，经D/A 变换器变成模拟信号。此信号经低通滤波器，滤除不需要的高频分量，然后输出平滑的模拟信号y(t).

1.2.2信号是如何被处理的？

我们在实际中所遇到的信号大部分是模拟信号。它们的幅度随时间连续地变化。通常用由无源或者有源器件所组成的电网络对他们作处理。这中方法称为模拟信号处理。例如无线电和电视接收机。

模拟信号：x(t)y(t)

它们也可以用包括加法器、乘法器和逻辑元件或用微处理器组成的数字硬件来处理。但是必须先把模拟信号变成适合于数字硬件来处理的形式。这样的信号形式成为数字信号。它在一定的时间间隔内只取有限数目的数值之一，因而它能用二进制或位来表示。数字信号的处理称为DSP,它可以用框图来表示。

图中各个方框的意义介绍如下：

PrF:

这是一个前置滤波器或反泄漏滤波器，它限制了模拟信号以防止泄漏。ADC:模-数转换器，它是由模拟信号产生出一个二进制流。

数字信号处理器：这是DSP 的心脏，它可以代表一个通用计算机或专用处理器，或数字硬件等等。

DAC:这是ADC 的逆运算，称为数-模转换器，它由一个二进制流产生一个阶梯波形，作为Prf

ADC 数字DSP 数字

DAC

Pof 模拟模拟模拟信号处理器抗混叠

滤波器A/D 数字信号处理D/A 低通滤波

X(x)Y(n)Y(t)

z(t)

形成模拟信号的第一步。

PoF:这是一个后置滤波器，它把阶梯波形平滑成为预期的模拟信号。

1.2.3DSP胜于ASP的优点

ASP的主要缺点是在做复杂信号处理时只有有限的能力。这造成了处理的不灵活性和系统时间的复杂性。这也造成了产品的昂贵。另一方面，采用DSP的方法，它可以把一个便宜的个人计算机变成一个强大的信号处理器。DSP的一些重要的优点如下：

(1)用DSP方法的系统的开发可以用通用计算机上的软件来进行。因此DSP比较容易开发和测试，而软件是可移植的。

(2)DSP运算是单纯的基于加法和乘法，这导致特别稳定的处理性能——例如不受温度影响。

(3)DSP运算可以方便地作实时的修改，通常只要改变程序，或者对寄存器重新加载。

(4)由于超大规模集成电路（VLSI）降低了寄存器、门、微处理器的价格，DSP比较便宜。

DSP的主要缺点是运算速度低，特别是对高频的信号。主要是由于它的上述有限，DSP 现在已经是很多技术的首选方法，比如消费电子、通信、无线电话和医学图象等。

第二章用MATLAB实现信号处理

*注：本文中的所有程序均在matlab 6.5下调试运行

信号经过采样和量化运算后变成离散信号。这些离散信号由数字信号处理器处理后，再用一个重构运算后变成模拟信号。用傅立叶分析可以从频域的观点来描述采样运算，分析其效果，然后重新把信号还原出来。我们假定量化的等级足够的多因而离散信号量化误差可以忽略不计。

第一节信号的取样：

取样定理论述了在一定的条件下,一个连续时间的信号完全可以用该信号在等时间间隔上的瞬时值表示,这些样本值完全包含了该连续时间信号的全部信息,利用这些样本值可以恢复原信号.可以说,取样定理在连续时间信号和离散时间信号之间架起了一座桥梁.可以说取样定理为连续时间信号与离散时间信号的相互转换提供了理论依据.

所谓的”取样”就是利用取样脉冲序列s(t)从连续时间信号f(t)中”抽取”一系列离散样本值的过程.这样得到的离散信号称为取样信号.如果图 2.1所示:

图2.1取样的模型

式中取样脉冲序列s(t)也称开关函数。如果其各脉冲间隔的时间相同，均为Ts,就称为均匀取样。Ts 称为取样周期，fs=1/Ts 称为取样频率或取样率。Ws=2*pi*fs=2*pi/Ts.

如果令f(t)和s(t)的傅立叶变换为F(jw)和S(jw),则由频域卷积定理，得取样信号fs(t)的频谱函数:

Fs(jw)=(1/2*pi)*F(jw)*S(jw)(式2.1)

2.1.1冲激取样：

如果取样序列s(t)是周期为Ts 的冲激函数序列x(t),则称为冲激取样。且冲激序列x(t)的频谱函数也是周期序列，即：

F[s(t)]=F(x(t))=F[∑x(t-n*Ts)]=ws ∑x(w-n*ws)

其中n 从-∞到+∞(式 2.2)

把(式2.2)代入（式2.1）中，可以得到：

Fs(jw)=(1/2*pi)*F(jw)*ws*∑x(w-nws)=1/Ts*∑F(j(w-nws))

其中n 从-∞到+∞(式 2.3)

由（式 2.3）可以知道，取样信号fs(t)的频谱由原信号频谱F(jw)的无限个频移项组成，其频谱的角频率为nws,其幅值为原频谱的1/Ts。2.1.2矩形脉冲取样：

如果取样脉冲序列s(t)是幅度为1，脉宽为L (L

S(jw)=P(jw)=F(p(t))=(2*pi*L)/Ts *∑Sa(n*ws*L/2)x(w-nws)(式2.4)

上式代入式（式2.1）可以得到取样信号fs(t)的频谱函数为：

Fs(jw)=1/(2*pi)*F(jw)*(2*pi*L/Ts)*∑Sa(nwsL/2)*x(w-nws)=L/T*∑Sa(nws*L/2)*F

[j(w-nws)]（式 2.5）

比较式2.5和式 2.1可以发现，经过冲激取样或矩形脉冲取样后，其取样信号fs(t)的频谱很相似。因此当ws>2wm 时取样信号的频谱Fs(jw)也不会出现混叠，从而能从取样信号fs(t)中恢复信号f(t)。f(t)s(t)

fs(t)

可以看出，如果ws>2wm,那么各相邻频谱不会发生互相重叠，这时候就设法（比如用低通滤波器），从取样信号的频谱Fs(jw)中得到原信号的频谱，即从取样信号fs(t)中恢复原信号f(t),如果ws<2wm，那么频移后的各相邻频谱将互相重叠，如图所示。这样就无法将它们分开，因而也就不能再恢复原信号。频谱重叠的这种现象叫做混叠现象。可见，为了不发生混叠现象，ws>2wm。

第二节信号的重构

从采样定理和上述的例子可以清楚的看到，如果对有限带宽信号xa(t)以高于2fc的频率进行采样，就可以从其采样序列x(n)重构原模拟信号。重构可以考虑为两个步骤：

1.先把样本集转换为一个加权脉冲串列。

2.然后把这个脉冲串列通过一个带宽为[-Fs/2,Fs/2]的低通滤波器进行滤波。

现在以冲激取样为例，研究如何从取样信号fs(t)恢复信号f(t)。

设有冲激信号fs(t),其取样角频率ws>2wm(wm为原信号的最高频率)。Fs(jw)为fs(t)的频谱函数，为了从Fs(jw)中无失真地恢复F(jw),选择一个理想的低通滤波器，其

频率响应的幅度为Ts,截止角频率为wc(wm

Ts,|w|

H(jw)={0,|w|>wc

因此：F(jw)=Fs(jw)*H(jw)（式 2.6）

即恢复了原信号的频谱函数F(jw)。

根据卷积定理，（式2.6）的相应于时域为：

f(t)=fs(t)*h(t)(式2.7)

由于冲激取样信号

fs(t)=f(t)·s(t)=f(t)·∑x(t-n*Ts)=∑f(n*Ts)x(t-nTs)(式2.8)其中n从-∞到+∞

为简便，选wc=ws/2,则Ts=2*pi/ws=pi/wc，得：

h(t)=Sa(ws*t/2)(式 2.9)

把（式 2.8）和（式 2.9）代入（式 2.7）中可以得到：

f(t)=∑

f(nTs)x(t-nTs)*Sa(wst/2)=∑f(nTs)Sa[ws/2*(t-nTs)]=∑f(nTs)Sa(ws*t/2-n*pi)

上式表明，连续信号f(t)可以展开成正交取样函数（Sa函数）的无穷级数，该级数的系数等于取样值f(nTs)。也就是说，若在取样信号fs(t)的每个样点出，画一个最大峰值为f(nTs)的波形函数，那么其合成波形就是原信号f(t).因此，只要以知各取样值f(nTs),就能唯一的确定出原信号f(t)。

2.2.1时域取样定理

一个频谱在区间（-wm,wm）以外为零的频带有限信号f(t),可以由其样点值f(nTs)确定。需要注意的是，为了能从取样信号fs(t)中恢复原信号f(t),需满足两个条件：（1）f(t)必须是带限信号，其频谱函数|w|>wm各处为零；（2）取样频率不能过低，必须fs>2fm,或者说取样间隔不能太长，必须Ts<1/(2fm),否则将发生混叠。

第三节MATLAB对信号重建和还原的实现

2.3.1离散傅立叶变换

因为信号的采样中需要用到离散傅立叶变换，故这里单独拿出来，并且将它编成一个子程序。严格的说MATLAB不能直接处理无限长连续的信号，但是根据离散傅立叶变换的两个性质：周期性和对称性，我们时间间隔dt或者频率间隔dw取的很小，然后来估计表达式在[0,pi]上，或者一个时间段里的特性。

假定x(n)是有限长的，我们在[0，pi]间按等间隔频点来估计X(exp(jw)),则可以用矩阵—向量相乘的运算来实现。假设

x(n)在[n1,nN]之间，且wk=(pi/M)*k,k=0,1,……M

那么X(exp(jwk))=∑exp(-j*(pi/M)*k*nl)x(nl),其中l范围为[1，N]

在MATLAB中，我们把序列和下标排成行向量，得到：

X=x*exp(-j*w*n’)

编写一个离散傅立叶变换子程序dtft.m:

function X=dtft(x,w)

X=x*exp(-j*[1:length(x)’*w]);

下面我们以函数x(t)=A*exp(-a*t)*sin(b*t)（其中a=b=50*sqrt(2)*pi）为例来模拟信号的取样和还原，并对上面的理论进行验证。

2.3.2采样

同样，因为MATLAB不能直接处理无限长连续的信号，我们把dt取的很小，dt<

Xa(jΩ)=∑xc(m)e(-j*Ωmdt)dt，在这个程序中，分别设定四个采样频率：fs=10kHz, 1kHz,400Hz,200Hz.

程序如下：

clear all;

close all;

fs=10000;fs1=1000;fs2=400;fs3=200;

t=0:1/fs:0.1;

A=444.128;a=50*sqrt(2)*pi;b=a;

xa=A*exp(-a.*t).*sin(b.*t);

k=0:511;f=fs*k/512;

Xa=dtft(xa,2*pi*k/512);

T1=1/fs1;t1=0:T1:0.1;

x1=A*exp(-a.*t1).*sin(b.*t1);

X1=dtft(x1,2*pi*k/512);

T2=1/fs2;t2=0:T2:0.1;

x2=A*exp(-a.*t2).*sin(b.*t2);

X2=dtft(x2,2*pi*k/512);

T3=1/fs3;t3=0:T3:0.1;

x3=A*exp(-a.*t3).*sin(b.*t3);

X3=dtft(x3,2*pi*k/512);

figure(1);

subplot(4,2,1);plot(t,xa)

subplot(4,2,3);stem(t1,x1,'.');

subplot(4,2,5);stem(t2,x2,'.');

subplot(4,2,7);stem(t3,x3,'.');

subplot(4,2,2);plot(f/pi,abs(Xa)/max(abs(Xa)));

subplot(4,2,4);plot(f/pi,abs(X1)/max(abs(X1)));

subplot(4,2,6);plot(f/pi,abs(X2)/max(abs(X2)));

subplot(4,2,8);plot(f/pi,abs(X3)/max(abs(X3)));

结果为

-200

200

00.050.100.5100.5100.5101000200030004000

0.51上图中，自上到下的四个图分别为模拟信号，采样频率分别为1kHz,400Hz,200Hz 的图以及他们的幅度谱。可以看出，fs=1kHz 的失真较小，fs=400Hz 的失真较大，而fs=200Hz 失真比较严重。

2.3.3重建

根据前面的公式xa(t)=∑x(n)sinc[Fs(t-nTs)];设定一个ti 的值求xa(ti)的问题，可归结为一个行向量x(n)和一个同长的由n’构成的列向量g(ti-n’T)相乘。这与前面的dtft 的思想比较类似。“t-n*T ”这个的实现可以把t 设成行向量，nT’为列向量。我们的目的是把它构成一个行数与n 一样，列数与t 一样的矩阵。因此可以把t 右乘ones(length(n),1);把nT 左乘向量ones(1,length(t))就可以了。

可以令TNT=ones(length(n),1)t-n’*T*ones(1,length(t))

程序如下：

clear all;

close all;

A=444.128;a=50*sqrt(2)*pi;b=a;

F1=400;T1=1/F1;dt1=T1/3;

tp=0.03;

t1=0:dt1:tp;

n1=0:tp/T1;

TMN=ones(length(n1),1)*t1-n1'*T1*ones(1,length(t1));

x1=A*exp(-a*n1*T1).*sin(b*n1*T1);

xa1=x1*sinc(F1*TMN);

F2=1000;T2=1/F2;dt2=T2/3;

tp=0.03;

t2=0:dt2:tp;

n2=0:tp/T2;

TMN=ones(length(n2),1)*t2-n2'*T2*ones(1,length(t2));

x2=A*exp(-a*n2*T2).*sin(b*n2*T2);

xa2=x2*sinc(F2*TMN);

t=0:1/800000:0.03

x0=A*exp(-a*t).*sin(b*t);

subplot (1,3,1);stem (t1,xa1,'.')

subplot (1,3,2);stem (t2,xa2,'.')

subplot (1,3,3);plot (t,x0);

结果如下：

00.0050.010.0150.020.0250.03

00.0050.010.0150.020.0250.03

00.0050.010.0150.020.0250.03-1000

100

200

其中自上到下的三个图为Fs=400Hz,Fs=1000Hz 的还原图以及模拟信号的初始图形。可以看出，Fs=1000Hz 的信号还原的比较好。

第四章总结

通过前面对MATLAB 和DSP 的一些总体介绍，以及文章的第三部分在理论上证明出信号经过采样以后能够还原成原始信号，并且运用MATLAB 这个强有力的软件对这一过程（采样—重建）的模拟和证明。加深了我对采样原理的理解，牢牢的记住了采样原理的重要前提条件： 1.信号的带宽是有限的，即(-fm,fm); 2.采样信号的频率fs 必须大于2*fm。

可以说，MATLAB 是一个功能强大的软件，但是她又是非常的简单易学。而这篇文章所用到的命令，仅仅是她的冰山一角：像MATLAB 的工具箱，画图等非常有用的功能还没有用

到，或者说只用到了一个非常小的部分。

当然，这也告诉我自己所学的只是很小的一部分，将来要更加努力的学习。

参考文献

[1].Ragael C.Gonzalez,Richard E.Woods,Digital Image Processing;,1987

[2].吴大正，信号与线性系统分析，高等教育出版社，1985

[3].Augustine H.Gray,John D.Markel,"Distance Measures for speech processing",IEEE,Trans. Acoust,Speech,Signal,Processing,Vol.Assp-24,p.p.380-391,October1976.

[4].L.R.Rabiner.R.W.schafer,"Digital Processing of Speech Signals,1978.

[5].陈国通，数字通信原理，哈尔滨工业大学出版社，1996

[6].John R Deller,JR,John G.Proakis,John H.L.Hansen,Discrete-Time Processing Of Speech Signals,276-292,(1990)

[7].John R Deller,JR、John G.Proakis、John H.L.Hansen,Discrete-Time Processing Of Speech Signals,276-292,(1990)

[8].Matlab User Manual"The Math Works Inc.1993

MATLAB上机指导书

MATLAB上机指导书 电子信息科学与技术专业 张焕明孙明编 佛山科学技术学院 2005年9月

目录 前言 实验一 MATLAB基础知识 1 实验二矩阵与数组 5 实验三基本操作命令 8 实验四高级操作命令 10 实验五 MATLAB的M函数编程 12

前言 MATLAB的名称源自Matrix Laboratory,是一门计算语言，它专门以矩阵的形式处理数据．MATLAB将计算与可视化集成到一个灵活的计算机环境中，并提供了大内置函数,可以在广泛的工程问题中直接利用这些函数获得数值解．此外，用MATLAB 编写程序，犹如在一张草稿纸上排列公式和求解问题一样效率高，因此被称为“演算纸式的”科学工程算法语言．在我们高等数学的学习过程中，可以结合 MATLAB 软件，做一些简单的编程应用，在一定程度上弥补我们常规教学的不足，同时，这也是我们探索高职高专数学课程改革迈出的一步．

实验一 MATLAB 基础知识 一、实验目的 1、MATLAB 的使用初步练习 2、MATLAB 的窗口组成 二、实验内容 1、掌握表达式的输入方法 2、MATLAB 的常量及其表示方法 3、分号、百分比号、逗号及省略号的用法 4、向量和矩阵的处理方式；常用的数学函数；搜索路径的概念；MATLAB 的帮 助功能。 三、实验仪器、设备和材料 1、微型计算机，能正常运行Matlab 6.0或以上版本 2、Matlab6.0或以上版本 四、实验原理 略（参考教材的相关部分） 五、实验步骤 1、MATLAB 文件的编辑、存储和执行 MATLAB 提供了两种运行方式，即命令行和M 文件方式． A ．命令行方式 直接在命令窗口输入命令来实现计算或作图功能． 例如，若要求表达式 9 .248.26107 sin 369.12÷?+π的值，我们可在MATLAB 命令窗口中键入下面的命令： >> 1.369^2+sin(7/10*pi)*sqrt(26.48)/2.9 （回车） 观测运行结果并解释原因 也可将计算的结果赋给某一个变量，例如输入 : >> a=1.369^2+sin(7/10*pi)*sqrt(26.48)/2.9 （回车） 观测运行结果并解释原因 B ．M 文件的运行方式 1）文件编辑 在MATLAB 窗口中单击File 菜单依次选择NewM-File,打开M 文件输入运行界面，如下图所示。此时屏幕上会出现所需的窗口，在该窗口中输入程序文件，可以进行调试和运行．与命令行方式相比，M 文件方式的优点是可以调试，可重复应用． 2）文件存储 单击File 菜单，选择Save 选项，可将自己所编写的程序存在一个后缀为m 的文件中． 3）运行程序 在M 文件窗口中选择Debug 菜单中的run 选项，即可运行此M 文件；也可在MATLAB 命令窗口中直接输入所要执行的文件名后回车即可．但需要的是该程序文件必须存在MATLAB 默认的路径下．用户可以在MATLAB 窗口中单击File 菜单选择Set Path 将要执行的文件所在的路径添加到MATLAB 默认的路径序列中． 2、MATLAB 基本运算符及表达式 表1-1 基本运算符

数字信号处理Matlab实现实例(推荐给学生)

数字信号处理Matlab 实现实例 第１章离散时间信号与系统 例1-1 用MATLAB计算序列{-2 0 1 –1 3}和序列{1 2 0 -1}的离散卷积。 解 MATLAB程序如下： a=[-2 0 1 -1 3]; b=[1 2 0 -1]; c=conv(a,b); M=length(c)-1; n=0:1:M; stem(n,c); xlabel('n'); ylabel('幅度'); 图1.1给出了卷积结果的图形，求得的结果存放在数组c中为：{-2 -4 1 3 1 5 1 -3}。 例1-2 用MATLAB计算差分方程 当输入序列为时的输出结果。 解 MATLAB程序如下： N=41; a=[0.8 -0.44 0.36 0.22]; b=[1 0.7 -0.45 -0.6]; x=[1 zeros(1,N-1)];

k=0:1:N-1; y=filter(a,b,x); stem(k,y) xlabel('n');ylabel('幅度') 图 1.2 给出了该差分方程的前41个样点的输出，即该系统的单位脉冲响应。 例1-3 用MATLAB 计算例1-2差分方程 所对应的系统函数的DTFT 。 解 例1-2差分方程所对应的系统函数为： 123 123 0.80.440.360.02()10.70.450.6z z z H z z z z -------++= +-- 其DTFT 为 23230.80.440.360.02()10.70.450.6j j j j j j j e e e H e e e e ωωωω ωωω--------++= +-- 用MATLAB 计算的程序如下： k=256; num=[0.8 -0.44 0.36 0.02]; den=[1 0.7 -0.45 -0.6]; w=0:pi/k:pi; h=freqz(num,den,w); subplot(2,2,1); plot(w/pi,real(h));grid title('实部') xlabel('\omega/\pi');ylabel('幅度')

基于Matlab的数字水印设计——基于DCT域的水印实现

摘要 数字水印（Digital Watermark）技术是指用信号处理的方法在数字化的多媒体数据中嵌入隐蔽的标记，这种标记通常是不可见的，只有通过专用的检测器或阅读器才能提取。数字水印是信息隐藏技术的一个重要研究方向。随着数字水印技术的发展，数字水印的应用领域也得到了扩展，数字水印的基本应用领域是版权保护、隐藏标识、认证和安全不可见通信。 当数字水印应用于版权保护时，潜在的应用市场在于电子商务、在线或离线地分发多媒体内容以及大规模的广播服务。数字水印用于隐藏标识时，可在医学、制图、数字成像、数字图像监控、多媒体索引和基于内容的检索等领域得到应用。数字水印的认证方面主要ID卡、信用卡、ATM卡等上面数字水印的安全不可见通信将在国防和情报部门得到广泛的应用。 本文主要是根据所学的数字图象处理知识，在MATLAB环境下，通过系统编程的方式，建立并实现基于DCT域的数字水印加密系统。该系统主要包含数字水印的嵌入与提取，仿真结果表明，数字水印算法具有有效性、可靠性、抗攻击性、鲁棒性和不可见性，能够为数字媒体信息在防伪、防篡改、认证、保障数据安全和完整性等方面提供有效的技术保障。 关键词：数字水印；MATLAB；DCT

目录 1 课程设计目的 (1) 2 课程设计要求 (2) 3 数字水印技术基本原理 (3) 3.1 数字水印基本框架 (3) 3.2 算法分类 (3) 3.2.1 DCT法 (4) 3.2.2 其他方法 (4) 3.3 实际需要考虑的问题 (4) 3.3.1 不可见性 (4) 3.3.2 鲁棒性 (5) 3.3.3 水印容量 (5) 3.3.4 安全性 (5) 4 基于DCT变换仿真 (6) 4.1 算法原理 (6) 4.1.1 准备工作 (6) 4.1.2 选取8*8变换块 (7) 4.1.3 边界自适应 (7) 4.1.4 DCT变换与嵌入 (7) 4.1.5 恢复空域 (8) 4.2 嵌入算法扩展 (8) 4.2.1 RGB彩色图像三个矩阵的划分 (8) 4.2.2 八色彩色水印 (8) 4.3 水印的提取 (9) 4.4 仿真程序 (9) 5 结果分析 (14) 结束语 (16) 参考文献 (17)

MATLAB上机实验(答案)

MATLAB工具软件实验(1) （1）生成一个4×4的随机矩阵，求该矩阵的特征值和特征向量。程序： A=rand(4) [L,D]=eig(A) 结果： A = 0.9501 0.8913 0.8214 0.9218 0.2311 0.7621 0.4447 0.7382 0.6068 0.4565 0.6154 0.1763 0.4860 0.0185 0.7919 0.4057 L = -0.7412 -0.2729 - 0.1338i -0.2729 + 0.1338i -0.5413 -0.3955 -0.2609 - 0.4421i -0.2609 + 0.4421i 0.5416 -0.4062 -0.0833 + 0.4672i -0.0833 - 0.4672i 0.4276 -0.3595 0.6472 0.6472 -0.4804 D = 2.3230 0 0 0 0 0.0914 + 0.4586i 0 0 0 0 0.0914 - 0.4586i 0 0 0 0 0.2275 （2）给出一系列的a值，采用函数 22 22 1 25 x y a a += - 画一组椭圆。 程序： a=0.5:0.5:4.5; % a的绝对值不能大于5 t=[0:pi/50:2*pi]'; % 用参数t表示椭圆方程 X=cos(t)*a; Y=sin(t)*sqrt(25-a.^2); plot(X,Y) 结果： （3）X=[9,2,-3,-6,7,-2,1,7,4,-6,8,4,0,-2], （a）写出计算其负元素个数的程序。程序： X=[9,2,-3,-6,7,-2,1,7,4,-6,8,4,0,-2]; L=X<0; A=sum(L) 结果： A =

数字信号处理MATLAB中FFT实现

MATLAB中FFT的使用方法 说明：以下资源来源于《数字信号处理的MATLAB实现》万永革主编 一.调用方法 X=FFT(x)； X=FFT(x，N)； x=IFFT(X); x=IFFT(X,N) 用MATLAB进行谱分析时注意： （1）函数FFT返回值的数据结构具有对称性。 例： N=8; n=0:N-1; xn=[43267890]; Xk=fft(xn) → Xk= 39.0000-10.7782+6.2929i0-5.0000i 4.7782-7.7071i 5.0000 4.7782+7.7071i0+5.0000i-10.7782-6.2929i Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。 （2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。 二.FFT应用举例 例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。

clf; fs=100;N=128;%采样频率和数据点数 n=0:N-1;t=n/fs;%时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);%信号 y=fft(x,N);%对信号进行快速Fourier变换 mag=abs(y);%求得Fourier变换后的振幅 f=n*fs/N;%频率序列 subplot(2,2,1),plot(f,mag);%绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2));%绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);%信号 y=fft(x,N);%对信号进行快速Fourier变换 mag=abs(y);%求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag);%绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2));%绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果：

数字水印技术及基于MATLAB的快速实现

2011年3月刊计算机工程应用技术信息与电脑 China Computer&Communication 1. 引言 多媒体及网络的迅速发展使得多媒体信息的交流和传输变得更 加简单和快捷，然而，这也使盗版者能以低廉的成本复制及传播未经 授权的数字产品，这种对数字产品保护和信息安全的迫切需求，导致 了数字水印技术成为多媒体信息安全领域的一个热点问题。数字水印 技术是在不影响宿主媒体主观质量的情况下，在宿主媒体（文本、图 像、视频、音频）中嵌入不易被人察觉的标识信息，用以证明原创作 者对其作品的所有权，并作为鉴定、起诉非法侵权的证据。 2. 数字水印的特征 一般认为数字水印应具有以下特征： (1) 安全性。数字水印应该是安全、难以被篡改的。当数字作品 发生变化时，数字水印应当也相应发生变化；同时，未经授权的个人 不得修改水印，理论上是未经授权的用户不能检测到产品中是否含有 水印。 (2) 鲁棒性。当被保护的数据在经过攻击后，嵌入的水印信息仍 能保持好的完整性并能以一定的正确概率被检测到。这些可能的攻击 包括噪声、滤波、剪切、旋转和编码等。 (3) 不可感知性。数字水印的嵌入不应使得原始作品发生可以感 知的变化，也不能使得被保护数据在质量上发生可以感觉到的失真。 (4) 可证明性。在多媒体作品的实际应用过程中可能需要多次加 入水印，这时水印系统必须能够允许水印被多次嵌入到被保护的数 据，而且每个水印均能独立地被证明。 (5) 无歧义性。恢复出的水印或对水印判决结果能够表明版权的 惟一，不会发生多重版权纠纷问题。 3. 数字水印的基本原理 通用的水印技术包含两个方面：水印的嵌入和水印的提取或检 测，如图1和图2所示。 图1 水印信号嵌入 图2 水印信号提取或检测 4. 数字水印的研究现状 4.1 文本水印 文本水印就是将代表著作人身份的信息(水印)嵌入到电子出版物 中，在产生版权纠纷时来验证版权的归属。其主要分为三大类：基于 文档结构的水印方法、基于自然语言处理技术的水印方法、基于传统 图像的水印方法。 基于文档结构的各种水印方法都只是提留在文本的表层，无法抵 抗对于文本结构和格式的攻击，简单的重新录入攻击就能使之失效， 因此这些水印方法普遍存在鲁棒性差的缺点。自然语言文本水印方法 相对提高了抗攻击的能力，但普遍存在容量不足的问题。基于传统图 像的文本水印普遍存在鲁棒性不高、操作复杂的缺点。 4.2 图像水印 根据水印的实现过程，图像水印算法可分为空域算法和变换域算 法。空域算法是通过直接改变原始图像的像素值来嵌入水印，通常具 有较快的速度，但鲁棒性差，且水印容量也会受到限制；变换域算法 是通过改变某些变换系数来嵌入水印，通常具有很好的鲁棒性和不可 见性。其实现一般是基于图像变换，如DCT、DFT、DWT等。重点介 绍一下变换域算法。 4.2.1 离散傅里叶变换 (DFT) 该方法是利用图像的DFT来嵌入信息。通信理论中调相信号的抗 干扰能力比调幅信号的抗干扰能力强，同样在图像中利用相位信息嵌 入的水印也比用幅值信息嵌入的水印更稳健。实验表明该方法的抗压 缩能力比较弱。 4.2.2 离散余弦变换 (DCT) DCT能把空间域的图像转换到变换域上进行研究，从而能很容易 了解到图像的各空间频域成分，进行相应处理。基于DCT的水印方法 与基于DFT的水印方法相比有较好的鲁棒性，但是无法做到对图像信 号内容的自适应，因此往往会造成对图像特征的明显损害，不可感知 性不是最佳。 4.2.3 离散小波变换 (DWT) DWT是一种时间---频率信号的多分辨率分析方法，在时频两域 都具有表征信号局部特征的能力。实验表明，与DCT、DFT变换相比 较，基于DWT的水印算法的鲁棒性最优，且与JPEG2000、MPEG4压 缩标准兼容，利用DWT产生的水印具有良好的视觉效果和抵抗多种 攻击的能力，且不可感知性最好。 4.3 音频水印 音频水印利用音频文件的冗余信息和人耳听觉系统的特点来嵌入 水印，其可以保护声音数字产品不被随意复制和篡改，如CD唱片， 广播电台的节目内容等。音频水印的三种基本方法：扩频嵌入方 法、回声隐藏方法和相位编码方法。 4.4 视频水印 视频水印是通过对视频载体的时间和空间冗余来嵌入水印，其既 不影响视频质量，又能达到保护节目制作者的合法权益和控制数字产 品的复制。视频水印从算法要求上同图像水印有许多相似之处，但视 频水印也有一些独特之处，如能够在压缩和未压缩的格式下实时完成 水印的检测，对MPEG压缩、A/D和D/A转换等都有较好的稳健性。 数字水印技术涉及到通信理论、编码理论、噪声理论、视听觉 感知理论、扩频技术、信号处理技术、数字图像处理技术、多媒体技 术、模式识别技术、算法设计等理论，用到经典的DFT离散傅立叶变数字水印技术及基于MATLAB的快速实现 张 巍1 时宏伟2 （1.78179部队，四川成都 610011；2. 川大智胜，四川成都 610045） 摘要：数字水印是近几年来出现的数字产品版权保护技术，是当前国际学术界的研究热点．该文论述了数字水印的提出及研究现状、水印的基本原理和算法、水印的分类等情况，并介绍了一种可以快速上手的高效的实用语言——MATLAB，同时给出了一个用MATLAB工具在静止图像上嵌入水印的实例。 关键词：数字水印；MATLAB；DCT 中图分类号：TP39 文献标识码：A 文章编号：1003-9767（2011）03-0130-02

南京理工大学数字信号处理matlab上机完美版

1.已知3阶椭圆IIR数字低通滤波器的性能指标为：通带截止频率0.4π，通带波纹为0.6dB，最小阻带衰减为32dB。设计一个6阶全通滤波器对其通带的群延时进行均衡。绘制低通滤波器和级联滤波器的群延时。 %Q1_solution %ellip(N,Ap,Ast,Wp) %N--->The order of the filter %Ap-->ripple in the passband %Ast->a stopband Rs dB down from the peak value in the passband %Wp-->the passband width [be,ae]=ellip(3,0.6,32,0.4); hellip=dfilt.df2(be,ae); f=0:0.001:0.4; g=grpdelay(hellip,f,2); g1=max(g)-g; [b,a,tau]=iirgrpdelay(6,f,[0 0.4],g1); hallpass=dfilt.df2(b,a); hoverall=cascade(hallpass,hellip); hFVT=fvtool([hellip,hoverall]); set(hFVT,'Filter',[hellip,hoverall]); legend(hFVT,'Lowpass Elliptic filter','Compensated filter'); clear; [num1,den1]=ellip(3,0.6,32,0.4); [GdH,w]=grpdelay(num1,den1,512); plot(w/pi,GdH); grid xlabel('\omega/\pi'); ylabel('Group delay, samples'); F=0:0.001:0.4; g=grpdelay(num1,den1,F,2); % Equalize the passband Gd=max(g)-g; % Design the allpass delay equalizer [num2,den2]=iirgrpdelay(6,F,[0,0.4],Gd); [GdA,w] = grpdelay(num2,den2,512); hold on; plot(w/pi,GdH+GdA,'r');

实验一 基于Matlab的数字信号处理基本

实验一 基于Matlab 的数字信号处理基本操作 一、 实验目的：学会运用MA TLAB 表示的常用离散时间信号；学会运用MA TLAB 实现离 散时间信号的基本运算。 二、 实验仪器：电脑一台，MATLAB6.5或更高级版本软件一套。 三、 实验内容： （一） 离散时间信号在MATLAB 中的表示 离散时间信号是指在离散时刻才有定义的信号，简称离散信号，或者序列。离散序列通常用)(n x 来表示，自变量必须是整数。 离散时间信号的波形绘制在MATLAB 中一般用stem 函数。stem 函数的基本用法和plot 函数一样，它绘制的波形图的每个样本点上有一个小圆圈，默认是空心的。如果要实心，需使用参数“fill ”、“filled ”，或者参数“.”。由于MATLAB 中矩阵元素的个数有限，所以MA TLAB 只能表示一定时间范围内有限长度的序列；而对于无限序列，也只能在一定时间范围内表示出来。类似于连续时间信号，离散时间信号也有一些典型的离散时间信号。 1. 单位取样序列 单位取样序列)(n δ，也称为单位冲激序列，定义为 ) 0() 0(0 1)(≠=?? ?=n n n δ 要注意，单位冲激序列不是单位冲激函数的简单离散抽样，它在n =0处是取确定的值1。在MATLAB 中，冲激序列可以通过编写以下的impDT .m 文件来实现，即 function y=impDT(n) y=(n==0); %当参数为0时冲激为1,否则为0 调用该函数时n 必须为整数或整数向量。 【实例1-1】 利用MATLAB 的impDT 函数绘出单位冲激序列的波形图。 解：MATLAB 源程序为 >>n=-3:3; >>x=impDT(n); >>stem(n,x,'fill'),xlabel('n'),grid on >>title('单位冲激序列') >>axis([-3 3 -0.1 1.1]) 程序运行结果如图1-1所示。 图1-1 单位冲激序列

基于MATLAB的数字水印算法实现

数字水印作为一门新的学科, 自 1993 年 Tirkel 等人正式提出到现在十几年里, 国内外对数字水印的研究都引起了极大的关注, 从最初的版权保护, 已扩展到多媒体技术, 广播监听, in-ternet 等多个领域。数字水印是永久镶嵌在其他数据( 主要指宿主数据) 中具有可鉴别性的数字信号或数字模式, 其存在不能影响宿主数据的正常使用。为了使数字水印技术达到一定的设计要求, 当前水印数据一般应具备不可感知性(imperceptible) 、鲁棒性(Robust) 、可证明性、自恢复性和安全保密性等特点。在数字水印技术中, 水印的数据量和鲁棒性构成了一对基本矛盾。理想的水印算法应该既能隐藏大量数据, 又可以抗各种信道噪声和信号变形。然而在实际中, 这两个指标往往不能同时实现, 实际应用往往只偏重其中的一个方面。如果是为了隐蔽通信, 数据量显然是最重要的, 由于通信方式极为隐蔽, 遭遇敌方篡改攻击的可能性很小, 因而对鲁棒性要求较为不高。但对保证数据安全来说, 情况恰恰相反, 各种保密的数据随时面临着被盗取和篡改的危险, 对鲁棒性的要求很高, 而对隐藏数据量的要求则居于次要地位。典型的数字水印系统至少包含两个组成部分- - 水印嵌入单元和水印检测与提取单元。将水印信息进行预处理后加入到载体中, 称为嵌入。从水印化数据中提取出水印信息或者检测水印信息的存在性称为水印的提取和检测。数字水印算法主要

是指水印的嵌入算法, 而提取算法往往被看成是嵌入算法的逆变换。 当前典型的嵌入算法主要被分为空间域水印算法和变换域水印算法。DCT 变换域算法是数字水印算法的典型代表, 也是数字水印中较为常用的一种稳健的算法。其算法思想是选择二值化灰度图像作为水印信息, 根据水印图像的二值性来选择不同的嵌入系数, 并将载体图像 ( 原始图像) 进行 8×8 的分块, 再将灰度载体图像( 原始图像) 进行 DCT变换。然后, 将数字水印信息的灰度值直接植入到载体灰度图像的 DCT 变换域中, 实现水印的嵌入。而后, 将嵌入了水印信息灰度图像进行 IDCT( 逆离散的余弦变换) 变换, 得到含有了嵌入水印信息的图像, 嵌入过程完毕。水印的提取、检测过程为嵌入过程的逆过程, 其方法和嵌入方法有所雷同不再进行介绍。 下面以 MATLAB 为工具, 给出一个在频域嵌入和提取黑白二值水印图像的实现过程。(1) 水印图像的预处理: 将水印信息图像进行灰度处理, 然后再将转换后的图像进行二值转换。而这些都是为了提高水印信息的安全性对图像所做的处理。(2) 读取原始公开图像(大小为 256×256) 和黑白水印图像(大小为 32×32, 模式为灰度) 到二维数组 I 和 J。(3) 将原始公开图像I 分割为互不覆盖的图像块, 每块大小为 8×8, 共分为 32×32 块。然后对分割后的每个小块Block- dct(x,y) 进行 DCT 变换, 得到变换后的小块 Block-dct(x, y)。(4) 取黑白水印图像中的一个元素 J(p, q) , 通过嵌入算法嵌入到原始公开图像块的中频系数中。(5) 对嵌入水印信息后的图像块Block- dct (x, y) 进行逆DCT 变换, 得到图像块 Block(x′, y′)。

青岛理工大学临沂年数字信号处理及MATLAB试卷

A卷

一、[15分] 1、10 2、f>=2fh

3、()()()y n x n h n =* 4、1 -az -11a 或者-z z ，a 1 -z 或1-1-az -1z 5、对称性 、 可约性 、 周期性 6、191点，256 7、典范型、级联型、并联型 8、T ω = Ω，)2 tan(2ω T = Ω或)2arctan(2T Ω=ω。 二、[20分] 1、C 2、 A 3、 C 4、C 5、B 6、D 7、B 8、A 9、D 10、A (CACCB DBADA) 三、[15分] 1、（5分） 混叠失真：不满足抽样定理的要求。 改善方法：增加记录长度 频谱泄漏：对时域截短，使频谱变宽拖尾，称为泄漏 改善方法：1）增加w (n )长度 2）缓慢截短 栅栏效应：DFT 只计算离散点（基频F0的整数倍处）的频谱，而不是连续函数。 改善方法：增加频域抽样点数N （时域补零），使谱线更密 2、（5分） 3、 （5分） IIR 滤波器: 1）系统的单位抽样相应h (n )无限长 2）系统函数H (z )在有限z 平面（ ）上有极点存在 3）存在输出到输入的反馈，递归型结构 Fir 滤波器: ? 1）系统的单位冲激响应h (n )在有限个n 处不为零； ? 2）系统函数 在||0 z >处收敛，在 处只有零点，即有限z 平面只有零点，而全部极点都在z =0处； ? 3）机构上主要是非递归结构，没有输入到输出的反馈，但有些结构中也包含有反馈的递归部分。 四、计算题（40分） 1、（12分）解： 解： 对上式两边取Z 变换，得： ()H z ||0z >

MATLAB上机答案

一熟悉Matlab工作环境 1、熟悉Matlab的5个基本窗口 思考题： （1）变量如何声明，变量名须遵守什么规则、是否区分大小写。 答：变量一般不需事先对变量的数据类型进行声明，系统会依据变量被赋值的类型自动进行类型识别，也就是说变量可以直接赋值而不用提前声明。变量名要遵守以下几条规则： 变量名必须以字母开头，只能由字母、数字或下划线组成。 变量名区分大小写。 变量名不能超过63个字符。 关键字不能作为变量名。 最好不要用特殊常量作为变量名。 （2）试说明分号、逗号、冒号的用法。 分号：分隔不想显示计算结果的各语句；矩阵行与行的分隔符。 逗号：分隔欲显示计算结果的各语句；变量分隔符；矩阵一行中各元素间的分隔符。 冒号：用于生成一维数值数组；表示一维数组的全部元素或多维数组某一维的全部元素。 （3）linspace（）称为“线性等分”函数，说明它的用法。 LINSPACE Linearly spaced vector.线性等分函数 LINSPACE(X1,X2)generates a row vector of100linearly equally spaced points between X1and X2. 以X1为首元素，X2为末元素平均生成100个元素的行向量。 LINSPACE(X1,X2,N)generates N points between X1and X2. For N<2,LINSPACE returns X2. 以X1为首元素，X2为末元素平均生成n个元素的行向量。如果n＜2，返回X2。 Class support for inputs X1,X2: float:double,single 数据类型：单精度、双精度浮点型。 （4）说明函数ones（）、zeros（）、eye（）的用法。 ones（）生成全1矩阵。 zeros（）生成全0矩阵。 eye（）生成单位矩阵。 2、Matlab的数值显示格式

数字水印技术DCT算法MATLAB源代码

%Name: Chris Shoemaker %Course: E ER-280 - Digital Watermarking %Project: Block DCT Based method, using comparision between mid-band coeffcients % Watermark Embeding clear all; % save start time start_time=cputime; k=50; % set minimum coeff difference blocksize=8; % set the size of the block in cover to be used for each bit in watermark % read in the cover object file_name='_lena_std_bw.bmp'; cover_object=double(imread(file_name)); % determine size of cover image Mc=size(cover_object,1); %Height Nc=size(cover_object,2); %Width % determine maximum message size based on cover object, and blocksize max_message=Mc*Nc/(blocksize^2); % read in the message image file_name='_copyright.bmp'; message=double(imread(file_name)); Mm=size(message,1); %Height Nm=size(message,2); %Width % reshape the message to a vector message=round(reshape(message,Mm*Nm,1)./256); % check that the message isn't too large for cover if (length(message) > max_message) error('Message too large to fit in Cover Object') end % pad the message out to the maximum message size with ones message_pad=ones(1,max_message); message_pad(1:length(message))=message; % generate shell of watermarked image watermarked_image=cover_object;

数字信号处理指导书matlab版

实验1 时域离散信号的产生 一、实验目的 学会运用MATLAB 产生常用离散时间信号。 二、实验涉及的matlab 子函数 1、square 功能：产生矩形波 调用格式： x=square(t)；类似于sin （t ），产生周期为2*pi ，幅值为+—1的方波。 x=square(t ，duty)；产生制定周期的矩形波，其中duty 用于指定脉冲宽度与整个周期的比例。 2、rand 功能：产生rand 随机信号。 调用格式： x=rand （n ，m ）；用于产生一组具有n 行m 列的随机信号。 三、实验原理 在时间轴的离散点上取值的信号，称为离散时间信号。通常，离散时间信号用x （n ）表示，其幅度可以在某一范围内连续取值。 由于信号处理所用的设备主要是计算机或专用的信号处理芯片，均以有限的位数来表示信号的幅度，因此，信号的幅度也必须“量化”，即取离散值。我们把时间和幅度上均取离散值的信号称为时域离散信号或数字信号。 在MATLAB 中，时域离散信号可以通过编写程序直接生成，也可以通过对连续信号的等间隔抽样获得。 下面介绍常用的时域离散信号及其程序。 1、单位抽样序列 ? ? ?≠==000 1)(k k k δ MATLAB 源程序为

1) function [x,n] = impuls (n0,n1,n2) % Generates x(n) = delta(n-n0); n=n0 处建立一个单位抽样序列% [x,n] = impuls (n0,n1,n2) if ((n0 < n1) | (n0 > n2) | (n1 > n2)) error('arguments must satisfy n1 <= n0 <= n2') end n = [n1:n2]; x = [zeros(1,(n0-n1)), 1, zeros(1,(n2-n0))]; 将上述文件存为：impuls.m，在命令窗口输入 n0=0,n1=-10,n2=11; [x,n]=impuls (n0,n1,n2); stem(n,x,’filled’) 2)n1=-5;n2=5;n0=0; n=n1:n2; x=[n==n0]; stem(n,x,'filled','k'); axis([n1,n2,1.1*min(x),1.1*max(x)]); title('单位脉冲序列'); xlabel('时间(n)'); ylabel('幅度x(n)'); 3)n1=-5;n2=5;k=0; n=n1:n2; nt=length(n); %求n点的个数 nk=abs(k-n1)+1; %确定k在n序列中的位置 x=zeros(1,nt); %对所有样点置0 x(nk)=1; %对抽样点置1 stem(n,x,'filled','k'); axis([n1,n2,0,1.1*max(x)]); title('单位脉冲序列'); xlabel('时间(n)'); Ylabel('幅度x(n)');

MATLAB上机考试题(一)

(1)在MATLAB的命令窗口中执行_____命令，将命令窗口的显示内容清空。（） A.clear B.clc C.echo off D.cd (2)在MATLAB的命令窗口中执行_____命令，使数据输出显示为十六进制表示。（） A.format long B.format rat C.format hex D.format short e (3)下列变量名中_____是合法的。（） A.x*y,a,1 B.x\y,a1234 C.end,1 bcx D.char_1,i,j (4)已知x=0:5,则x有_____个元素。（） A.5 B.6 C.7 D.8 (5)一下运算符中哪个的优先级最高_____。（） A./ B.^ C.~= D.& (6)使用检测函数isnumeric(10)的结果是_____。（） A.1 B.0 C.false D.true (7)三维图形中默认视角是_____度。（） A.方位角=0 俯仰角=90 B.方位角=90 俯仰角=0 C.方位角=37.5 仰俯角=30 D.方位角=0 仰俯角=180 (8)将符号表达式化简为因式分解因式分解因式分解因式分解形式，使用_____函数。（） A.collect B.expand C.horner D.factor (9)运行以下命令，则_____描述是正确的。（）>>syms a b c d >>A=[a b;c d] A.A占用的内存小于100B B.创建了5个符号变量 C.A占用的内存是a b c d的总和 D.不存在 （10）已知数组a=[1 2 3;4 5 6;7 8 9],则a(:,end)是指_____元素。 （11）运行命令bitor（8,7）的结果是_____。 （12）运行以下命令： >>x=0:10; >>y1=sin(x); >>y2=5*sin(x); >>y3=[10*sin(x );20*sin(x)]; >>plot(x,y1,x,y2,x,y3) 则在一个图形窗口中，可以看到_____条曲线。 （13）符号表达式“g=sym(sin(a*z)+cos(w*v))”中的自由符号变量是_____。 （14）运行以下命令： >>syms t >>f1=1/t >>limitf1_r=limit(f1,'t','0','right'); 则函数limitf1_r趋向0的右极限为_____。 15.在MATLAB的命令窗口中执行______命令，使数值5.3显示为5.300000000000000e+000 A. format long B. format long e C. format short D. format short e 16.下列变量名中______是合法的。A．char_1,i,j B．1_1, a.1 C．x\y,a1234 D.end,1bcx 17.已知x=0:9,则x有_____个元素。 A.12 B.11 C.10 D.9 18.产生对角线上为全1其余为0的2行3列矩阵的命令是______ A. ones(2,3) B. ones(3,2) C. eye(2,3) D. eye(3,2) 19.已知数组a= [1 2 3 4 5 6 7 8 9] ，则运行a(:,1)=[]命令后______ A. a变成行向量 B. a数组为2行2列 C. a 数组为3行2列 D. a数组中没有元素3 20.按含义选出各个函数名：表示4舍5入到整数的是____，表示向最接近0取整的是____,表示向最接近-∞取整的是____,表示向最接近∞取整的是_____ A. round(x) B. fix(x) C. floor(x) D. ceil(x) 21.已知a=0:5,b=1:6,下面的运算表达式出错的为______ A. a+b B. a./b C. a’*b D. a*b 22.已知s=’显示”hello”’,则s的元素个数是______ A. 12 B. 9 C. 7 D.18

数字信号处理的MATLAB实现

昆明理工大学信息工程与自动化学院学生实验报告 （2011—2012 学年第二学期） 课程名称：数字信号处理开课实验室：信自楼111 2012 年 5 月 31 日年级、专业、班生医学号姓 名 成绩 实验项目名称数字信号处理的matlab 实现指导教师 教 师 评语教师签名： 年月日 一．实验目的 熟练掌握matlab的基本操作。 了解数字信号处理的MATLAB实现。 二．实验设备 安装有matlab的PC机一台。 三．实验内容 .1.求信号x(n)=cos(6.3Пn/3)+cos(9.7Пn/30)+cos(15.3Пn/30),0≤n≤29的幅度频谱. 2. 用冲击响应不变法设计一个Butterworth低通数字滤波器,要求参数为: Wp=0.2Пαp=1dB Ws=0.3Пαs=15dB 3.用双线性变换法设计一个Chebyshev高通IIR滤波器,要求参数为: Wp=0.6Пαp=1dB Ws=0.4586Пαs=15dB 4.用窗函数法设计一个低通FIR滤波器,要求参数为: Wp=0.2Пαp=0.3dB Ws=0.25Пαs=50dB 5.用频率抽样法设计一个带通FIR滤波器,要求参数为: W1s=0.2П W1p=0.35П W2p=0.65П W2s=0.8П αs=60dB αp=1dB 6.根据 4 点矩形序列，( n ) = [1 1 1 1] 。做 DTFT 变换，再做 4 点 DFT 变换。然后分别补零做 8 点 DFT 及 16 点 DFT。 7.调用filter解差分方程，由系统对u(n)的响应判断稳定性 8编制程序求解下列系统的单位冲激响应和阶跃响应。 y[n]+ 0.75y[n -1]+ 0.125y[n -2] = x[n]- x[n -1] 四．实验源程序 1. n=[0:1:29]; x=cos(6.3*pi*n/30)+cos(9.7*pi*n/30)+cos(15.3*pi*n/30);

Matlab 上机题及答案

1 一个三位整数各位数字的立方和等于该数本身则称该数为水仙花数。输出全部水仙花数。 for m=100:999 m1=fix(m/100); %求m的百位数字 m2=rem(fix(m/10),10); %求m的十位数字 m3=rem(m,10); %求m的个位数字 if m==m1*m1*m1+m2*m2*m2+m3*m3*m3 disp(m) end end 2.从键盘输入若干个数，当输入0时结束输入，求这些数的平均值和它们之和。 sum=0; n=0; val=input('Enter a number (end in 0):'); while (val~=0) sum=sum+val; n=n+1; val=input('Enter a number (end in 0):'); end if (n > 0) sum mean=sum/n end 3. 若一个数等于它的各个真因子之和，则称该数为完数，如6=1+2+3，所以6是完数。求[1,500]之间的全部完数。 for m=1:500 s=0; for k=1:m/2 if rem(m,k)==0 s=s+k; end end if m==s disp(m); end end 4. 从键盘上输入数字星期，在屏幕上显示对应英文星期的单词。 function week n=input('input the number:'); if isempty(n) errror('please input !!')

end if n>7|n<1 error('n between 1 and 7') end switch n case 1 disp('Monday') case 2 disp('Tuesday') case 3 disp('Wednesday') case 4 disp('Thursday') case 5 disp('Friday') case 6 disp('Saturday') case 7 disp('Sunday') end 5. 某公司销售电脑打印机的价格方案如下： （）如果顾客只买一台打印机，则一台的基本价格为$150。 （）如果顾客购买两台以上打印机，则第二台价格为$120。 （）第三台以后，每台$110。 写一段程序分别计算出购买1--10台打印机所需的钱数。打印机台数可以在程序开始处指定，或通过input命令读入。运行程序，计算出购买10台打印机的总价格。 写出程序，生成分别购买1--10台打印机所需价格的图表（使用fprintf命令输出图表，不允许手算）。 x=input('请输入购买的打印机台数：'); for m=1:x if m<=1 y(m)=150*m; elseif m<=2 y(m)=150+120*(m-1); else y(m)=150+120+110*(m-2); y(1,m)=y(m); end end y(x) plot(1:m,y,'r*--')