close all;
clear all;
clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%---------輸入信號-------------------------------
N=512;
n=0:1:N-1;
f1=50;
f2=10;
f3=500;
fs=1400;

x1=sin(f1*2*pi*n./fs); % x1=cos(2*pi*n*f1/1000+5*pi*n.*n/1000);%
x2=sin(f2*2*pi*n./fs); % x2=[(n-100)==0];
x3=sin(f3*2*pi*n./fs); % x3=[(n-20)>=0].*[(n-22)<=0];
xn=sin(f1*2*pi*n./fs)+sin(f2*2*pi*n./fs)+sin(f3*2*pi*n./fs);
%xn=1*cos(2*pi*n*10/1000+0.01*pi*n.*n/1000);
figure
subplot(4,2,1);plot(n,xn/max(abs(xn)));axis([0 512 -1 1]);title('時域信號');xlabel('xn=x1+x2+x3');grid
xk=fft(xn);
subplot(4,2,2);plot(n,xk/max(abs(xk)));axis([0 512 -1 1]);title('頻域信號');xlabel('xk');

subplot(4,2,3);plot(n,x1/max(abs(x1)));axis([0 512 -1 1]);xlabel('x1');grid
xk1=fft(x1);
subplot(4,2,4);plot(n,xk1/max(abs(xk1)));axis([0 512 -1 1]);xlabel('xk1');

subplot(4,2,5);plot(n,x2/max(abs(x2)));axis([0 512 -1 1]);xlabel('x2');grid
xk2=fft(x2);
subplot(4,2,6);plot(n,xk2/max(abs(xk2)));axis([0 512 -1 1]);xlabel('xk2');

subplot(4,2,7);plot(n,x3/max(abs(x3)));axis([0 512 -1 1]);xlabel('x3');grid
xk3=fft(x3);
subplot(4,2,8);plot(n,xk3/max(abs(xk3)));axis([0 512 -1 1]);xlabel('xk3');
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%----------polyphase filter---------------------
N=512;   % 濾波器階數
M=32;    % 32個信道
K=N/M;   % 子信道階數
Bw=800;  % 800MHz
Bw0=round(1.1*Bw/32);
Fstart=0;
Fstop=Fstart+Bw0;

% hh=fir1(N-1,Fstop/Bw,hamming(N));
hh=remez(N-1,[0 96 160 192*32]/(192*32),[1 1 0 0],[1 13]);
% freqz(hh0,1,N)
% pause

figure(2)
Q=1024;
t=-0.5:1/Q:.5-1/Q;
subplot(2,1,1)
plot(hh/max(abs(hh)),'r');
hold on;
plot(n,xn/max(abs(xn)),'b');

subplot(2,1,2)
plot(t*64,fftshift(20*log10(abs(fft(hh,Q)))),'r');
axis([-32 32 -80 20]);
title('Impulse Response: Prototype Filter');
xlabel('Frequency (kHz)');
ylabel('Log-Magnitude (dB)');
hold on

nn=-N/2+0.5:N/2-0.5;
hh1=hh.*exp(j*2*pi.*nn*1/M);
plot(t*64,fftshift(20*log10(abs(fft(hh1,Q)))),'r--');
hold on
hh2=hh.*exp(j*2*pi.*nn*2/M);
plot(t*64,fftshift(20*log10(abs(fft(hh2,Q)))),'r--');
hold on
hh3=hh.*exp(j*2*pi.*nn*3/M);
plot(t*64,fftshift(20*log10(abs(fft(hh3,Q)))),'r--');
hold on
hh16=hh.*exp(j*2*pi.*nn*16/M);
plot(t*64,fftshift(20*log10(abs(fft(hh16,Q)))),'r--');
hold on

hh31=hh.*exp(j*2*pi.*nn*31/M);
plot(t*64,fftshift(20*log10(abs(fft(hh31,Q)))),'r--');
hold on

ww=kaiser(1024,8)';
ww=ww/sum(ww);
plot(t*64,fftshift(20*log10(abs(fft(xn,Q).*ww))),'b')
pause
%%%%%%%%%%%%%%%%%%%---------------------------------------------------

figure(3)
subplot(2,1,1);
plot(hh);hold on
grid
title('Impulse Response: Prototype Filter')
xlabel('Normalized time nT/T')
ylabel('Amplitude')
n=1:16;
nn=-K/2:K/2-1;

hh0=hh(1+(n-1).*M); plot(hh0,'g');hold on
hh1=hh(2+(n-1).*M); % plot(hh1,'g');hold on
hh2=hh(3+(n-1).*M); % plot(hh2,'g');hold on
hh3=hh(4+(n-1).*M); % plot(hh3,'g');hold on
hh4=hh(5+(n-1).*M); % plot(hh4,'g');hold on
hh5=hh(6+(n-1).*M); % plot(hh5,'g');hold on
hh6=hh(7+(n-1).*M); % plot(hh6,'g');hold on
hh7=hh(8+(n-1).*M); % plot(hh7,'g');hold on
hh8=hh(9+(n-1).*M); % plot(hh8,'g');hold on
hh9=hh(10+(n-1).*M); % plot(hh9,'g');hold on
hh10=hh(11+(n-1).*M); % plot(hh10,'g');hold on
hh11=hh(12+(n-1).*M); % plot(hh11,'g');hold on
hh12=hh(13+(n-1).*M); % plot(hh12,'g');hold on
hh13=hh(14+(n-1).*M); % plot(hh13,'g');hold on
hh14=hh(15+(n-1).*M); % plot(hh14,'g');hold on
hh15=hh(16+(n-1).*M); % plot(hh15,'g');hold on
hh16=hh(17+(n-1).*M); % plot(hh16,'g');hold on
hh17=hh(18+(n-1).*M); % plot(hh17,'g');hold on
hh18=hh(19+(n-1).*M); % plot(hh18,'g');hold on
hh19=hh(20+(n-1).*M); % plot(hh19,'g');hold on
hh20=hh(21+(n-1).*M); % plot(hh20,'g');hold on
hh21=hh(22+(n-1).*M); % plot(hh21,'g');hold on
hh22=hh(23+(n-1).*M); % plot(hh22,'g');hold on
hh23=hh(24+(n-1).*M); % plot(hh23,'g');hold on
hh24=hh(25+(n-1).*M); % plot(hh24,'g');hold on
hh25=hh(26+(n-1).*M); % plot(hh25,'g');hold on
hh26=hh(27+(n-1).*M); % plot(hh26,'g');hold on
hh27=hh(28+(n-1).*M); % plot(hh27,'g');hold on
hh28=hh(29+(n-1).*M); % plot(hh28,'g');hold on
hh29=hh(30+(n-1).*M); % plot(hh29,'g');hold on
hh30=hh(31+(n-1).*M); % plot(hh30,'g');hold on
hh31=hh(32+(n-1).*M); % plot(hh31,'g');hold on
% 
% subplot(2,1,2);
% plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(hh,1024)))),'b');axis([-16 16 -100 1]);hold on;
% plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(hh0,1024)))),'r');hold on;
% plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(hh1,1024)))),'m--');hold on;
% plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(hh2,1024)))),'m--');hold on;
% 
% plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(hh15,1024)))),'g--');hold on;

%--------------------------------------------------------------------------
hh0=hh0;
hh1=hh1.*exp(j*2*pi.*nn*1/M);
hh2=hh2.*exp(j*2*pi.*nn*2/M);
hh3=hh3.*exp(j*2*pi.*nn*3/M);
hh4=hh4.*exp(j*2*pi.*nn*4/M);
hh5=hh5.*exp(j*2*pi.*nn*5/M);
hh6=hh6.*exp(j*2*pi.*nn*6/M);
hh7=hh7.*exp(j*2*pi.*nn*7/M);
hh8=hh8.*exp(j*2*pi.*nn*8/M);
hh9=hh9.*exp(j*2*pi.*nn*9/M);
hh10=hh10.*exp(j*2*pi.*nn*10/M);
hh11=hh11.*exp(j*2*pi.*nn*11/M);
hh12=hh12.*exp(j*2*pi.*nn*12/M);
hh13=hh13.*exp(j*2*pi.*nn*13/M);
hh14=hh14.*exp(j*2*pi.*nn*14/M);
hh15=hh15.*exp(j*2*pi.*nn*15/M);
hh16=hh16.*exp(j*2*pi.*nn*16/M);
hh17=hh17.*exp(j*2*pi.*nn*17/M);
hh18=hh18.*exp(j*2*pi.*nn*18/M);
hh19=hh19.*exp(j*2*pi.*nn*19/M);
hh20=hh20.*exp(j*2*pi.*nn*20/M);
hh21=hh21.*exp(j*2*pi.*nn*21/M);
hh22=hh22.*exp(j*2*pi.*nn*22/M);
hh23=hh23.*exp(j*2*pi.*nn*23/M);
hh24=hh24.*exp(j*2*pi.*nn*24/M);
hh25=hh25.*exp(j*2*pi.*nn*25/M);
hh26=hh26.*exp(j*2*pi.*nn*26/M);
hh27=hh27.*exp(j*2*pi.*nn*27/M);
hh28=hh28.*exp(j*2*pi.*nn*28/M);
hh29=hh29.*exp(j*2*pi.*nn*29/M);
hh30=hh30.*exp(j*2*pi.*nn*30/M);
hh30=hh31.*exp(j*2*pi.*nn*31/M);

subplot(2,1,2);
plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(hh,1024)))),'b');axis([-16 16 -100 1]);hold on;
plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(hh0,1024)))),'r');hold on;
plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(hh1,1024)))),'m--');hold on;
plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(hh2,1024)))),'m--');hold on;

plot((-0.5:1/1024:.5-1/1024)*32,fftshift(20*log10(abs(fft(hh15,1024)))),'g--');hold on;
pause
%%%%%%%%%%%%%%%%%---------------------濾波后信號---------------------------

y1=conv(x1,hh);yk1=fft(y1,N);
y2=conv(x2,hh);yk2=fft(y2,N);
y3=conv(x3,hh);yk3=fft(y3,N);
yn=conv(xn,hh);yk=fft(yn,N);

figure(4)
subplot(4,2,1);plot(yn);title('512階濾波后信號');subplot(4,2,2);plot(abs(yk));title('512階濾波后信號頻譜');
subplot(4,2,3);plot(y1);subplot(4,2,4);plot(abs(yk1));
subplot(4,2,5);plot(y2);subplot(4,2,6);plot(abs(yk2));
subplot(4,2,7);plot(y3);subplot(4,2,8);plot(abs(yk3));

pause

y1=conv(x1,hh0);yk1=fft(y1,N);
y2=conv(x2,hh0);yk2=fft(y2,N);
y3=conv(x3,hh0);yk3=fft(y3,N);
yn=conv(xn,hh0);yk=fft(yn,N);

figure(5)
subplot(4,2,1);plot(yn);title('16階濾波后信號');subplot(4,2,2);plot(abs(yk));title('16階濾波后信號頻譜');
subplot(4,2,3);plot(y1);subplot(4,2,4);plot(abs(yk1));
subplot(4,2,5);plot(y2);subplot(4,2,6);plot(abs(yk2));
subplot(4,2,7);plot(y3);subplot(4,2,8);plot(abs(yk3));

pause

y1=conv(x1,hh5);yk1=fft(y1,N);
y2=conv(x2,hh5);yk2=fft(y2,N);
y3=conv(x3,hh5);yk3=fft(y3,N);
yn=conv(xn,hh5);yk=fft(yn,N);

figure(6)
subplot(4,2,1);plot(abs(yn));title('16階濾波后信號');subplot(4,2,2);plot(abs(yk));title('16階濾波后信號頻譜');
subplot(4,2,3);plot(abs(y1));subplot(4,2,4);plot(abs(yk1));
subplot(4,2,5);plot(abs(y2));subplot(4,2,6);plot(abs(yk2));
subplot(4,2,7);plot(abs(y3));subplot(4,2,8);plot(abs(yk3));