matlab信号处理相关函数程序

序列的傅里叶变换及其逆变换定义：

)(|)(|)()(ω?ωωω

j j n n j j e e X e

n x e X ==∑∞-∞=- 1(n)(e )e 2jw jwn x X dw π

ππ-=?

其幅度特性为|)(|ωj e X ，在Matlab 中采用abs 函数；相位特性为)(ω?，在Matlab 中采用angle 函数。

为了方便，考虑在两个周期，例如[ππ2,2-]中2M+1个均匀频率点上计算FT ，并且观察其周期性和对称性。为此给出function 文件如下，求解FT 变换：

function [X,w]=ft(x,n,k)

%X ：序列x(n)的傅里叶变换

%w ：X 的自变量

%x ：要进行傅里叶变换的序列x(n)

%n ：序列x(n)的位置向量

%k ：求和区间

w=(pi/abs(max(k)/2))*k;

X=x*(exp(-1i*pi/abs(max(k)/2))).^(n'*k);

使用方法如下：

n=-5:5;%序列区间

x=(-0.9).^n;%序列表达式

k=-200:200;%求和区间

[Xw,w]=ft(x,n,k);%求傅里叶变换

magX=abs(Xw);%求幅度

angX=angle(Xw);%求相位

realX=real(Xw);

imagX=imag(Xw);

subplot(2,2,1)

plot(w/pi,magX)%绘制幅度曲线

grid on

title('幅度曲线')

xlabel('\omega/\pi')

ylabel('幅度')

xmin=0;

xmax=2;

set(gca,'xlim',[xmin,xmax],'ylimmode','auto','zlimmode','auto'); %xmin xmax 为范围 subplot(2,2,2)

plot(w/pi,angX/pi) %绘制相位曲线

grid on

title('相位曲线')

xlabel('\omega/\pi')

ylabel('相位')% angle(X)/pi

xmin=0;

xmax=2;

set(gca,'xlim',[xmin,xmax],'ylimmode','auto','zlimmode','auto'); %xmin xmax 为范围 subplot(2,2,3)

plot(w/pi,realX) %绘制实部曲线

grid on

title('实部曲线')

xlabel('\omega/\pi')

ylabel('实部')

xmin=0;

xmax=2;

set(gca,'xlim',[xmin,xmax],'ylimmode','auto','zlimmode','auto'); %xmin xmax 为范围 subplot(2,2,4)

plot(w/pi,imagX) %绘制虚部曲线

grid on

title('虚部曲线')

xlabel('\omega/\pi')

ylabel('虚部')

xmin=0;

xmax=2;

set(gca,'xlim',[xmin,xmax],'ylimmode','auto','zlimmode','auto'); %xmin xmax 为范围

序列的DFT 及IDFT 定义：

∑∑-=-=-==10102)()()(N n N n nk N j nk N e

n x W n x k X π

21100

11(n)(k)W (k)e N N j nk nk N N k k x X X N N π---====∑∑ 离散傅里叶变换的的性质：

（1）DFT 的共轭对称性。若)()()(n x n x n x op ep +=，[])()(n x DFT k X =，则：)()]([k X n x DFT R ep =， )()]([k jX n x DFT I op =。

（2）实序列DFT 的性质。若)(n x 为实序列，则其离散傅里叶变换)(k X 为共轭对称，即10),()(*-≤≤-=N k k N X k X 。

（3）实偶序列DFT 的性质。若)(n x 为实偶序列，则其离散傅里叶变换)(k X 为实偶对称，即10),()(-≤≤-=N k k N X k X 。

（4）实奇序列DFT 的性质。若)(n x 为实奇序列，则其离散傅里叶变换)(k X 为纯虚奇对称，即10),()(-≤≤--=N k k N X k X 。

离散傅立叶变换函数

function [Xk,k]=dft(xn,N)

n=0:1:N-1;

k=0: N-1;

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

nk=n'*k;

WNnk=WN.^(nk);

Xk=xn*WNnk; %采用矩阵相乘的方法

magX=abs(Xk);

k=(0:length(magX)-1)*N/length(magX);

离散傅立叶反变换函数

function [xn]=idft(Xk,N)

n=0:1:N-1;

k=0:1:N-1;

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

nk=n'*k;

WNnk=WN.^(-nk);

xn=(Xk*WNnk)/N;

使用方法如下：

1、序列的傅里叶变换及离散傅里叶变换计算

N=5;

n=0:4;

x=[ones(1,5)]; %产生矩形序列

k=0:999;w=(pi/500)*k;

X=x*(exp(-j*pi/500)).^(n'*k); %计算序列的傅立叶变换

Xe=abs(X);

subplot(3,2,1);stem(n,x);ylabel('x(n)');

subplot(3,2,2);plot(w/pi,Xe);ylabel('|X(ejw)|'); %画出序列的傅立叶变换X=dft(x,N); %进行5点DFT

magX=abs(X);

k=(0:length(magX)'-1)*N/length(magX);

subplot(3,2,3);stem(n,x);ylabel('x(n)');

subplot(3,2,4);stem(k,magX);

axis([0,5,0,5]);ylabel('|X(k)|');

2、序列产生及其DFT

N=20;

n=0:N-1;

x=2*cos(0.25*pi*n)+cos(0.65*pi*n);

subplot(2,1,1);stem(n,x);title('N=20时的信号');

Y=dft(x,N);

k1=0:N-1;w1=2*pi/N*k1;

subplot(2,1,2);stem(w1/pi,abs(Y));title('信号的频谱');

FFT 算法计算矢量xn 的离散傅立叶变换

一维快速傅立叶变换函数 fft

Xk=fft(xn,N) xn 为时域序列向量，N 是DFT 变换区间长度。

当N 大于xn 的长度时，fft 函数自动在xn 后面补零。函数返回xn 的N 点DFT 变换结果向量Xk 。当N 小于xn 的长度时，fft 函数计算xn 的前面N 个元素构成的N 长序列的N 点DFT ，忽略xn 后面的元素。

当x 为矩阵时，Xk 为矩阵xn 的每一列的FFT ,fft 函数按类似的方法处理列长度。 当xn 的长度为2的幂次方时，则fft 采用基2的FFT 算法，否则采用稍慢的混合基算法。 一维快速傅立叶反变换 ifft

xn=ifft(Xk,N) Xk 的N 点IFFT

xn=ifft(Xk) 用于计算矢量Xk 的IFFT

使用方法如下：

M=26;N=32;n=0:M;

xa=0:M/2;

xb=ceil(M/2)-1:-1:0;

xn=[xa,xb];%产生M 长三角波序列x(n)

Xk=fft(xn,512);%512点FFT

X32k=fft(xn,32);%32点FFT

x32n=ifft(X32k);%32点IFFT 得到x32(n)

X16k=X32k(1:2:N);%隔点抽取X32k 得到X16k

x16n=ifft(X16k,N/2); %16点IFFT 得到x16(n)

连续信号的傅里叶变换及其逆变换定义：

(w)(t)e jwt F f dt ∞

--∞=

?

1(t)(w)e 2jwt f F dw π∞-∞=?

Fw=fourier(ft,t,w)%求时域函数ft 的Fourier 变换Fw

Ft=ifourier(Fw,w,t)%求频域函数Fw 的Fourier 逆变换ft

使用方法如下：

syms t w

ut=heaviside(t);%单位阶跃函数

Uw=fourier(ut,t,w);%傅里叶变换

SUt=simple(Uw);%化简

Inv_ut=ifourier(Uw,w,t);%傅里叶逆变换

ezplot(w,Uw)%绘制傅里叶变换曲线

title('傅里叶变换曲线')

grid on

因果序列的Z 变换及其逆变换定义：

X(z)x(n)z

n n ∞-=-∞=

∑ 11x(n)(z)z 2n c X dz j π-=?

Xz=ztrans(xn,n,z)%求时域序列xn 的Z 变换Xz

xn=iztrans(Xz,z,n)%求频域序列Xz 的Z 逆变换xn

使用方法如下：

syms n w T z

xn=sin(w*n*T);%定义序列

Xz=ztrans(xn,n,z);%Z 变换

Inv_xn=iztrans(Xz,z,n);%Z 逆变换

连续信号的单边拉普拉斯变换及其逆变换定义：

0_(s)(t)e st F f dt ∞

-=

? 1(t)(s)e ,02j st j f F ds t j σσπ+∞

-∞

=>? Fs=laplace(ft,t,s)%求时域函数ft 的Laplace 变换Fs

ft=ilaplace(Fs,s,t)%求频域函数Fs 的Laplace 逆变换ft

使用方法如下：

syms t s

syms a positive

ft=dirac(t);%单位冲激函数

Fs=laplace(ft,t,s);%拉普拉斯变换

Inv_ft=ifourier(Fs,s,t);%拉普拉斯逆变换

线性常系数差分方程的递推求解

00(n i)(n i)N M i i

i i a y b x ==-=-∑∑ 0,1a =

求差分方程的零状态响应和全响应函数

yn=filter(B,A,xn) 计算系统对输入信号向量xn 的零状态响应输出信号向量yn ，yn 与xn 长度相等，其中，B 和A 是差分方程的系数向量，即01[b ,b ,,b ]M B = ,其中01a =，如果01a ≠，

则filter 用0a 对系数向量B 和A 归一化。

yn=filter(B,A,xn,xi) 计算系统对输入信号向量xn 的全响应输出信号yn 。其中，xi 是等效初始条件的输入序列，所以xi 是由初始条件确定的。

xi=filtic(B,A,ys,xs) 其中，ys 和xs 是初始条件向量：ys=[y(-1),y(-2),y(-3),…,y(-N)]，xs=[x(-1),x(-2),y(-3),…,x(-M)]。如果xn 是因果序列，则xs=0，调用时可缺省xs 。 使用方法如下：

%调用filter 解差分方程y(n)-ay(n-1)=x(n)

a=0.8;ys=1; %设差分方程系数a=0.8，初始状态：y(-1)=1 xn=[1,zeros(1,30)]; %x(n)=单位脉冲序列，长度N=31

B=1;A=[1,-a]; %差分方程系数

xi=filtic(B,A,ys); %由初始条件计算等效初始条件的输入序列xi yn=filter(B,A,xn,xi); %调用filter 解差分方程，求系统输出函数y(n)

n=0:length(yn)-1;

stem(n,yn,'.')

xlabel('n');ylabel('y(n)')

离散系统函数定义：

00(z)(z)(z)M i i i N i

i

i b z Y H X a z -=-===∑∑ 系统的因果性：其单位脉冲响应h(n)是因果序列，满足n<0时，h(n)=0

或 其系统函数H(z)的收敛域包含无穷点

系统的稳定性：其单位脉冲响应h(n)绝对可和，即

|h(n)|n ∞=-∞<∞∑

或 其系统函数H(z)的收敛域包含单位圆 系统稳定性判定函数 stab(A)

function stab(A)

%stab ：系统稳定性判定函数，A 是H(z)的分母多项式系数向量

disp('系统极点为：')

P=roots(A)%求H(z)的极点，并显示

disp('系统极点模的最大值为：')

M=max(abs(P))%求所有极点模的最大值，并显示

if M<1

disp('系统稳定')

else

disp('系统不稳定')

end

零、 极点及频率响应曲线

zplane(B,A) 绘制出系统函数H(z)的零极点图。其中B 和A 为系统函数H(Z)的分子和分母

多项式系数向量。 假设11(z)(1)(2)z (M 1)z (z)(z)(1)(2)z (N 1)z M N B B B B H A A A A ----++++==++++

则

B=[B(1),B(2),...,B(M+1)], A=[A(1),A(2),...,A(N+1)]

H=freqz(B,A,w) 计算由向量w 指定的数字频率点上数字滤波器H(z)的频率响应H(jw e )，结果存于H 向量中。

[H,w]=freqz(B,A,M) 计算出M 个点上的频率响应，存于H 向量中，M 个频率响应存于w 向量中，频率范围自动均匀设在[0,π]上。

[H,w]=freqz(B,A,M,’whole ’) 频率范围自动均匀设在[0,2π]上。

freqz(B,A) 自动选取512个频率点计算，频率范围自动均匀设在[0,π]上。不带输出向量的freqz 函数将自动绘出幅频响应和相频响应曲线，频率和相位是线性坐标，幅频响应为对数坐标。

使用方法如下：

B=[1,0,0,0,0,0,0,0,-1];A=1;%设置系统函数系数向量B 和A

subplot(2,2,1)

zplane(B,A);%绘制零极点图

title('零极点图')

[H,w]=freqz(B,A);%计算频率响应

subplot(2,2,2)

plot(w/pi,abs(H)/max(abs(H)))%绘制幅频响应曲线

xlabel('\omega/\pi')

ylabel('幅度')

title('幅频特性曲线')

axis([0,1,0,2.5])

grid on

subplot(2,2,4)

plot(w/pi,angle(H))%绘制相频响应曲线

xlabel('\omega/\pi')

ylabel('相位')

title('相频特性曲线')

grid on

实序列分解函数

% 实序列可分解为共扼对称分量 ]n ))-x ((N [x (n )*(1/2)x e c N +=

% 和共扼反对称分量 ]n ))-x ((N -[x (n )*(1/2)xoc N =

function [xec,xoc]=circevod(x)

N=length(x);

n=0:(N-1);

xec=0.5*(x+x(mod(-n,N)+1)); % mod 取余，例如 mod(4,5)=4， mod(-4,5)= mod(5-4,5)=1

xoc=0.5*(x-x(mod(-n,N)+1));

序列的卷积定义：

1212(n)(m)(n m)x (n)*x (n)y x x ∞

-∞=-=∑

序列的卷积函数 convu(x1,n1,x2,n2)

function [y,ny]=convu(x1,n1,x2,n2)

%convu 通用卷积函数

%y 为卷积结果序列向量，ny 是y 的位置向量

%x1和x2是有限长序列

%n1和n2分别是x1和x2的位置向量

nys=n1(1)+n2(1);

nyf=n1(end)+n2(end);

y=conv(x1,x2);

ny=nys:nyf;

用重叠相加法计算序列线性卷积

y=fftfilt(h,x,M) 其中，h 是系统单位脉冲响应向量，x 是输入序列向量，y 是系统的输出序列向量（h 与x 的卷积结果），M 是由用户选择的输入序列x 的分段长度，缺省M 时，默认输入序列x 的分段长度M=512。

使用方法如下：

Lx=41;N=5;M=10; %Lx 为序列x(n)长度

hn=ones(1,N);hn1=[hn,zeros(1,Lx-N)]; %产生h(n),其后补零是为了绘图好看

n=0:L-1;

xn=cos(pi*n/10)+cos(2*pi*n/5); %产生x(n)的Lx 个样值

yn=ftfilt(hn,xn,M) %用重叠相加法计算卷积

连续信号的卷积定义：

1212(t)()f (t )d (t)*(t)f f f f τττ∞

-∞=

-=?

连续信号的卷积函数 sconv(f1,k1,f2,k2,p)

function [f,k]=sconv(f1,k1,f2,k2,p)

%计算连续信号卷积f(t)=f1(t)*f2(t)

%f ：卷积f(t)对应的非零样值向量

%k ：f(t)对应的时间向量

%f1：f1(t)的非零样值向量

%k1：f1(t)对应的时间向量

%f2：f2(t)的非零样值向量

%k2：f2(t)对应的时间向量

%p ：抽样时间间隔

f=conv(f1,f2); %计算f1(t)和f2(t)的卷积和f(t)

f=f*p;

k0=k1(1)+k2(1); %计算f(t)非零样值得起点位置

k3=length(f1)+length(f2)-2;%计算卷积和f(t)的非零样值的宽度

k=k0:p:k3*p;%确定卷积和f(t)非零样值的时间向量

subplot(2,2,1)

plot(k1,f1) %在子图1绘制f1(t)的时域波形图

title('f1(t)')

xlabel('t')

ylabel('f1(t)')

grid on

subplot(2,2,2)

plot(k2,f2)%在子图2绘制f2(t)的时域波形图

title('f2(t)')

xlabel('t')

ylabel('f2(t)')

grid on

subplot(2,2,3)

plot(k,f);%在子图3绘制卷积和f(t)的时域波形图

h=get(gca,'position') ;

h(3)=2.5*h(3);

set(gca,'position',h)%将第三个子图的横坐标范围扩大为原来的2.5倍

title('f(t)=f1(t)*f2(t)')

xlabel('t')

ylabel('f(t)')

grid on

序列的循环移位函数 N m n x n y ))(()(-=

function y=cirshftt(x,m,N)

if length(x)>N

error('N mustbe >= the length of x') %要求移位周期大于信号长度

end

x=[x zeros(1,N-length(x))];

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

n=mod(n-m,N);

y=x(n+1);

时域循环卷积定理：

设有限长序列为)(1n x 和)(2n x ，它们的N 点DFT 分别为)(1k X 和)(2k X ，如果)()()(21k X k X k X ?=，则其IDFT 为两序列的循环卷积1

120()IDFT[()]()(())()N N N m x n X k x m x n m R n -===-∑。

计算两序列的循环卷积函数 1

120(n)[(m)x ((n m))]R (n)L c L L m y x -==-∑

function y=circonv(x1,x2,N)

%循环卷积的长度N 应大于等于max(length(x1),length(x2))

if length(x1)>N

error('N must be >= the length of x1')

end

if length(x2)>N

error('N must be >= the length of x2')

end

x1=[x1 zeros(1,N-length(x1))]; %将序列补0成为N 长序列

x2=[x2 zeros(1,N-length(x2))];

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

x3=x2(mod(-m,N)+1); %该语句的功能相当于序列翻褶，延拓，取主值序列

H=zeros(N,N);

for n=1:1:N %得到x3序列的循环移位矩阵

H(n,:)=cirshftt(x3,n-1,N);

end

y=x1*H' %用矩阵相乘的方法得到结果

信噪比

设纯净信号为)(n s ，噪声信号为)(n v ，带噪信号为)(n x ，则信噪比的定义如式（4-2）所示，单位为dB 。

21

21

0|)(|1|)(|1log 10∑∑-=-==N n N n n v N n s N SNR （4-2） function y=snr(x0,xn)

%计算信噪比，单位dB

%x0是纯净信号,xn 是带噪信号

p0=sum(abs(x0).^2); %sum 函数为求和

y=10*log10(p0/(sum(abs(xn-x0).^2)));

序列的基本运算

%1.单位取样序列x(n)=delta(n-n0) 要求n1<=n0<=n2 function [x,n]=impseq(n0,n1,n2)

n=[n1:n2];

x=[(n-n0)==0];

%2.单位阶跃序列x(n)=u(n-n0) 要求n1<=n0<=n2 function [x,n]=stepseq(n0,n1,n2)

n=[n1:n2];

x=[(n-n0)>=0];

%3.信号加y(n)=x1(n)+x2(n)

function [y,n]=sigadd(x1,n1,x2,n2)

%find函数：找出非零元素的索引号

%x1:第一个序列的值，n1:序列x1的位置向量

%x2:第二个序列的值，n2:序列x2的位置向量

n=min(min(n1),min(n2)):max(max(n1),max(n2));

y1=zeros(1,length(n));

y2=y1;

y1(find((n>=min(n1))&(n<=max(n1))==1))=x1;

y2(find((n>=min(n2))&(n<=max(n2))==1))=x2;

y=y1+y2;

%4.信号乘y(n)=x1(n)*x2(n)

function [y,n]=sigmult(x1,n1,x2,n2)

n=min(min(n1),min(n2)):max(max(n1),max(n2));

y1=zeros(1,length(n));

y2=y1;

y1(find((n>=min(n1))&(n<=max(n1))==1))=x1;

y2(find((n>=min(n2))&(n<=max(n2))==1))=x2;

y=y1.*y2;

%5.移位y(n)=x(n-n0)

function [y,n]=sigshift(x,m,n0)%m：序列x的位置向量n=m+n0;

y=x;

%6.翻褶y(n)=x(-n)

function [y,n]=sigfold(x,n)

y=fliplr(x);

n=-fliplr(n);

特殊的连续函数

单位冲激函数dirac(t)

使用方法如下：

syms t n

dirac(t)

delta=sym('kroneckerDelta(n,0)');%定义单位冲激

单位阶跃函数heaviside(t)

使用方法如下：

syms t

heaviside(t)

录音函数

function record(n,fs,file) %录音函数

%设置n录音时间，采样率fs ,保存文件名file

fprintf('Press any key to start %g seconds of recording...\n',n);

pause;

fprintf('Recording...\n');

y=wavrecord(n*fs,fs,'int16'); %录音

fprintf('Finished recording.\n');

wavwrite(y,fs,file); %保存为文件为file

fprintf('The recording is saved as ');fprintf(file);fprintf('\n');

读取音频

y = wavread(filename) 加载音频文件文件名指定的字符串,返回y的抽样数据。如果不包括一个扩展文件名,wavread默认文件后缀为. wav。

[y, Fs] = wavread(filename) 返回的采样率(Fs)赫兹用于编码中的数据文件。

[y, Fs, nbits] = wavread(filename) 返回每个抽样的比特数(nbits)。

写入音频

wavwrite(y,filename) 把在变量y中的存储数据写入文件名为filename的文件。文件名输入是一个用单引号括起来的字符串。默认数据采样率为8000 Hz和被认为是16位。每一列代表一个单独的数据通道。因此,立体数据应该指定为一个矩阵有两个列。

wavwrite(y,Fs,filename) 把在变量y中的存储数据写入文件名为filename的文件。数据的采样率为Fs赫兹和被认为是16位。

wavwrite(y,Fs,N,filename) 把在变量y中的存储数据写入文件名为filename的文件。数据采样率为Fs赫兹和N位,其中N是8,16、24、32。

播放音频

sound(y,Fs) 以采样率Fs播放音频信号y。如果你不指定采样率,默认为8192赫兹。单通道(单声道)音频,y是一个米长、1列向量,m是音频样本的数量。如果你的系统支持立体声播放,y是一个m-by-2矩阵,第一列对应于左声道,和第二列对应于右声道。声音功能假设y 包含-1和1之间的浮点数,和剪辑范围之外的值。

sound(y,Fs,bits) 指定的位深度(即精度)的样本值。位深度的可能的值取决于可用的音频硬件系统上。大多数平台支持8位或16位的深度。如果你不指定位,声音功能在一个8位深度。

wavplay(y,Fs) 播放存储的音频信号向量y。Fs是整数采样率，单位是Hz。Fs的默认值是11025赫兹。wavplay只支持1 或2声道音频(单声道或立体声)信号。在立体声,y必须是一个两列矩阵。

使用方法如下：

[y,fs]=wavread('sound.wav');%读取音频文件

L=length(y);%采样点数

T=L/fs;%采样时间

fprintf('按任意键播放原录音及波形\n');

pause;

plot(y); %画图

title('原信号8000Hz波形')

wavplay(y,fs);

fprintf('录音频率：%g Hz\n',fs);

fs1=16000;

L1=T*fs1;

B1=L/L1;

y1=0*(1:L1);

for i=1:L1

y1(i)=y(ceil(i*B1));%重新采样

end

wavwrite(y1,fs1,'sound1.wav'); %保存为文件为sound1.wav

y1=wavread('sound1.wav');%读取音频文件

fprintf('按任意键播放采样频率%g Hz录音及波形\n',fs1);

pause;

plot(y1); %画图

title('采样信号16000Hz波形')

wavplay(y1,fs1);

fs2=8000;

L2=T*fs2;

B2=L/L2;

y2=0*(1:L2);

for i=1:L2

y2(i)=y(ceil(i*B2));%重新采样

end

wavwrite(y2,fs2,'sound2.wav'); %保存为文件为sound2.wav

y2=wavread('sound2.wav');%读取音频文件

fprintf('按任意键播放采样频率 %g Hz 录音及波形\n',fs2);

pause;

plot(y2); %画图

title('采样信号8000Hz 波形')

wavplay(y2,fs2);

示例

1、令||1000)(t a e t x -=，绘制其傅立叶变换)(Ωj X a 。用不同频率对其进行采样，分别画出)(ωj e X 。

Dt=0.00005; %步长为0.00005s

t=-0.005:Dt:0.005;

xa=exp(-1000*abs(t));%取时间从-0.005s 到0.005s 这段模拟信号 Wmax=2*pi*2000; %信号最高频率为2 *2000

K=500; %频域正半轴取500个点进行计算

k=0:1:K;

W=k*Wmax/K;%求模拟角频率

Xa=xa*exp(-1j*t'*W)*Dt;%计算连续时间傅立叶变换（利用矩阵运算实现） Xa=real(Xa);%取实部

W=[-fliplr(W),W(2:501)];%将角频率范围扩展为从-到+

Xa=[fliplr(Xa),Xa(2:501)];

subplot(3,2,1);

plot(t*1000,xa); %画出模拟信号，横坐标为时间（毫秒），纵坐标为幅度 xlabel('time(millisecond)');ylabel('xa(t)');

title('anolog signal');

subplot(3,2,2);

plot(W/(2*pi*1000),Xa*1000);%画出连续时间傅立叶变换

xlabel('frequency(kHZ)'); %横坐标为频率（kHz ）

ylabel('xa(jw)');%纵坐标为幅度

title('FT');

%下面为采样频率1kHz时的程序

T1=0.001; %采样间隔为

n=-25:1:25;

x1=exp(-1000*abs(n*T1)); %离散时间信号

K=500;k=0:1:K;w1=pi*k/K; %w为数字频率

X1=x1*exp(-1i*n'*w1);%计算离散时间傅立叶变换（序列的傅立叶变换）

X1=real(X1);

w1=[-fliplr(w1),w1(2:K+1)];

X1=[fliplr(X1),X1(2:K+1)];

subplot(3,2,3);

stem(n*T1*1000,x1,'.','MarkerSize',15); %画出采样信号（离散时间信号）xlabel('time(millisecond)');

ylabel('x1(n)');

title('discrete signal');

subplot(3,2,4);

plot(w1/pi,X1);%画出离散时间傅立叶变换

xlabel('frequency(radian)');%横坐标为弧度

ylabel('x1(jw)');title('DTFT');