DSP using MATLAB中的函数学习

为了解决matlab自带的函数conv中没有时间的信息，所以改写了conv_m以达到保留序列的时间（顺序）的目的

function [y,ny] = conv_m(x,nx,h,nh)

% Modified convolution routine for signal processing
% -------------------------------------------------
% [y,ny] = conv_m(x,nx,h,nh)
%  y = convolution result
% ny = support of y
%  x = first signal on support nx
% nx = support of x
%  h = second signal on support nh
% nh = support of h
%
nyb = nx(1)+nh(1); nye = nx(length(x)) + nh(length(h));


ny = [nyb:nye];


y = conv(x,h);

这段代码中我认为比较经典的是在

nyb = nx(1)+nh(1); nye = nx( length(x) + nh(length(h)) );

这一句中通过length求得长度，这样就相当于在没有序列输入的时候将值置位0，在最后不显示

或者说是允许在输入序列的时候序列的标号长度与序列值的长度不同。

这一段代码是傅里叶系数的获取

function [Xk] = dfs(xn,N)

% Computes Discrete Fourier Series Coefficients
% ---------------------------------------------
% [Xk] = dfs(xn,N)
% Xk = DFS coeff(系数的意思啦). array over 0 <= k <= N-1
% xn = One period of periodic signal over 0 <= n <= N-1
%  N = Fundamental period of xn
%


n = [0:1:N-1];                       % row vector for n


k = [0:1:N-1];                       % row vecor for k


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


nk = n'*k;                           % creates a N by N matrix of nk values


WNnk = WN .^ nk;                     % DFS matrix


Xk = xn * WNnk;                      % row vector for DFS coefficients

这里就摘抄一下书上的原话吧

”进行傅里叶级数“系数，为了计算每一个样本X~(k)，能够作为一个for...end 循环的和来实现；为了计算所有的DFS系数就需要另一个for...end循环，这就会得到一种嵌套两个for...end循环的实现。在MATLAB中，这显然效率很低的。MATLAB的一种高效实现是采用矩阵向量乘法；在实现离散时间傅丽叶变换的数值近似中，早先已用过这一办法。

这里指的就是上面程序中的WNnk的获得啦。

还有就是MATLAB毕竟是矩阵实验室，所以对矩阵还是比较感兴趣，对循环，那是C++的事情啦。

下面是对系统的差分方程输出的求解过程。主要是通过LTI的线性时不变性拆分成零输入和零状态。

clear all;close all; 

num=[2 -1 0]; %系统函数分子的系数
den=[2 -1 -3]; %系统函数分母的系数
n=0:50; 
nl=length(n); 
%求零输入响应
y01=[1 3]; %y的初始状态
x01=[0 0]; %x 的初始状态
x1=zeros(1,nl); 
zi1=filtic(num,den,y01,x01); 
%为filter函数准备初始值
y1=filter(num,den,x1,zi1); 

 
%求零输入响应
subplot(311); 
stem(n,y1,'r.'); 
title('零输入响应'); 
grid on; 
%求零状态响应
y02=[0 0]; 
x02=[0 0]; 
x2=0.5.^n; 
zi2=filtic(num,den,y02,x02); 
y2=filter(num,den,x2,zi2); 
subplot(312); 
stem(n,y2,'r.');  

上面是百度文库的一段例子，感觉可以很好的完成差分方程的求解啦。

这里主要注意的是filtic的运用。因为以前看到的都是直接用filter

它山之石