数字信号处理学习笔记（MATLAB)

1、求解线性常系数差分方程
B、A分别是差分方程的系数向量
xn是输入信号
xi是等效初始条件的输入序列，由初始条件确定

计算零状态响应：yn=filter(B,A,xn)
计算完全相应：yn=filter(B,A,xn,xi)
xi是等效初始条件的输入序列，由初始条件确定：xi=filtic(B,A,ys,xs)
ys和xs是初始条件向量，如果xn是因果序列，则xs=0

a=0.8;ys=1;%设置差分方程系数a=0.8;初始值ys=y(-1)=1
xn=[1,zeros(1,30)];%输入序列x(n)为单位脉冲序列，长度N=31
B=1;A=[1,-a];%差分方程系数
xi=filtic(B,A,ys);%计算等效初始条件的输入序列
yn=filter(B,A,xn,xi);
n=0:length(yn)-1;
stem(n,yn,'o','LineWidth',2);
axis([0,30,0,2]);
title('y(-1)=1时的响应y(n)');
xlabel('n');
ylabel('y(n)');

2、信号采样和重建
重构目的：离散时间信号->连续时间信号
重构方式：频率滤波、时域内插

% Analog signal
Dt = 0.00005;
t = - 0.005:Dt:0.005;
xa = exp(-1000 * abs(t));
 
subplot(3,1,1);
plot(1000*t,xa);
title('Analog signal');
xlabel('t in msec');
ylabel('xa');
 
% Discrete-time signal
Ts = 0.0002;
Fs = 1/Ts;
n = -25:25;
nTs = n*Ts;
x = exp(-1000*abs(nTs));
subplot(3,1,2)
stem(n*Ts*1000,x);
title('Discrete-time signal');

% Analog signal reconstruction
xa_r = x * sinc( Fs * ( ones(length(n),1) * t - nTs'* ones(1,length(t) ) ));
%实际完成了求和运算

subplot(3,1,3);
plot(1000*t,xa_r);
title('Analog signal reconstruction');
xlabel('t in msec');
ylabel('xa after reconstruction');
hold on 
stem(n*Ts*1000,x)
hold off
%check
error = max(abs(xa_r - xa ))

3、傅里叶变换三维视角

%   傅里叶变换的理解
clear all
close all
x = -6:0.01:6;  %时间
tin=1:60;       %频域的范围
count = 30;      %显示谐波的数量为2*count+1
y = zeros(count, length(x));
figure(1)
for N = 1:count
    k = 2*N-1;
    HarmonicSig = 0.5 + sinc(k/2)*cos(k*pi/2*x);    %方波的傅里叶级数展开
    yin=k*ones(1,length(x));    %奇次谐波的y坐标
    yin1=2*N*ones(1,length(x)); %偶次谐波的y坐标
    %y为所有谐波的和，结果就是原始的方波信号
    if N==1
        y(N,:) = 0.5 + sinc(k/2)*cos(k*pi/2*x);
    else
        y(N,:) = y(N-1,:) + sinc(k/2)*cos(k*pi/2*x);
    end
    zin(N,:)=HarmonicSig;
    plot3(yin1,x,zeros(1,length(x)));%偶次谐波都为0
    hold on
    plot3(yin,x,zin(N,:));
    %下面是每次谐波的傅里叶变换。所有谐波的和即为方波的频域表示
    if N==1
        spec=abs(fft(zin(N,:)-mean(zin(N,:))));
    else
        spec=spec+abs(fft(zin(N,:)-mean(zin(N,:))));
    end
end
yin1=(25)*ones(1,length(x));
plot3(yin1,x,y(N,:));   %方波的傅里叶级数展开的和
stem3((tin/6+0.333)*2-1,7+zeros(1,60),spec(1:60)/max(spec(1:60)))
% axis([0 10 0 1])
hold off
axis tight
x1=xlabel('频率方向——>');
x2=ylabel('时间方向——>');
% set(x1,'Rotation',15); %x轴名称旋转
% set(x2,'Rotation',-15); %y轴名称旋转
title('傅里叶的理解')

%%  下面是傅里叶级数展开的理解
f = x;
T = 4;
for N = 1:length(f)
    temp = rem(abs(x(N)),T);
    if temp>1 && temp<3
        f(N) = 0;
    else
        f(N) = 1;
    end
end
figure(2)
N=count;
plot(x,f,'k');
hold on
h = plot(x,y(N,:),'r');
hold off
% figure(3)
% plot(tin/6-0.667,spec/max(spec))
% axis([0 10 0 1])

4、因式分解

syms z
f=z^2+z-6;
fnew=factor(f)

5、零极点图、频率响应
zplane(z,p)：绘制列向量z的零点和列向量的极点，同时画出参考单位圆，并在多阶零点和极点的右上角标出其阶数。如果z和p为矩阵，则zplane以不同的颜色分别绘出z和p各列中的零点和极点。
[H,w]=freq(B,A,N)：B和A分别为离散系统的系统函数的分子、分母多项式的系数向量，N为正整数，返回向量H则包含了离散系统频响在0—pi范围内N个频率等分点。调用中若N默认，默认值为512。