gibbs.m

Gibbs现象在MATLAB上的仿真程序。

close all
clear all
T1=2;    T01=4*T1;T02=8*T1;T03=16*T1;    w01=2*pi/T01;w02=2*pi/T02;w03=2*pi/T03;

%%%%%%%%%%%%%%%%%%%%%周期方波的频谱图%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k=-80:80;
a1=sin(2*k*pi*T1/T01)./(pi*k);    a1(1,81)=2*T1/T01;  %%%%%矩阵的点除./

a2=sin(2*k*pi*T1/T02)./(pi*k);    a2(1,81)=2*T1/T02;

a3=sin(2*k*pi*T1/T03)./(pi*k);    a3(1,81)=2*T1/T03;

stem(k,a1);   grid;xlabel('kw0');ylabel('ak');title('不同基波周期的周期方波的谱线图');axis([-25,25,-0.2,0.5]);%gtext('T0=4T1');
figure;
stem(k,a2);   grid;xlabel('kw0');ylabel('ak');title('不同基波周期的周期方波的谱线图');axis([-20,20,-0.2,0.5]);%gtext('T0=8T1');
figure;
stem(k,a3);   grid;xlabel('kw0');ylabel('ak');title('不同基波周期的周期方波的谱线图');axis([-20,20,-0.2,0.5]);%gtext('T0=16T1');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Gibbs现象%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
figure;
T1=2;    T0=4*T1;  w0=2*pi/T0;

[u,t] = gensig('square',T01,30,0.01);     %产生方波信号数据u
%plot(t-3*T1,u);     axis([-3*T1,3*T1,-0.5,1.5]);

N1=1;N2=3;N3=7;N4=19;N5=79;    %N项近似！！！
%M1=-N1:N1;M2=-N2:N2;M3=-N3:N3;M4=-N4:N4;M5=N5:N5;    %用矩阵的运算代替for循环

xN1=a1(1,81)+2*a1(1,82)*exp(j*w0*t);    %%%%%%%%%%%%%%%%求N1项近似的函数和xN（t）%%%%%%%%%%
for i=2:N1 xN1=xN1+2*a1(1,81+i)*exp(j*i*w01*t);
end;
plot(t-3*T1,u,'m',t-T0*T1,xN1,'g',t-3*T1,max(xN1),'r');axis([-3*T1,3*T1,-0.5,1.5]);grid;title('Gibbs现象:N=1');xlabel('t');ylabel('x(t)和xN(t)');
disp('y1=(max(xN1)-1)/1.0');y1=(max(xN1)-1)/1.0

figure;
xN2=a1(1,81)+a1(1,82)*exp(-j*w0*t)+a1(1,82)*exp(j*w0*t);    %%%%%%%%%%%%%%%%求N2项近似的函数和xN（t）%%%%%%%%%%
%xN2=xN2+2*a1(1,82)*exp(j*w0*t)+2*a1(1,83)*exp(j*2*w0*t)+2*a1(1,84)*exp(j*3*w0*t);
for i=2:N2 xN2=xN2+a1(1,81+i)*exp(-j*i*w0*t)+a1(1,81+i)*exp(j*i*w0*t);
end;
plot(t-3*T1,u,'m',t-T0*T1,xN2,'g',t-3*T1,max(xN2),'r');axis([-3*T1,3*T1,-0.5,1.5]);grid;title('Gibbs现象:N=3');xlabel('t');ylabel('x(t)和xN(t)');
disp('y2=(max(xN2)-1)/1.0');y2=(max(xN2)-1)/1.0

figure;
xN3=a1(1,81)+a1(1,82)*exp(-j*w0*t)+a1(1,82)*exp(j*w0*t);    %%%%%%%%%%%%%%%%求N3项近似的函数和xN（t）%%%%%%%%%%
for i=2:N3 xN3=xN3+a1(1,81+i)*exp(-j*i*w0*t)+a1(1,81+i)*exp(j*i*w0*t);
end;
plot(t-3*T1,u,'m',t-T0*T1,xN3,'g',t-3*T1,max(xN3),'r');axis([-3*T1,3*T1,-0.5,1.5]);grid;title('Gibbs现象:N=7');xlabel('t');ylabel('x(t)和xN(t)');
disp('y3=(max(xN3)-1)/1.0');y3=(max(xN3)-1)/1.0

figure;
xN4=a1(1,81)+a1(1,82)*exp(-j*w0*t)+a1(1,82)*exp(j*w0*t);    %%%%%%%%%%%%%%%%求N4项近似的函数和xN（t）%%%%%%%%%%
for i=2:N4 xN4=xN4+a1(1,81+i)*exp(-j*i*w0*t)+a1(1,81+i)*exp(j*i*w0*t);
end;
plot(t-3*T1,u,'m',t-T0*T1,xN4,'g',t-3*T1,max(xN4),'r');axis([-3*T1,3*T1,-0.5,1.5]);grid;title('Gibbs现象:N=19');xlabel('t');ylabel('x(t)和xN(t)');
disp('y4=(max(xN4)-1)/1.0');y4=(max(xN4)-1)/1.0

figure;
xN5=a1(1,81)+a1(1,82)*exp(-j*w0*t)+a1(1,82)*exp(j*w0*t);    %%%%%%%%%%%%%%%%求N5项近似的函数和xN（t）%%%%%%%%%%
for i=2:N5 xN5=xN5+a1(1,81+i)*exp(-j*i*w0*t)+a1(1,81+i)*exp(j*i*w0*t);
end;
plot(t-3*T1,u,'m',t-T0*T1,xN5,'g',t-3*T1,max(xN5),'r');axis([-3*T1,3*T1,-0.5,1.5]);grid;title('Gibbs现象:N=79');xlabel('t');ylabel('x(t)和xN(t)');
disp('y5=(max(xN5)-1)/1.0');y5=(max(xN5)-1)/1.0