function test_pfb_fft_business
%put freq and phase here
phase=300/180*pi;
%freq = 86119800; %2600*33123
%freq = 110822400; %2600*64*666
freq = 86138000; %2600*33130
%generated signal for ADC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% amp_diff and DC here is just try to simulate the situation when the I and
% Q stream have different property. Right now, we can leave the amp_diff =1
% and DC = 0, which mean there is not amplitude and phase difference (other than I and Q)
amp_diff=1;
DC = 0;
step = 1/(340787200);
x = [0:step:(2^19-1)*step];
ADC_noise =60;
y = 2048*0.7*sin(x*2*pi*freq+phase);
z = 2048*0.7*cos(x*2*pi*freq+phase);
y_new = 2048*amp_diff*0.7*sin(x*2*pi*freq+phase)+2048*amp_diff*DC;
z_new = 2048*0.7*cos(x*2*pi*freq+phase)+2048*DC;
%z = sin(x*2*pi/4x
y = round(y);
z = round(z);
y_new = round(y_new);
z_new = round(z_new);
%plot(20*log10(abs(fft(complex(y,-z)))))
y_nois = awgn(y,ADC_noise,'measured');
z_nois = awgn(z,ADC_noise,'measured');
y_new_nois = awgn(y_new,ADC_noise,'measured');
z_new_nois = awgn(z_new,ADC_noise,'measured');
y_nois = round(y_nois);
z_nois = round(z_nois);
y_new_nois = round(y_new_nois);
z_new_nois = round(z_new_nois);
%hold on
%plot(z,'red')
c_nois=complex(z_nois,-1*y_nois);
c_new_nois=complex(z_new_nois,-1*y_new_nois);
size(c_new_nois)
%size is 1 131072*4
temp_c_new_nois = c_new_nois(1,1:2^17);
temp_c_new_nois_fft = ifft(temp_c_new_nois);
figure
plot(20*log10(abs(temp_c_new_nois_fft)),'.r')
title('matlab 2^17 FFT output')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ADD PFB_FIR EFFECT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PFBSize = 11;
TotalTaps = 4;
WindowType = 'hamming';
n_inputs = 0;
nput = 0;
fwidth =1;
a1 = 1;
a2=2;
a3=3;
a4=4;
alltaps = TotalTaps*2^PFBSize;
windowval = transpose(window(WindowType, alltaps));
% plot(windowval)
total_coeffs = windowval .* sinc(fwidth*([0:alltaps-1]/(2^PFBSize)-TotalTaps/2));
%in this case, buf = total_coeff
for i=1:alltaps/2^n_inputs,
buf(i)=total_coeffs((i-1)*2^n_inputs + nput + 1);
end
coeff1 = buf((a1-1)*2^(PFBSize-n_inputs)+1 : a1*2^(PFBSize-n_inputs));
coeff2 = buf((a2-1)*2^(PFBSize-n_inputs)+1 : a2*2^(PFBSize-n_inputs));
coeff3 = buf((a3-1)*2^(PFBSize-n_inputs)+1 : a3*2^(PFBSize-n_inputs));
coeff4 = buf((a4-1)*2^(PFBSize-n_inputs)+1 : a4*2^(PFBSize-n_inputs));
size(coeff1)
for i=1:1:64
for j=1:1:2048
c_new_nois_fir(1,(j+2048*(i-1))) = c_new_nois(1,(j+2048*(i-1)))*coeff1(j)+ c_new_nois(1,(j+2048*(i-1))+2048)*coeff2(j)+c_new_nois(1,(j+2048*(i-1))+4096)*coeff3(j)+c_new_nois(1,(j+2048*(i-1))+6144)*coeff4(j);
end
i
end
size(c_new_nois_fir)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c_new_nois_original = c_new_nois;
c_new_nois = c_new_nois_fir;
%c=y;
%plot(20*log10(abs(fft(c_new))))
%figure
%plot(20*log10(abs(fft(c_new_nois))),'.r')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do the first FFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%maybe add lower frequency ( or so call 1/f) type noise later on
fft_matrix_new = zeros(64,2048);
for i=1:1:64
fft_first_new = (ifft(c_new_nois(1,(2048*(i-1)+1):2048*i),2^11));
fft_matrix_new(i,:) = fft_first_new;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%transpose and 2nd FFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fft_second_new = zeros(1,131072);
for i=1:1:2048
fft_second_new(1,(64*(i-1)+1):64*i)=ifft(fft_matrix_new(:,i),64);
end
figure
plot(20*log10(abs(fft_matrix_new(1,:))),'.r')
title('first PFB output')
figure
plot(20*log10(abs(fft_second_new)),'.b')
title('final FFT output')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bin_pos = freq/2600 + 1;
sprintf('original phase is')
temp_c_new_nois_fft(1,bin_pos)
angle(temp_c_new_nois_fft(1,bin_pos))/pi*180
sprintf('pfb + fft phase is')
fft_second_new(1,bin_pos)
angle(fft_second_new(1,bin_pos))/pi*180
angle(fft_second_new(1,bin_pos+64))/pi*180
% freq =86119800, phase is set to 60, direct 2^17 matlab fft get -60 and PFB+FFT get -93.7598
% freq = 86119800, phase is set to 20, direct 2^17 matlab fft get -20 and PFB+FFT get -53.7598
% freq = 86138000, phase is set to 20, direct 2^17 matlab fft get -20 and PFB+FFT get -132.5052
% freq = 86138000, phase is set to 300, direct 2^17 matlab fft get 60 and PFB+FFT get -52.5052
%%%%% for future simulation %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%select Bin
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FIR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%decimation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% xx = abs(fft(c,2^17))/(2^17);
% for i=1:1:2^17
% if xx(i)==0
% xx(i)=1/(2^17);
% end
% end
%
% xx_new = abs(fft(c_new,2^17))/(2^17);
% for i=1:1:2^17
% if xx_new(i)==0
% xx_new(i)=1/(2^17);
% end
% end
%
%
% yy = real(fft(c,2^17))/(2^17);
% for i=1:1:2^17
% if yy(i)==0
% yy(i)=1/(2^17);
% end
% end
%
% yy_new = real(fft(c_new,2^17))/(2^17);
% for i=1:1:2^17
% if yy_new(i)==0
% yy_new(i)=1/(2^17);
% end
% end
%
%
%
%
% zz = imag(fft(c,2^17))/(2^17);
% for i=1:1:2^17
% if zz(i)==0
% zz(i)=1/(2^17);
% end
% end
%
% zz_new = imag(fft(c_new,2^17))/(2^17);
% for i=1:1:2^17
% if zz_new(i)==0
% zz_new(i)=1/(2^17);
% end
% end
%
%
%
% subplot(6,1,1);
% plot(20*log10((xx)))
% subplot(6,1,2);
% plot(20*log10((xx_new)),'red')
%
%
%
%
% subplot(6,1,3);
% semilogy(yy)
%
% subplot(6,1,4);
% semilogy(yy_new,'red')
%
% subplot(6,1,5);
% semilogy(zz)
% subplot(6,1,6);
% semilogy(zz_new,'red')
%
%
% figure
% semilogy(yy_new)
% hold on
%
% semilogy(zz_new,'red')
end