a=[1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0];f=[1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0];
I=8;n1=200;Fs=25.0;fs=Fs*2*I;kf=12.5;
for k=1:I
    for r=1:(n1+I)
        m(k,r)=a(k)*(1.0+0.5*cos(2*pi*f(k)/Fs*(r-1)));
    end 
end
[n0,f0,m0,w]=remezord([8 12.5],[1 0],[0.001 0.001],fs);
b=remez(287,f0,m0,w);
figure(1)
plot(20*log10(abs(fft(b))));
grid on;
for r=1:36
    for k=1:I
        h(k,r)=b((r-1)*I+k);
    end
end
for r=1:(n1+I)
    for k=1:I
        mk(k)=m(k,r);
    end
    mfft=fft(mk);
    for k=1:I
        x0(k,r)=mfft(k)*exp(j*pi/(2*I)*(k-1));
    end
end
for r=1:(n1+I)
    for k=1:I
        x00(k,(2*r-1))=x0(k,r);
        x00(k,2*r)=0;
    end
end
for k=1:I
    for r=1:(n1+I)
        x00k(r)=x00(k,r);
    end
    for r=1:36
        hk(r)=h(k,r);
    end
    y0=conv(x00k,hk);
    for r=1:n1
        y0(r)=y0(r)*exp(j*pi/2*(r-1));
    end
    for r=1:n1
        y(k,r)=y0(r+36);%*(-1)^(k-1)*exp(j*pi/I*(k-1));
    end
end
for k=1:I
    for n=1:n1*I
        if mod((n-1),I)==0
            y00(k,n)=y(k,(n-1)/I+1);
        else y00(k,n)=0.0;
        end
    end
end
for n=I:(n1*I-I)
    yout(n-I+1)=y00(1,n)+y00(2,n-1)+y00(3,n-2)+y00(4,n-3)+y00(5,n-4)+y00(6,n-5)+y00(7,n-6)+y00(8,n-7);
end
point=512;
yy(1:point)=yout(101:(100+point));
for n=1:point
    yy(n)=(yy(n)+0.001*randn)*(0.42323-0.49755*cos(2*pi*(n-1)/point)+0.07922*cos(4*pi*(n-1)/point));
   l(n)=fs/point*(n-1);
end
yy1=real(yy);
pp1=abs(fft(yy1));
ppm1=max(pp1);
figure(3)
plot(l(1:256),20*log10(pp1(1:256)/ppm1));
grid on;
pp=abs(fft(yy));
ppm=max(pp);
figure(2)
plot(l,20*log10(pp/ppm));