%CIC濾波器設計
R_CIC = 40; % Decimation factor--CIC濾波器抽取因子
M    = 1;  % Differential delay--微分延遲，可為任何正整數，一般限制為1或2,SDR書上為M=1的情況
Nsecs= 5;  % Number of sections--CIC濾波器級數，級數越高，旁瓣抑制性能越好
IWL  = 14; % Input word length--輸入字長
OWL  = 18; % Output word length--輸出字長

% If the output wordlength is specified when creating a CIC filter then the
% "FilterInternals" property is set to "MinWordLengths" automatically.
% Therefore, the minimum word sizes are used between each section.
hcic1 = mfilt.cicinterp(3,M,Nsecs,IWL,IWL);
hcic2 = mfilt.cicdecim(R_CIC,M,Nsecs,IWL,OWL);

% hcic = mfilt.cicdecim(R_CIC,M,Nsecs,IWL,OWL);
% info(hcic)
% 
Fs_in =57.6e6;
% hcic=mfilt.cascade(hcic1,hcic2);
% h = fvtool(hcic,'Fs',Fs_in);
% set(gcf, 'Color', 'blue');

%增益歸一化，將幅頻響應最大值變為0dB
% hgain = dfilt.scalar(1/gain(hcic)); % Define gain
hgain1=dfilt.scalar(1/gain(hcic1));
hcicnorm1=cascade(hgain1,hcic1);
hgain2=dfilt.scalar(1/gain(hcic2));
hcicnorm2=cascade(hgain2,hcic2);
hcicnorm=cascade(hcicnorm1,hcicnorm2);
% hgain=hgain1*hgain2;
% hcicnorm = cascade(hgain,hcic);

% Replace the CIC in FVTool with a normalized CIC.
h=fvtool(hcicnorm,'Fs',Fs_in);
set(gcf,'Color','blue');
% setfilter(h,hcicnorm,'Fs',Fs_in);
% axis([0 0.2 -0.8 0]);

%CIC補償濾波器設計
%Compensation FIR Decimator
% Filter specifications
Fs_hcfir = Fs_in/R_CIC*3; % Sampling frequency 57.6MHz/20--CIC補償濾波器輸入采樣速率
Apass  = 0.01;     % dB  --通帶截止頻率處衰減值
Astop  = 60;       % dB  --阻帶截止頻率處衰減值
Aslope = 60;       % 60 dB slope over half the Nyquist range
Fpass  = 80e3;    % Hz passband-edge frequency  --通帶截止頻率
Fstop  = 640e3;    % Hz stopband-edge frequency  --阻帶截止頻率，值越大，所要求的濾波器階數越少

% Design decimation filter. D and Nsecs have been defined above as the
% differential delay and number of sections, respectively.
% 補償濾波器設計，抽取因子=2
d = fdesign.decimator(2,'ciccomp',M,Nsecs,Fpass,Fstop,Apass,Astop,Fs_hcfir);
hcfir = design(d,'equiripple',...
               'StopbandShape', 'linear',...
               'StopbandDecay', Aslope);

% Now we have to define the fixed-point attributes of our multirate filter.
% By default, the fixed-point attributes of the accumulator and multipliers
% are set to ensure that full precision arithmetic is used, i.e. no
% quantization takes place.
set(hcfir,...
    'Arithmetic',      'fixed',...
    'CoeffWordLength',  16,...
    'InputWordLength',  16);
info(hcfir)

hcas1 = cascade(hcicnorm,hcfir);
set(h,'Filters', [hcicnorm,hcfir,hcas1],'Fs',[Fs_in,Fs_in/R_CIC,Fs_in]);
axis([0 0.12 -0.8 0.8]);  %查看補償后的效果
legend(h,'hcic','hcfir','cascade');

%最后一級FIR濾波器設計，抽取因子=2（通過fvtool可查看所設計濾波器的系數）
%Third Stage FIR Decimator
N = 225;       % 63 taps  ，FIR濾波器階數
Fs_fir = Fs_hcfir/2;  % FIR濾波器輸入采樣速率
Fpass = 80e3;         % 通帶截止頻率
Fstop = 100e3;        % 阻帶截止頻率

% FIR濾波器設計，抽取因子=2
d = fdesign.decimator(1,'lowpass','N,Fp,Fst',N,Fpass,Fstop,Fs_fir);
hpfir = design(d,'equiripple','Wpass',2);  % Give more weight to passband
set(hpfir,...
    'Arithmetic',      'fixed',...
    'CoeffWordLength',  16,...
    'InputWordLength',  20,...
    'InputFracLength', -12);
set(hpfir,...
    'FilterInternals', 'specifyPrecision',...
    'outputWordLength', 20,...
    'outputFracLength',-12,...
    'RoundMode',       'round',... % = nearest in SL
    'OverflowMode',    'Saturate');
hcasnorm = cascade(hcicnorm,hcfir,hpfir);
set(h,'Filters',hcasnorm,'Fs',Fs_in,'NumberofPoints',8192*3);
axis([0 1 -200 10]);  % Zoom-in
axis([0 0.1 -0.2 0.2]);

%取出濾波器系數
hpfir_coef=get(hpfir,'Numerator');
hpfir_N=16;
max_coef=max(abs(hpfir_coef));
hpfir_coef=hpfir_coef/max_coef;
hpfir_coef=round(hpfir_coef* 2^(hpfir_N-1)-1);
input_width=16;
output_width=input_width+ceil(log2(sum(abs(hpfir_coef))))


%取出濾波器系數
hcfir_coef=get(hcfir,'Numerator');
hcfir_N=16;
max_coef=max(abs(hcfir_coef));
hcfir_coef=hcfir_coef/max_coef;
hcfir_coef=round(hcfir_coef* 2^(hcfir_N-1)-1);
input_width=16;
output_width=input_width+ceil(log2(sum(abs(hcfir_coef))))