

N=512;
t=0:N-1;

p=3;
q=5;

rolloff=.1;
rejection=10^-3;

r=p/q;
NN=ceil(r*N);
tt=0:NN-1;

err1=zeros(N/2-1,1);
err2=zeros(N/2-1,1);

reject1=zeros(N/2-1,1);
reject2=zeros(N/2-1,1);

for n = 0:N/2-1,
  phi0=2*pi*rand;
  f0=n/N;
  x=sin(2*pi*f0*t' + phi0);
  y1=resample(x,p,q);
  [y2 h2]=my_resample(x,p,q);
  t2=ceil((length(h2)-1)/2/q);

  xx=sin(2*pi*f0/r*tt' + phi0);
  
  err1(n+1)=norm(y1-xx(1:length(y1)));
  err2(n+1)=norm(y2-xx);
  reject1(n+1)=max(abs(y1));
  idx=t2+1:NN-t2;
  reject2(n+1)=max(abs(y2(idx)));  
end;

f=(0:N/2-1)/N;

clf;

subplot(2,2,1);
plot(f/r,err1,'rx-;OctaveForge;');
hold on;
plot(f/r,err2,'g+-;my resample;');
m=max([err1;err2]);
plot((1-rolloff/2)*[.5 .5],[0 m],'--');
plot([.5 .5],[0 m],'--');
plot((1+rolloff/2)*[.5 .5],[0 m],'--');
hold off;
axis([0 0.55]);
title("(a) reconst error and specif for roll-off of anti-aliasing filter");
xlabel("frequency");
ylabel("RMS error");


subplot(2,2,2);
semilogy(f/r,reject1,'rx-;OctaveForge;');
hold on;
semilogy(f/r,reject2,'g+-;my resample;');
semilogy([.5 .5/r],[rejection rejection]);
m=max([reject1;reject2]);
plot((1-rolloff/2)*[.5 .5],[0 m],'--');
plot([.5 .5],[0 m],'--');
plot((1+rolloff/2)*[.5 .5],[0 m],'--');
hold off;
axis([0.5 0.5/r]);
xlabel("frequency");
ylabel("rejection");
title("(b) specif for rejection of anti-aliasing filter");

n=27;
phi0=2*pi*rand;
f0=n/N;
x=sin(2*pi*f0*t' + phi0);
y1=resample(x,p,q);
y2=my_resample(x,p,q);
xx=sin(2*pi*f0/r*tt' + phi0);

subplot(2,2,3);
plot(y1,'rx-;OctaveForge;');
hold on;
plot(y2,"g+;my resample;");
plot(xx,"b;exact;");
hold off;
title(sprintf("(c) f=%2.2f",f0/r));
xlabel("time");
ylabel("data");
axis([1 100])

n=200;
phi0=2*pi*rand;
f0=n/N;
x=sin(2*pi*f0*t' + phi0);
y1=resample(x,p,q);
y2=my_resample(x,p,q);

subplot(2,2,4);
plot(y1,'rx-;OctaveForge;');
hold on;
plot(y2,"g+;my resample;");
hold off;
title(sprintf("(d) f=%2.2f",f0/r));
xlabel("time");
ylabel("data");

print("-deps","resample_test.eps");