数字信号处理的matlab实现

数字信号处理的matlab实现

简谐振动的特性完全取决于振幅、频率、初相位角。

x1 = (0.5).^t;

x1 = 0.5* sin(2*pi*f*t+pi/4);

x1 = [(n - n0) >= 0]; %阶跃信号

x1 = [(n-n0) == 0]; %脉冲信号

y = sin(pi * x +eps)./(pi * x +eps);%sinc function，eps是matlab系统的精度，这里防止被零除

n = [-10:10];

alpha = -0.1 + 0.3 * j;

x = exp(alpha * n);%复指数信号

rx = real(x);

ix = imag(x);

mx = abs(x); %振幅

px = (180/pi) * angle(x); % 相位，并转换为度

x = rand(1,10);

x = randn(1,10); %guass series numbers

xp = [x x x x x];

x2 = fliplr(x);%左右折叠

第二章信号

%test sampling rule

dt = 0.01; %samping frequence for draw 100Hz

t = 0:dt:1;

f = 10;

x = sin(2 * pi * f * t + 0.3);

dt1 = 0.1; t1 = 0:dt1:1;%10Hz

x1 = sin(2 * pi * f * t1 + 0.3);

dt2 = 0.05; t2 = 0:dt2:1;%10Hz

x2 = sin(2 * pi * f * t2 + 0.3);

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

title('f = 10Hz, fs = 100Hz');

subplot(3,1,2), plot(t1, x1), hold on, stem(t1, x1), plot(t, x);

title('f = 10Hz, fs = 10Hz');

subplot(3,1,3), plot(t2, x2), hold on, stem(t2, x2), plot(t, x);

title('f = 10Hz, fs = 20Hz');

基本信号：

x1 = (0.5).^t; %指数信号

x1 = 0.5* sin(2*pi*f*t+pi/4); %余弦信号

x1 = [(n - n0) >= 0]; %阶跃信号

x1 = [(n-n0) == 0]; %脉冲信号

y = sin(pi * x +eps)./(pi * x +eps);%sinc function，eps是matlab系统的精度，这里防止被零除，sinc函数信号

n = [-10:10];

alpha = -0.1 + 0.3 * j;

x = exp(alpha * n);%复指数信号

rx = real(x);

ix = imag(x);

mx = abs(x); %振幅

px = (180/pi) * angle(x); % 相位，并转换为度x = rand(1,10); %随机信号

x = randn(1,10); %高斯随机序列

xp = [x x x x x]; %复制

信号运算：

时移，倒置，

尺度改变x2 = fliplr(x);%左右折叠

信号加

信号微分和积分

% test2 微分和积分

clear all; clc;

dt = 0.01; t = 0:dt:4*pi;

y = sin(t);

y1 = diff(y)/dt; % 微分，dt为采样间隔

for ii = 1:length(y) %积分，dt为采样间隔

y2(ii) = sum(y(1:ii)) * dt;

end

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

subplot(3,1,2), plot(t, [0 y1]); %微分后信号比原信号少一个元素，用零补

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

grid on;

xlabel('Time/s');

ylabel('Amplitude');

信号乘

注意：若参与信号乘的两个信号长度不一样，则必须进行转换之后才能在matlab中进行操作。

% test3 信号乘，滤波器滤波。注意振幅谱的绘制方法

clear all; clc;

dt = 0.02; df = 1/(6000 * dt);% 采样间隔为0.02s，则2min内数据个数为6000，df为信号频率分辨率n = 0:2999; %折叠频率之前（取前3000个）数据进行操作

f = n * df; %给出频率序列

sig = rand(1, length(n)); % 运用随机函数产生信号振幅谱

filt = [ones(1, 5/df), zeros(1, (length(n) - 5/df))]; % 理想滤波器幅频响应函数

subplot(3,1,1), plot(n * df, sig);

subplot(3,1,2), plot(n * df, filt, 'LineWidth', 3);

subplot(3,1,3), plot(n * df, sig .* filt);

grid on;

xlabel('Time/s');

ylabel('Amplitude');

信号奇偶性

第三章Fourier变换

一个带相位的余弦函数可以看成一个不带相位的正弦函数与一个不带相位的余弦函数的合成。

谐波函数是周期函数中最简单的函数，它描述的也是最简单的周期现象，在实际中遇到的周期现象往往比它复杂的多。但这些复杂函数均在一定近似程度上可分解为不同频率的正弦函数和余弦函数。可以将一种复杂的函数分解为一系列不同频率的正弦函数和余弦函数。

1 周期函数的Fourier变换

连续傅里叶级数

2 离散Fourier变换

离散傅里叶级数

对于一个有限长观测数据序列来说，使用Fourier级数法求得的各种周期或频率都是有限的。????=1

????=???????=，T为所取数据总的时间长度????Nyquist频率由取样间隔???来决定，为取样间隔2倍的倒数。

Ck表示了k次谐波的振幅大小。观察哪些分量的振幅大，哪些分量的振幅小，看出所给波列种含有哪种频率成分，如果我们想滤掉某种频率的波，直接去掉该频率的谐波即可。我们通常用频率作为横坐标，用振幅作

为纵坐标绘出图形来研究频率与振幅的关系，简称为振幅谱。

3 信号的Fourier分解与合成举例

% test3 DFT

clear all

N = 256; dt = 0.05; % data numbers and sampling intervel,sampling frequence is 20Hz

n = 0:N-1; t = n * dt; % 序号序列和时间序列

x1 = sin(2 * pi * t);

x2 = 0.5 * sin(2 * pi * 5 * t);

x = sin(2 * pi * t) + 0.5 * sin(2 * pi * 5 * t); %signals add

m = floor(N/2) + 1; %down for integer

a = zeros(1, m);

b = zeros(1, m);

for k = 0:m-1

for ii = 0:N-1

a(k+1) = a(k+1) + 2/N * x(ii+1) * cos(2 * pi * k * ii/N);%matlab's array index must be increase from 1 b(k+1) = b(k+1) + 2/N * x(ii+1) * sin(2 * pi * k * ii/N);

end

c(k+1) = sqrt(a(k+1).^2 + b(k+1).^2);

end

if (mod(N, 2) ~= 1) a(m) = a(m)/2;end

for ii = 0:N-1

xx(ii + 1) =a(1)/2;

for k = 1:m-1;

xx(ii + 1) = xx(ii + 1) + a(k+1) * cos(2 * pi * k * ii/N) + b(k+1) * sin(2 * pi * k * ii/N);

end

end

subplot(4,1,1), plot(t, x1); title('orignal signal One '),

xlabel('Time/s');

subplot(4,1,2), plot(t, x2); title('orignal signal Two'),

xlabel('Time/s');

%subplot(4,1,1), plot((0:N-1) * dt, xx);title('Composed

sitgnal');

subplot(4,1,3), plot(t, x); title('orignal signal'),

xlabel('Time/s');

subplot(4,1,4), plot((0:m-1)/(N * dt), c), title('Fourier

transform'),xlabel('Frequence/Hz'),

ylabel('Amplitude'); 此处的1Hz和5Hz的振幅与原来信号振幅不完全一致，是由于数据采样点较少导致的，即N较小。离散有限信号的频谱为周期谱。

方波及其Fourier分解系数和合成与原始方波信号对比

% test4 方波的DFT分解和合成

clear all; close all

N = 200; dt = 4/N;

for n = 1:N %方波信号

if (n * dt >= 2)

x(n) = 0.8;

else

x(n) = -0.8;

end

end

figure(1)

subplot(2,1,1), plot((1:N) * dt, x), hold on;

plot((1:N) * dt, zeros(1, N), 'k'), xlabel('Time/s'), title('Orignal signal');

a = zeros(1, N),

b = zeros(1, N);

nn = floor(N/2) + 1;

for k = 0:nn - 1%DFT分解方波信号，得到a和b

a(k+1) = 0;

b(k+1) = 0;

for ii = 0:N-1

a(k+1) = a(k+1) + 2/N * x(ii + 1) * cos(2 * pi * k * ii /N);

b(k+1) = b(k+1) + 2/N * x(ii + 1) * sin(2 * pi * k * ii /N);

end

c(k+1) = sqrt(a(k+1).^2 + b(k+1).^2);

end

subplot(2,1,2), plot((0:nn - 1) / (N * dt), c);

title ('DFT'), xlabel('Frequence/Hz'),

ylabel('Amplitude');

m = input('Input the max k value: ');%观察m的值不到N/2的时候的合成情况

if (m > floor(N/2) + 1)

error('Too large!');

end

if (mod(N, 2) ~= 1) a(nn) = a(nn)/2;

end

for ii = 0:N-1

xx(ii+1) = a(1)/2;

for k = 1:m

xx(ii+1) = x(ii+1) + a(k+1) * cos(2 * pi * k * ii /N) + b(k+1) * sin(2 * pi * k * ii /N);%合成信号end

end

figure(2)

plot((1:N) * dt,xx, (0:N-1) * dt, x)

hold on

plot((1:N) * dt, zeros(1, N), 'k'), xlabel('Time/s'), title('composed signal');

最基本的Fourier变换的形式。都是指复指数的形式。如果已知时间域内的等间隔数据，要求频率域内的振

???1幅使用????=?? ?????=0??????12??????????1,??=0,1,2,…,???1，称之为Fourier变换公式。那么，????= ?????=0??????2????????，??= 0,1,2,…,???1, 称之为逆Fourier变换公式。（已知频域求时间域的值）注意：在作Fourier变换时，所求的结果为复数形式。我们求其模，就得到振幅的相对值。真正的振幅大小与ck的关系为：

222 ???? =2 ????= ????=??? ????+????

相位的计算方法是：????=arctan?(?????????) ????????

% dfs的函数，特别注意WN，nk，WNnk的关系及怎么转换的，如何通过矩阵来计算

function [Xk] = dfs(xn, N)

n = [0:1:N-1];

k = [0:1:N-1];

WN = exp(-j * 2 * pi/N);

nk = n' * k;

WNnk = WN .^ nk;

Xk = xn * WNnk;

% idfs函数

function [xn] = idfs(Xk, N)

n = [0:1:N-1];

k = [0:1:N-1];

WN = exp(-j * 2 * pi/N);

nk = n' * k;

WNnk = WN .^ (-nk);

xn = (Xk * WNnk)/N;

%test 6 FT example

clear all

clf;

N = 100 ; dt = 1;

n = 0:N-1; t = n*dt;

xn = cos(2 * pi * 0.24 * t) + cos(2 * pi * 0.26 * t);

Xk = dfs(xn, N);

magXk = abs(Xk); phaXk = angle(Xk);

subplot(2,2,1), plot(t, xn), hold on, stem(t, xn); xlabel('Time/s'), title('orignal signal( N=100 )');

xx = idfs(Xk, N);

x = real(xx);

subplot(2,2,2), plot(t, x), xlabel('Time/s'),

title('Composed signal');

k = 0:length(magXk) - 1;

subplot(2,1,2), plot(k/(N*dt), magXk * 2/N);

xlabel('Frequence/Hz'), ylabel('Amplitude');

Fourier变换的性质：

线性:两个振动的Fourier变换的叠加在频域轴上与合成振动的Fourier 变换是一致的。

时移定理：将信号移位后，即推迟或提前几个时间单位，并不影响Fourier变换的振幅，但要影响其相位。若想要相位谱不发生变化，则需要添加Xk2 = Xk2 ./ exp(-2 * pi * j * k * 20 / N);

%test 6 FT时移性质

clear all; clf;

N = 128; dt =1;

n = 0:N-1; t = n * dt;

xn1 = cos(0.14 * pi * t);

subplot(3,2,1), plot(t, xn1), title('Orignal signal');

xn2 = cirshift(xn1, 20, N);

subplot(3,2,2), plot(t, xn2), title('After 20shift moves signal');

Xk1 = dfs(xn1, N)

magXk1 = abs(Xk1); phaXk1 = angle(Xk1);

k = 0:length(magXk1) - 1;

subplot(3,2,3), plot(k/(N * dt), magXk1 * 2/N),

ylabel('Amplitude'), title('Graph of amplitude');

subplot(3,2,4), plot(k/(N * dt), unwrap(phaXk1)),

ylabel('Phase/rad'), title('Phase graph');

% unwrap(p)用于展开弧度相位角p，即将不在-π~π的相位角归算到该区间内便于比较

Xk2 = dfs(xn2, N);

%Xk2 = Xk2./exp(- 2 * j * pi * k * 20 / N);

magXk2 = abs(Xk2); phaXk2 = angle(Xk2);

k = 0:length(magXk2) - 1;

subplot(3,2,5), plot(k/(N * dt), magXk2 * 2 /N),

ylabel('Amplitude'), title('Graph of amplitude') subplot(3,2,6), plot(k/(N * dt), unwrap(phaXk2)),

ylabel('Phase/rad'), title('Phase graph');

频移定理:如果x(n)的Fourier变换为X(k)，将X(k)移位i，即X(k-i)的Fourier逆变换为x(n)ej2πi/N。偶函数和奇函数与Fourier变换后实部和虚部的关系：

偶函数的Fourier变换对应于频率域的实部，奇函数的Fourier变换对应于频率域的虚部。即信号可分为偶函数和奇函数，其中信号偶函数部分Fourier变换的虚部以及奇函数部分Fourier变换的实部均几乎为零。

1褶积定理：两个函数h(n)和x(n)在时间域内的褶积等于在频率域内乘积。即?????=0?? ?? ?(???)的Fourier变换为

X(k)H(k)。

快速Fourier变换

用Matlab进行谱分析时应注意：

(1) 函数FFT返回值的数据结构具有对称性

归一化频率是数字信号处理中的一个特定概念。由于通常数字信号处理技术（Fourier变换，滤波器等）不是对应于特定采样间隔的信号，而是

对于任何采样间隔的数据均可以处理，因此在分析数字信号处理技术时，通常定义采样频率的一半（即Nyquist频率）为1，而将小于Nyquist频率的频率归一化，这样的频率称为归一化频率。由采样定理可知，如果某种信号的归一化频率大于1，则会出现频率折叠或频率重复。

一般而言，对于N个点的xn序列的FFT为N个点的复数序列，其点n=N/2+1对应于Nyquist频率，做谱分析时仅取序列Xk的一半，即N/2点即可。Xk的后一半与前一半是对称的。

(2) 做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在快速Fourier逆变换时已经做

了处理。要得到真实的振幅值大小，只要将得到的变换后结果乘以2除以N即可。

由于在时间域内信号加零，致使振幅谱中出现很多其他成分。其振幅由于加了多个零而明显减小。则，在对信号做频谱分析时，数据样本应有足够的长度，一般应取的FFT程序中所用数据点数与原含有信号数据点数相同，这样的频谱图具有较高的质量，可减少因补零或截断而产生的影响。

%test 6 FFT

clear all

clf;

fs = 1000; N = 1024;

n = 0:N-1; t = n/fs;

f1 = 50; f2 = 120;

x = sin(2 * pi * f1 * t ) + 2 * sin(2 * pi * f2 * t);

x = x + 1 * randn(1, length(t));

y = fft(x, N);

mag = abs(y);

f = n * fs /N; %真实频率序列

plot(f(1:N/2), mag(1:N/2) * 2/N), xlabel('Frequence/Hz'), ylabel('Amplitude'), title('Signal with noise');

%test7 FFT comparable

clear all;

fs = 100; N = 128;

n = 0:N - 1; t = n/fs;

x = sin(2 * pi * 40 * t) + sin(2 * pi * 15 * t);

subplot(2,2,1), plot(t, x), title('orignal signal');

grid on;

y = fft(x, N);

mag = abs(y);

f = n * fs / N ;

subplot(2,2,2), plot(f(1:N/2), mag(1:N/2) * 2/N), title('FFT');

xifft = ifft(y);

realx = real(xifft);

ti = [0:length(xifft) - 1]/fs;

subplot(2,2,3), plot(ti, realx);

yif = fft(xifft, N);

mag = abs(yif);

f = [0:length(y) - 1]'* fs/length(y);

subplot(2,2,4), plot(f(1:N/2), mag(1:N/2) * 2/N);

如果采样数据过少，运用Fourier变换不能分辨出其中的频率成分。添加零后可增加频谱中的数据个数，谱的密度增高，但仍不能分辨其中的频率成分，即谱的分辨率没有提高。只有数据点数足够多时才能分辨出其中的频率成分。

%test 8 利用Fourier变换进行滤波

clear all, clf;

dt = 0.02; N = 512;

n = 0:N-1; t = n * dt; f = n/(N * dt);

f1 = 3; f2 = 10;

x = .5 * sin(2 * pi * f1 * t) + cos(2 * pi * f2 * t);

subplot(2,2,1), plot(t, x);

y = fft(x);

subplot(2,2,2), plot(f, abs(y) * 2/N);

xlim([0 50]);

f1 = 8; f2 = 15;

yy = zeros(1, length(y));

for m = 0:N - 1

if (m/(N * dt) > f1 & m/(N * dt) < f2) | (m/(N * dt) > (1/dt - f2) & m/(N * dt) < (1/dt - f1)) % 置在此频率范围内的振动幅度为零

%小于Nyquist频率的滤波范围，大于Nyquist频率的滤波范围

%1/dt为一个频率周期

yy(m+1) = 0.;

else

yy(m+1) = y(m+1);%其余频率范围的振动幅度不变

end

end

subplot(2,2,4), plot(f, abs(yy) * 2/N), xlim([0 50]); %滤波后的振幅图

gstext = sprintf('%4.1f - %4.1fHz的频率被过滤',f1, f2);

title(gstext);

subplot(2,2,3), plot(t, real(ifft(yy)));

title('通过IFFT回到时间域');

第四章系统

线性连续时间系统，离散时间系统

传递函数与输入和输出信号的特征无关，它描述了系统的一种特性。如果系统是一个放大器，则对应于放大器的放大倍数，跟输入信号的特征无关。

任何离散系统可以由z传递函数来描述。系统的频率响应是一种特殊的传递函数。

1卷积的matlab实现y=conv(x,h)，其中（两个序列都是从n=0开始的）

2 滤波方式的matlab实现y=filter(b, a, x)

%test 绘制幅频响应和相频响应

clear all, clf;

b = [1 2 1]; a = [1 0.5 0.25];

m = 0:length(b) - 1;

n = 0:length(a) - 1;

K = 512;

k = 1:K; w = pi * k/K;

num = b * exp(-j * m' * w);

den = a * exp(-j * n' * w);

H = num./den;

magH = 20 * log10(abs(H));

phaH = angle(H);

subplot(2,1,1), plot(w/pi, magH);

title('幅频图');

grid on;

subplot(2,1,2), plot(w/pi, phaH * 180/pi), title('相频图'); %test 绘出系统的输出信号波形

clear all;

b = 0.15; a = [1 -0.8];

N = 100;dt = 1; n = 0:N-1; t = n * dt;

x = 2 * sin(0.05 * pi * t) + 0.2 * randn(1, N);

imp = [1,zeros(1, N-1)];

h = filter(b, a, imp);

yc = conv(h, x);

y = yc(1:N);

y1 = filter(b, a, x);

subplot(3,1,1), plot(t, x), title('input signal');

subplot(3,1,2), plot(t, y), title('after conv');

subplot(3,1,3), plot(t, y1), title('after filter');

xlabel('Time/s');

运用卷积计算结果和滤波器输出结果的比较可知，这两种方法得到的结果完全相同。

第五章模拟滤波器设计

若滤波器的输入、输出都是离散时间信号，那么该滤波器的脉冲响应必然是离散的。我们称这样的滤波器为数字滤波器。当用硬件实现一个数字滤波器时，所需的元件是延迟器、乘法器和加法器。当在计算机上用软件实现时，它就是一段线性卷积的程序。模拟滤波器只能用硬件来实现，其元件是电阻、电容、电感及运算放大器。

模拟滤波器的设计原理

1 信号无失真的条件

系统的频率响应H(jω)满足下面的特性：

?? ???? =?? ?? ?? =∠?? ???? = ???????

滤波器应该具有无限宽的定值幅频域与线性相频。通常定义群延迟为信号系统的延迟时间。即

????(??)????=?群延迟为相频特性曲线的斜率。对于信号无失真传输，

td为常数，即群延迟为常数；否则它是频率的非线性函数。

常用的模拟原型滤波器包括：Butterworth、Chebyshev I、Chebychev Ⅱ、Elliptical、Bessel原型低通滤波器。模拟原型滤波器指的是截止频率为1的滤波器。后面所介绍的各类模拟滤波器和数字滤波器可通过这些低通原型滤波器变换得到。

1）Butterworth滤波器的特点：通带内具有最大平坦的频率特性，且随着频率增大平滑单调下降；阶数愈高，特性愈接近矩形，过渡带愈窄，传递函数无零点。

Matlab信号处理工具箱提供Butterworth模拟低通滤波器原型设计函数buttap，函数调用为：

[z, p , k] = buttap(n);其中，n为butterworth滤波器阶数；zpk分别为滤波器的零点，极点和增益。

[b, a] = zp2tf(z, p, k)为将模拟原型滤波器函数设计出的零点z，极点p，增益k形式转换为传递函数形式；其中，b为滤波器传递函数分子多项式系数，a为滤波器传递函数分母多项式系数，

[H, ω] = freqs(b,a, ω);求出传递函数形式（分子和分母多项式的系数为b，a）表示的滤波器的对应于频率点ω的复数频率响应H（包括实部和虚部），这里ω为一个矢量，表示对应的角频率。若该函数不写出输出变量，则执行后绘出该滤波器的幅频响应和相频响应图。

%test Butterworth

n = 0:0.01:2; %频率点

for ii = 1:4 %取四种滤波器

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上机实验(答案)

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数字信号处理

MATLAB 下的数字信号处理实现示例 附录一信号、系统和系统响应 1、理想采样信号序列 （1）首先产生信号x(n),0<=n<=50 n=0:50; A=444.128; a=50*sqrt(2.0)*pi; T=0.001; w0=50*sqrt(2.0)*pi; x=A*exp(-a*n*T).*sin(w0*n*T); close all subplot(3,1,1);stem(x); %定义序列的长度是50 %设置信号有关的参数 %采样率 %pi 是MATLAB 定义的π，信号乘可采用“.*”%清除已经绘制的x(n)图形 %绘制x(n)的图形 title(…理想采样信号序列?); （2）绘制信号x(n)的幅度谱和相位谱 k=-25:25; W=(pi/12.5)*k; X=x*(exp(-j*pi/12.5)).^(n?*k); magX=abs(X); %绘制x(n)的幅度谱subplot(3,1,2);stem(magX);title(…理想采样信号序列的幅度谱?); angX=angle(X); %绘制x(n)的相位谱subplot(3,1,3);stem(angX) ; title (…理想采样信号序列的相位谱?) （3）改变参数为：A = 1,? = 0.4, & 0 = 2.0734, T = 1 n=0:50; A=1; a=0.4; T=1; w0=2.0734; x=A*exp(-a*n*T).*sin(w0*n*T); close all subplot(3,1,1);stem(x); title(…理想采样信号序列?); k=-25:25; W=(pi/12.5)*k; X=x*(exp(-j*pi/12.5)).^(n?*k); magX=abs(X); %定义序列的长度是50 %设置信号有关的参数 %采样率 %pi 是MATLAB 定义的π，信号乘可采用“.*” %清除已经绘制的x(n)图形 %绘制x(n)的图形 %绘制x(n)的幅度谱 subplot(3,1,2);stem(magX);title(…理想采样信号序列的幅度谱?); angX=angle(X); %绘制x(n)的相位谱

南京理工大学数字信号处理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试卷

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的数值显示格式

数字信号处理指导书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*--')

数字信号处理MATLAB实验1

实验一熟悉MATLAB环境 一、实验目的 (1)熟悉MATLAB的主要操作命令。 (2)学会简单的矩阵输入和数据读写。 (3)掌握简单的绘图命令。 (4)用MATLAB编程并学会创建函数。 (5)观察离散系统的频率响应。 二、实验内容 认真阅读本章附录，在MATLAB环境下重新做一遍附录中的例子，体会各条命令的含义。在熟悉了MATLAB基本命令的基础上，完成以下实验。 上机实验内容： (1)数组的加、减、乘、除和乘方运算。输入A=[1234]，B=[345 6]，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.^B并用stem语句画出 A、B、C、D、E、F、G。 (2)用MATLAB实现以下序列。 a）x(n)=0.8n0≤n≤15 b）x(n)=e(0.2+3j)n0≤n≤15 c）x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π)0≤n≤15 (n)=x(n+16),绘出四个d）将c)中的x(n)扩展为以16为周期的函数x 16 周期。 (n)=x(n+10),绘出四个e）将c)中的x(n)扩展为以10为周期的函数x 10 周期。

(3)x(n)=[1,-1,3,5]，产生并绘出下列序列的样本。 a）x 1(n)=2x(n+2)-x(n-1)-2x(n) b）∑=-=5 1k 2) k n (nx (n) x (4)绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注。 a)x(t)=sin(2πt)0≤t≤10s b)x(t)=cos(100πt)sin(πt) 0≤t≤4s (5)编写函数stepshift(n0,n1,n2)实现u(n-n0),n1

MATLAB实现数字信号处理

《数字信号处理》课程设计实例： 声音信号的处理 一．摘要： 这次课程设计的主要目的是综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB或者ＤＳＰ开发系统作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。通过对声音的采样，将声音采样后的频谱与滤波。 MATLAB全称是Matrix Laboratory，是一种功能强大、效率高、交互性好的数值和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。。经过多年的发展，已经发展成为一种功能全面的软件，几乎可以解决科学计算中所有问题。MATLAB软件还提供了非常广泛和灵活的用于处理数据集的数组运算功能。 在本次课程设计中，主要通过MATLAB来编程对语音信号处理与滤波，设计滤波器来处理数字信号并对其进行分析。 二．课程设计目的： 综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。 三．设计容： 容：录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性

变换法设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；换一个与你性别相异的人录制同样一段语音容，分析两段容相同的语音信号频谱之间有什么特点；再录制一段同样长时间的背景噪声叠加到你的语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除。 四．设计原理： 4.1．语音信号的采集 熟悉并掌握MATLAB中有关声音（wave）录制、播放、存储和读取的函数，在MATLAB环境中，有关声音的函数有： a：y=wavrecord(N,fs,Dtype);利用系统音频输入设备录音，以fs为采样频率，默认值为11025，即以11025HZ进行采样。Dtype为采样数据的存储格式，用字符串指定，可以是：‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度进行采样，其它三种都是16位采样结果转换为指定的MATLAB数据； b：wavplay(y,fs)；利用系统音频输出设备播放，以fs为播放频率，播放语音信号y； c：wavwrite((y,fs,wavfile);创建音频文件； d：y=wavread(file);读取音频文件； 关于声音的函数还有sound()；soundsc()；等。 4.2滤波器： 4．21.IIR滤波器原理 冲激响应不变法是使数字滤波器在时域上模拟滤波器，但是它们的缺点是产生频率响应的混叠失真，这是由于从s平面到z平面是多值的映射关系所造成的。 双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了克服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里，再通过变换关系将此横带变换到整个z平面上去，这样就使得s平面与z平面是一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象。

数字信号处理基本知识点Matlab实现

数字信号处理（第二版） 绪论 1.4 MATLAB 在信号处理中的应用简介 MATLAB 是美国Mathworks 公司于1984年推出的一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算、信号处理、系统仿真和图形显示于一体，从而被广泛地应用于科学计算、控制系统、信息处理等领域的分析、仿真和设计工作。 MATLAB 软件包括五大通用功能：数值计算功能（Numeric ），符号运算功能（Symbolic ）；数据可视化功能（Graphic ），数据图形文字统一处理功能（Notebook ）和建模仿真可视化功能（Simulink ）。该软件有三大特点：一是功能强大；二是界面友善、语言自然；三是开放性强。目前，Mathworks 公司已推出30多个应用工具箱。MA TLAB 在线性代数、矩阵分析、数值及优化、数理统计和随机信号分析、电路与系统、系统动力学、信号和图像处理、控制理论分析和系统设计、过程控制、建模和仿真、通信系统、以及财政金融等众多领域的理论研究和工程设计中得到了广泛应用。 2.10 离散时间信号与系统的Matlab 表示 2.10.1 离散时间信号的表示和运算 1、基本序列的Matlab 表示 单位采样序列 在MA TLAB 中，单位采样序列可以通过编写以下的DTimpulse .m 文件来实现，即 function y=DTimpulse (n) y=(n==0); %当参数为0时冲激为1,否则为0 调用该函数时n 必须为整数或整数向量。 单位阶跃序列 在MA TLAB 中，单位阶跃序列可以通过编写DTu .m 文件来实现，即 function y=DTu (n) y=n>=0; %当参数为非负时输出1 调用该函数时n 必须为整数或整数向量。 矩形序列 用MA TLAB 表示矩形序列可根据公式()()()N R n u n u n N =--并利用DTu 函数生成，即 function y=DTR(n,N) y=DTu(n)-DTu(n-N); 调用该函数时n 必须为整数或整数向量，N 必须为整数。 实指数序列 用MA TLAB 表示实指数序列()(),n x n a u n n N a R =∈∈，即

基于MATLAB的数字信号处理

数字信号处理课程设计报告题目：语音数字信号处理与分析及 Matlab实现 系别通信工程 专业班级 学生姓名 学号 指导教师 提交日期

摘要 本次课程设计综合利用数字信号处理的理论知识进行语音信号的频谱分析，通过理论推导得出相应结论，再利用MATLAB作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。本次课程设计要求利用MATLAB对语音信号进行分析和处理，要求学生采集语音信号后，在MATLAB软件平台进行频谱分析；并对所采集的语音信号加入干扰噪声，对加入噪声的信号进行频谱分析，设计合适的滤波器滤除噪声，恢复原信号。待处理语音信号是一个在20Hz～20kHz 频段的低频信号。采用了高效快捷的开发工具——MATLAB，实现了语音信号的采集，对语音信号加噪声及设计滤波器滤除噪声的一系列工作。利用采样原理设计了高通滤波器、低通滤波器、带通滤波器、带阻滤波器。同学通过查阅资料自己获得程序进行滤波器的设计，能过得到很好的锻炼。 关键词：MATLAB滤波器数字信号处理

目录 第一章绪论 (1) 1.1设计的目的及意义 (1) 1.2设计要求 (1) 1.3设计内容 (1) 第二章系统方案论证 (3) 2.1设计方案分析 (3) 2.2实验原理 (3) 第三章信号频谱分析 (6) 3.1原始信号及频谱分析 (6) 3.2加入干扰噪声后的信号及频谱分析 (7) 第四章数字滤波器的设计与实现 (11) 4.1高通滤波器的设计 (11) 4.2低通滤波器的设计 (12) 4.3带通滤波器的设计 (15) 4.4带阻滤波器的设计 (16) 第五章课程设计总结 (19) 参考文献 (20) 附录Ⅰ..................................................................................I 附录Ⅱ................................................................................II