ex0806.m

通信课程仿真

%  [ex8.6] 用FFT分析周期矩形脉冲及其频谱
%  [CTFSfft.m]    这是计算周期矩形脉冲双边频谱的程序。
function [A,B,C,fn,t,x]=CTFSfft(T,M,Nf)
% 利用FFT，计算[0,T}区间上定义的时间波形的Fourier级数展开系数A,B和频谱C,fn 。
%	T		时间波形周期
%	M		用作2的幂次
%	Nf		输出谐波的阶次，决定A,B的长度为(Nf+1)。Nf不要超2^(M-1)。
%	A,B		分别是Fourier级数中cos,sin展开项的系数。A(1)是直流量。
%	C		是定义在[-fs/2,fs/2]上的频谱
%	t,x		是原时间波形数据对
T=5;
if (nargin<2 | isempty(M));M=8;end										
if nargin<3;Nf=input('plear Input 所需展开的最高谐波次数:Nf=');end	
N=2^M;     %  使总采样点是2的整数倍
f=1/T;     % 被变换函数的频率
w0=2*pi*f; 
dt=T/N;     % 时间分辨率
n=0:1:(N-1);  % 采样点序列
t=n*dt;       % 采样时间序列
x=time_fun(t,T);	    % 被变换时间函数的采样序列, 调用CTFSfft.m的子函数		
W=fft(x);     % 给出n=0,1,…..,N-1上的DFT数据值
cn=W/N;       % 据式(3.1.2)计算n=0,1,…,N-1上的FS系数
z_cn=find(abs(cn)<1.0e-10);    % 寻找有限字长运算而产生(原应为0)的"小"复数
cn(z_cn)=zeros(length(z_cn),1); % 强制那些"小"复数为0						
cn_SH=fftshift(cn);	% 据式(3.1.3)计算n=-N/2,…,-1,0,1,…..,(N/2)-1上的FS系数
C=[cn_SH cn_SH(1)];  %  形成关于0对称的（N+1）个FS系数
A(1)=C(N/2+1);
A(2:N/2+1)=2*real(C((N/2+2):end));
B(2:N/2+1)=-2*imag(C((N/2+2):end));
if Nf>N/2;error(['第三输入宗量 Nf 应小于 ' int2str(N/2-1)]);end
A(Nf+2:end)=[];
B(Nf+2:end)=[];
n1=-N/2*2*pi/T:2*pi/T:N/2*2*pi/T;   % 产生总点数为（N+1）个关于0对称的FS系数
fn=n1*f;         % 关于0对称的频率分度序列
y=time_fun_e(t)   % 调用CTFSfft.m的子函数
subplot(2,1,2),
stem(n1,abs(C));
title('周期矩形脉冲频谱')
 axis([-150,150,0,0.12])
hold off
%-------------------------
function   x=time_fun(t,T)
% 该函数是CTFSfft.m的子函数。
%	t		是时间数组
%	T		是周期
A=1;tao=1;T=5;
x=zeros(size(t));ii=find(t>=-tao /2& t<=tao/2);
x(ii)=ones(size(ii)).*A;x(t==0)=1;

%------------------------------
function y=time_fun_e(t)
% 该函数是CTFSfft.m的子函数。它由符号变量和表达式写成。
%	t		是时间数组
%	T		是周期   duty=tao/T=0.2
T=5;t=-2*T:0.01:2*T;tao=T/5;
y=rectpuls(t,1);     %产生一个宽度tao=1的矩形脉冲
subplot(2,2,2)
plot(t,y)
hold on
y=rectpuls(t-5,1);   %产生一个宽度tao=1的矩形脉,中心位置在t=5处
plot(t,y)   
hold on
y=rectpuls(t+5,1);  %产生一个宽度tao=1的矩形脉,中心位置在t=-5处
plot(t,y)
title('周期为T=5，脉宽tao=1的矩形脉冲')
axis([-10,10,0,1.2])