序列运算的MATLAB实现
序列之和
function [y,n] = seqadd(x1,n1,x2,n2)
% 序列之和
n = min(min(n1),min(n2)):max(max(n1),max(n2));
y1 = zeros(1,length(n));y2 = y1;
y1((n >= min(n1))&(n <= max(n1)) == 1) = x1;
y2((n >= min(n2))&(n <= max(n2)) == 1) = x2;
y=y1+y2;
输入:
[y, n] = seqadd([1,2,3],[0,1,2],[1,2,3],[1,2,3])
输出:
y = [1, 3, 5, 3]
n = [0, 1, 2, 3]
序列之乘积
function [y,n] = seqmult(x1,n1,x2,n2)
% 序列之积
n = min(min(n1),min(n2)):max(max(n1),max(n2));
y1 = zeros(1,length(n));y2 = y1;
y1((n >= min(n1))&(n <= max(n1)) == 1) = x1;
y2((n >= min(n2))&(n <= max(n2)) == 1) = x2;
y=y1.*y2;
输入:
[y, n] = seqmult([1,2,3],[0,1,2],[1,2,3],[1,2,3])
输出:
y = [0, 2, 6, 0]
n = [0, 1, 2, 3]
综合运用例子:
% 已知两序列为x1,x2有
% x1 = [1,3,5,7,6,4,2,1],起始位置ns1 = -3
% x2 = [4,0,2,1,-1,3],起始位置ns2 = 1
% 求它们的和ya以及乘积ymx1 = [1,3,5,7,6,4,2,1];ns1 = -3;
x2 = [4,0,2,1,-1,3];ns2 = 1;
nf1 = ns1 + length(x1)-1;
nf2 = ns2 + length(x2)-1;
n1 = ns1:nf1;
n2 = ns2:nf2;
[ya,t1] = seqadd(x1,n1,x2,n2);
[ym,t2] = seqmult(x1,n1,x2,n2);figure(1)
subplot(2,2,1);stem(n1,x1,'.');ylabel('x1(n)');grid on;
subplot(2,2,3);stem(n2,x2,'.');ylabel('x2(n)');grid on;
subplot(2,2,2);stem(t1,ya,'.');ylabel('x1(n)+x2(n)');grid on;
subplot(2,2,4);stem(t2,ym,'.');ylabel('x1(n)*x2(n)');grid on;
输出:
乘常数
function y = mult_a(x,a)
y = a*x;
输入:
y = mult_a([1,2,3],3)
输出:
y = [3, 6, 9]
移位
function [y,ny] = seqshift(x,nx,m)
% 实现y(n) = x(n-m)
% -----------------
% [y,n] = seqshift(x,nx,m)
%
ny = nx + m;y = x;
输入:
[y,ny] = seqshift([1,2,3],[0,1,2],2)
输出:
y = [1, 2, 3]
ny = [2, 3, 4]
翻褶
function [y,ny] = seqfold(x,nx)
% 序列翻褶
y = fliplr(x);ny = -fliplr(nx);
输入:
[y,ny] = seqfold([1,2,3],[0,1,2])
输出:
y = [3, 2, 1]
ny = [-2, -1, 0]
序列的能量
function ex = seqenergy(x)
% ex = sum(x.*conj(x));
ex = sum(abs(x).^2);
输入:
ex = seqenergy([1,2,3])
输出:
ex = 14
序列的功率
只需将序列在一个周期N内的能量除以N。
序列的卷积
不能提供位置信息的序列卷积
使用MATLAB内置的conv函数。
输入:
y = conv([1,2,3],[1,2,3])
输出:
y = [1, 4, 10, 12, 9]
可以提供位置信息的序列卷积
function[y,ny] = convwthn(x,nx,h,nh)
% 序列的卷积
ny1 = nx(1) + nh(1);
ny2 = nx(end) + nh(end);
y = conv(x,h);ny = ny1:ny2;
输入:
[y,ny] = convwthn([1,1,1,1],0:3,[1,2,3,4],-1:2)
输出:
y = [1, 3, 6, 10, 9, 7, 4]
ny = [-1, 0, 1, 2, 3, 4, 5]
综合运用的例子:
% 用convwthn求解x(n) = [1, 2, 3, -1, -2],nx = [-1, 3]
% 与h(n) = [2, 2, 1, -1, 4, -2],nh = [-3, 2]的卷积
clc;clear;x = [1, 2, 3, -1, -2];nx = -1:3;
h = [2, 2, 1, -1, 4, -2];nh = -3:2;
[y,ny] = convwthn(x,nx,h,nh);
stem(ny,y,'.');xlabel('n');ylabel('y(n)');grid on;
输出:
序列的相关
不能提供位置信息的序列相关
输入:
rxy = xcorr([1,2,3],[1,2,1])
输出:
rxy = [1, 4, 8, 8, 3]
可以提供位置信息的序列相关
% x(n) = [_2_,1,3,2,1,5,1];y(n) = [2,1,_3_,4];
% 采用x(m)*y(-m)公式,即用卷积代替相关
clc;clear;x = [2,1,3,2,1,5,1];y = [2,1,3,4];
nx = 0:length(x) - 1;ny = -2:1;
[yf,nyf] = seqfold(y,ny);
rxy = convwthn(x,nx,yf,nyf);k = nx(1) + nyf(1):nx(end) + nyf(end);
stem(k,rxy,'.');
xlabel('n');ylabel('rxy(n)');
axis([k(1),k(end),0,35]);grid on;
输出:
说明
文章中出现的代码参考:数字信号处理教程/程佩青. 第四版: 清华大学出版社。本人在其代码上有一定改进。本人写这一系列文章纯粹出于自己及广大学子学习MATLAB与数字信号处理的目的。欢迎指正!