%% HW6 Program Assignment
% By: 105060012 張育菘
%
%% P1
% Consider again the inverse DFT given in (8.2).
% x\left\lbrack n\right\rbrack =\frac{1}{N}\sum_{k=0}^{N-1} X\left\lbrack k\right\rbrack W_N^{-\textrm{kn}} ,\;n=0,1,\ldotp \ldotp \ldotp ,N-1 (8.2)
%
% (a) Replace k by {\left\langle -k\right\rangle }_N in (8.2) and show that the resulting summation is a DFT expression, that is, \mathrm{IDFT}\left\lbrace X\left\lbrack k\right\rbrack \right\rbrace =\frac{1}{N}\mathrm{DFT}\left\lbrace X\left\lbrack {\left\langle -k\right\rangle }_N \right\rbrack \right\rbrace.
% Ans.
% \textrm{IDFT}\left\lbrace X\left\lbrack k\right\rbrack \right\rbrace =x\left\lbrack n\right\rbrack =\frac{1}{N}\sum_{k=0}^{N-1} X\left\lbrack k\right\rbrack W_N^{-\textrm{kn}}
% \textrm{DFT}\left\lbrace X\left\lbrack {\left\langle -k\right\rangle }_N \right\rbrack \right\rbrace =\sum_{k=0}^{N-1} X\left\lbrack {\left\langle -k\right\rangle }_N \right\rbrack W_N^{\textrm{kn}} =\sum_{k=0}^{N-1} X\left\lbrack k\right\rbrack W_N^{-\textrm{kn}}
% \Longrightarrow \textrm{IDFT}\left\lbrace X\left\lbrack k\right\rbrack \right\rbrace =\frac{1}{N}\sum_{k=0}^{N-1} X\left\lbrack k\right\rbrack W_N^{-\textrm{kn}} =\frac{1}{N}\textrm{DFT}\left\lbrace X\left\lbrack {\left\langle -k\right\rangle }_N \right\rbrack \right\rbrace
%
% (b) Develop a MATLAB function x = IDFT(X,N) using the fft function that uses the above approach. Verify your function on signal x[n] = {1, 2, 3, 4, 5, 6, 7, 8}.

close all; clear;
fprintf('1(b)\n');

open IDFT.m;

x = [1 2 3 4 5 6 7 8]
X = fft(x); N = length(X);
x_verify = IDFT(X, N)

% (Comment)
% 由上面的結果顯示，從自己設計的 IDFT 把經由 fft 轉換得到的 X 轉回去 time domain 的 x_verify 與原本的 x 相同。
%
%% P2
% In this problem we will investigate differences in the speeds of DFT and FFT algorithms when stored twiddle factors are used.
% \left\lbrack \begin{array}{c}
% X\left\lbrack 0\right\rbrack \\
% X\left\lbrack 1\right\rbrack \\
% X\left\lbrack 2\right\rbrack \\
% X\left\lbrack 3\right\rbrack \\
% X\left\lbrack 4\right\rbrack \\
% X\left\lbrack 5\right\rbrack \\
% X\left\lbrack 6\right\rbrack \\
% X\left\lbrack 7\right\rbrack
% \end{array}\right\rbrack =\left\lbrack \begin{array}{cccccccc}
% 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\
% 1 & W_8  & W_8^2  & W_8^3  & W_8^4  & W_8^5  & W_8^6  & W_8^7 \\
% 1 & W_8^2  & W_8^4  & W_8^6  & 1 & W_8^2  & W_8^4  & W_8^6 \\
% 1 & W_8^3  & W_8^6  & W_8^{1\;}  & W_8^4  & W_8^7  & W_8^2  & W_8^5 \\
% 1 & W_8^4  & 1 & W_8^4  & 1 & W_8^4  & 1 & W_8^4 \\
% 1 & W_8^5  & W_8^2  & W_8^7  & W_8^4  & W_8  & W_8^6  & W_8^3 \\
% 1 & W_8^6  & W_8^4  & W_8^2  & 1 & W_8^6  & W_8^4  & W_8^2 \\
% 1 & W_8^7  & W_8^6  & W_8^5  & W_8^4  & W_8^3  & W_8^2  & W_8
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% x\left\lbrack 0\right\rbrack \\
% x\left\lbrack 1\right\rbrack \\
% x\left\lbrack 2\right\rbrack \\
% x\left\lbrack 3\right\rbrack \\
% x\left\lbrack 4\right\rbrack \\
% x\left\lbrack 5\right\rbrack \\
% x\left\lbrack 6\right\rbrack \\
% x\left\lbrack 7\right\rbrack
% \end{array}\right\rbrack  (8.8)
%
% (a) Write a function W = dft_matrix(N) that computes the DFT matrix W_N given in (8.8).
% Ans.
% W_N =e^{-j\frac{2\pi }{N}} , {\left(W_N \right)}_{\left(n+1\right)\times \left(m+1\right)} ={\left(e^{-j\frac{2\pi }{N}} \right)}^{nm} ,\;n,m=0,1,\ldotp \ldotp \ldotp ,N-1

close all; clear;
fprintf('2(a)\n');
open dft_matrix.m;

%
% (b) Write a function X = dftdirect_m(x,W) that modifies the dftdirect function using the matrix W from (a). Using the tic and toc functions compare computation times for the dftdirect and dftdirect_m function for N = 128, 256, 512, and 1024. For this purpose generate an N-point complex-valued signal as x = randn(1,N) + 1j*randn(1,N). (verify your code with fft first)

close all; clear;
fprintf('2(b)\n');

open dftdirect_m.m;

N_set = [128 256 512 1024];
for i = 1:1:length(N_set)
    fprintf('N = %d\n', N_set(i));
    x = randn(1,N_set(i)) + 1j*randn(1,N_set(i));
    
    fprintf('  by dftdirect_m\n');
    tic;
    W = dft_matrix(N_set(i));
    X1 = dftdirect_m(x, W)
    toc;
    
    fprintf('  by dftdirect\n');
    tic; 
    X2 = dftdirect(x)
    toc;
end

% Ans:
% 由output結果可得知，x 經過 W = dft_matrix(N) 和 X = dftdirect_m(x,W) 得到 X 的速度比直接做 DFT (X = dftdirect(x)) 快，這是因為 dftdirect 是把所有 W*x 的 term 算完後再做加總；而前者的方法是將 W 的 matrix 算完作為 input 至 dftdirect_m 得到 X，因此此方法會比較快。
%
% (c) Write a function X = fftrecur_m(x,W) that modifies the fftrecur function given on page 439 using the matrix W from (a). Using the tic and toc functions compare computation times for the fftrecur and fftrecur_m function for N = 128, 256, 512, and 1024. For this purpose generate an N-point complex valued signal as x = randn(1,N) + 1j*randn(1,N). (verify your code with fft first)
% Ans.
% X=W_8 x=W_4 \left\lbrack \begin{array}{cc}
% I & D_8 \\
% I & -D_8
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% x_E \\
% x_O
% \end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
% I & D_8 \\
% I & -D_8
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% W_4 x_E \\
% W_4 x_O
% \end{array}\right\rbrack ,
% W_4 x_E =W_2 \left\lbrack \begin{array}{cc}
% I & D_4 \\
% I & -D_4
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% x_{\textrm{EE}} \\
% x_{\textrm{EO}}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
% I & D_4 \\
% I & -D_4
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% W_2 x_{\textrm{EE}} \\
% W_2 x_{\textrm{EO}}
% \end{array}\right\rbrack ,
% W_2 x_{\textrm{EE}} =W_1 \left\lbrack \begin{array}{cc}
% I & D_2 \\
% I & -D_2
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% x_{\textrm{EEE}} \\
% x_{\textrm{EEO}}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
% I & D_2 \\
% I & -D_2
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% x_{\textrm{EEE}} \\
% x_{\textrm{EEO}}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
% 1 & 1\\
% 1 & -1
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% x_{\textrm{EEE}} \\
% x_{\textrm{EEO}}
% \end{array}\right\rbrack

close all; clear;
fprintf('2(c)\n');

open fftrecur_m.m;

N_set = [128 256 512 1024];
for i = 1:1:length(N_set)
    fprintf('N = %d\n', N_set(i));
    x = randn(1,N_set(i)) + 1j*randn(1,N_set(i));
    
    fprintf('  by fftrecur_m\n');
    tic;
    W = dft_matrix(N_set(i));
    X1 = fftrecur_m(x, W)
    toc;
    
    fprintf('  by fftrecur\n');
    tic; 
    X2 = fftrecur(x)'
    toc;
end

% 由output結果可得知，x 經過 fftrecur_m 得到 X 的速度比做 fftrecur 慢，這是因為 fftrecur 是直接算所需的 w；而前者的方法在每次迴圈都會把 W 的 matrix 作為 input 至下一層的 fftrecur_m ，等於是在每次迴圈都會傳一些不會用到的資料，因此此方法會比較慢。
%
%% P3
% Consider the flow graph in Figure 8.10 which implements a DIT-FFT algorithm with both input and output in natural order. Let the nodes at each stage be labeled as sm[k],0 ≤ m ≤ 3 with s0[k] = x[k] and s3[k] = X[k], 0 ≤ k ≤ 7.
%
%
% (a) Express s_m \left\lbrack k\right\rbrack in terms of s_{m-1} \left\lbrack k\right\rbrack for m = 1, 2, 3.
% Ans.
% s_1 \left\lbrack \mathrm{k1}\right\rbrack =x\left\lbrack \mathrm{k1}\right\rbrack +x\left\lbrack \mathrm{k1}+4\right\rbrack ,s_1 \left\lbrack \mathrm{k1}+4\right\rbrack =x\left\lbrack \mathrm{k1}\right\rbrack -x\left\lbrack \mathrm{k1}+4\right\rbrack ,\mathrm{k1}=0,1,2,3
% s_2 \left\lbrack \mathrm{k2}\right\rbrack =s_1 \left\lbrack \mathrm{k2}\right\rbrack +s_1 \left\lbrack \mathrm{k2}+2\right\rbrack ,s_2 \left\lbrack \mathrm{k2}+2\right\rbrack =s_1 \left\lbrack \mathrm{k2}+4\right\rbrack +W_8^2 s_1 \left\lbrack \mathrm{k2}+6\right\rbrack ,
% s_2 \left\lbrack \mathrm{k2}+4\right\rbrack =s_1 \left\lbrack \mathrm{k2}\right\rbrack -s_1 \left\lbrack \mathrm{k2}+2\right\rbrack ,s_2 \left\lbrack \mathrm{k2}+6\right\rbrack =s_1 \left\lbrack \mathrm{k2}+4\right\rbrack -W_8^2 s_1 \left\lbrack \mathrm{k2}+6\right\rbrack ,\;\mathrm{k2}=0,1
% s_3 \left\lbrack 0\right\rbrack =s_2 \left\lbrack 0\right\rbrack +s_2 \left\lbrack 1\right\rbrack ,s_3 \left\lbrack 1\right\rbrack =s_2 \left\lbrack 2\right\rbrack +W_8^1 s_2 \left\lbrack 3\right\rbrack ,s_3 \left\lbrack 2\right\rbrack =s_2 \left\lbrack 4\right\rbrack +W_8^2 s_2 \left\lbrack 5\right\rbrack ,s_3 \left\lbrack 3\right\rbrack =s_2 \left\lbrack 6\right\rbrack +W_8^3 s_2 \left\lbrack 7\right\rbrack
% s_3 \left\lbrack 4\right\rbrack =s_2 \left\lbrack 0\right\rbrack -s_2 \left\lbrack 1\right\rbrack ,s_3 \left\lbrack 5\right\rbrack =s_2 \left\lbrack 2\right\rbrack -W_8^1 s_2 \left\lbrack 3\right\rbrack ,s_3 \left\lbrack 6\right\rbrack =s_2 \left\lbrack 4\right\rbrack -W_8^2 s_2 \left\lbrack 5\right\rbrack ,s_3 \left\lbrack 7\right\rbrack =s_2 \left\lbrack 6\right\rbrack -W_8^3 s_2 \left\lbrack 7\right\rbrack
%
% (b) Write a MATLAB function X = fftalt8(x) that computes an 8-point DFT using the equations in part (a). Verify with sequence x[n] = {0,1,2,2,3,3,3,4}.

close all; clear;
fprintf('3(b)\n');

open fftalt8.m;

x = [0 1 2 2 3 3 3 4];
X_fft = fft(x,8)
X_fftalt8 = fftalt8(x)

%
% (c) Compare the coding complexity of the above function with that of MATLAB function fftditr2 shown in Figure 8.6, and comment on its usefulness.

close all; clear;
fprintf('3(c)\n');
x = [0 1 2 2 3 3 3 4];
X_fftditr2 = fftditr2(x)

% (Comment)
% 就 coding complexity 而言，fftalt8 顯得複雜，因為需要將前一級算完才能算下一級，而且架構相對 fftditr2 複雜，沒辦法像 fftditr2 用三層 for 迴圈就完成此 function，因次在 coding 的複雜度方面也就比 fftditr2 複雜。
%
%% P4
% Using the flow graph of Figure 8.13 and following the approach used in developing the fftditr2 function.
%
%
% (a) Develop a radix-2 DIF-FFT function X = fftdifr2(x) for power-of-2 length N.

close all; clear;
fprintf('4(a)\n');

open fftdifr2.m;

%
% (b) Verify your function for N = 2^v, where 2 ≤ ν ≤ 10. For this purpose generate an N-point complex-valued signal as x = randn(1,N) + 1j*randn(1,N).

fprintf('4(b)\n');

for v = 2:10
    N = 2^v;
    x = randn(1,N) + 1j*randn(1,N);
    fprintf('for N = 2^%d\n', v);
    X_fftdifr2 = fftdifr2(x)
    X_fft = fft(x)
end

% (Comment)
% 由上述結果我們發現 x 經過 fftdifr2 和 fft 得到的結果相同。
%
%% P5
% The filterfirdf implements the FIR direct form structure.
%
% (a) Develop a new MATLAB function y=filterfirlp(h,x) that implements the FIR linear-phase form given its impulse response in h. This function should first check if h is one of type-I through type-IV and then simulate the corresponding equations. If h does not correspond to one of the four types then the function should display an appropriate error message.
% Ams.
% For Type I: y\left\lbrack n\right\rbrack =\sum_{k=0}^{\frac{M}{2}-1} h\left\lbrack k\right\rbrack \left(x\left\lbrack n-k\right\rbrack +x\left\lbrack n-M+k\right\rbrack \right)+h\left\lbrack \frac{M}{2}\right\rbrack x\left\lbrack n-\frac{M}{2}\right\rbrack
% For Type II: y\left\lbrack n\right\rbrack =\sum_{k=0}^{\frac{M-1}{2}} h\left\lbrack k\right\rbrack \left(x\left\lbrack n-k\right\rbrack +x\left\lbrack n-M+k\right\rbrack \right)
% For Type III: y\left\lbrack n\right\rbrack =\sum_{k=0}^{\frac{M}{2}-1} h\left\lbrack k\right\rbrack \left(x\left\lbrack n-k\right\rbrack -x\left\lbrack n-M+k\right\rbrack \right)
% For Type IV: y\left\lbrack n\right\rbrack =\sum_{k=0}^{\frac{M-1}{2}} h\left\lbrack k\right\rbrack \left(x\left\lbrack n-k\right\rbrack -x\left\lbrack n-M+k\right\rbrack \right)

close all; clear;
fprintf('5(a)\n');

open filterfirlp.m;

%
% (b) Verify your function on each of the following FIR systems:
% h1[n] = {1,2,3,2,1}, h2[n] = {1,-2,3,3,-2,1}, h3[n] = {1,-5,0,5,-1}, h4[n] = {1,-3,-4,4,3,-1}, h5[n] = {1,2,3,-2,-1},
% For verification determine the first ten samples of the step responses using your function and compare them with those from the filter function.

fprintf('5(b)\n');

x = ones(1, 10);

h1 = [1 2 3 2 1];
[y1, type] = filterfirlp(h1, x)
y1_verify = filterfirdf(h1, x)

h2 = [1 -2 3 3 -2 1];
[y2, type] = filterfirlp(h2, x)
y2_verify = filterfirdf(h2, x)

h3 = [1 -5 0 5 -1];
[y3, type] = filterfirlp(h3, x)
y3_verify = filterfirdf(h3, x)

h4 = [1 -3 -4 4 3 -1];
[y4, type] = filterfirlp(h4, x)
y4_verify = filterfirdf(h4, x)

h5 = [1 2 3 -2 -1];
[y5, type] = filterfirlp(h5, x)
y5_verify = filterfirdf(h5, x)

%
%% P6
% Consider the IIR normal direct form II structure given in Figure 9.6 and implemented by (9.18) and (9.20).
%
%
% (a) Using the MATLAB function filterdf1 as a guide, develop a MATLAB function y=filterdf2(b,a,x) that implements the normal direct form II structure. Assume zero initial conditions.

close all; clear;
fprintf('6(a)\n');

open filterdf2.m;

%
% (b) Determine y[n], 0 ≤ n ≤ 500 using your function and filterdf1 function with following inputs:
% x\left\lbrack n\right\rbrack ={\left(\frac{1}{4}\right)}^n u\left\lbrack n\right\rbrack ,a=\left\lbrack \begin{array}{ccc}
% 1 & -1\ldotp 5 & 0\ldotp 5
% \end{array}\right\rbrack ,b=1
% Compare your results to verify that your filterdf2 function is correctly implemented.

fprintf('6(b)\n');

n = 0:1:500;
x = (1/4).^n; a = [1 -1.5 0.5]; b = [1];
y = filterdf2(b,a,x)'
y_verify = filterdf1(b,a,x)'

% (Comment)
% 由上面的結果顯示，從 filterdf1 和 filterdf2 得到的結果相同。
%
%% P7
% The following numerator and denominator arrays in MATLAB represent the system function of a discrete-time system in direct form:
% b = [1,-2.61,2.75,-1.36,0.27], a = [1,-1.05,0.91,-0.8,0.38].
% Determine and draw each of the following structures:

close all; clear;
fprintf('7\n');
num = [1 -2.61 2.75 -1.36 0.27]; 
den = [1 -1.05 0.91 -0.8 0.38];
figure; zplane(num, den);
[num,den] = eqtflength(num,den);
[zero,pole,k] = tf2zp(num,den)

% (Comment)
% H\left(z\right)=\frac{1-2\ldotp 61z^{-1} +2\ldotp 75z^{-2} -1\ldotp 36z^{-3} +0\ldotp 27z^{-4} }{1-1\ldotp 05z^{-1} +0\ldotp 91z^{-2} -0\ldotp 8z^{-3} +0\ldotp 38z^{-4} }=\frac{\left(1+0\ldotp 26z^{-1} +0\ldotp 72z^{-2} \right)\left(1-1\ldotp 31z^{-1} +0\ldotp 53z^{-2} \right)}{\left(1-1\ldotp 41z^{-1} +0\ldotp 63z^{-2} \right)\left(1-1\ldotp 2z^{-1} +0\ldotp 43z^{-2} \right)}
%
% (a) Cascade form with second-order sections in normal direct form I,
% Ans.
% y\left\lbrack n\right\rbrack =-\sum_{k=1}^N a_k y\left\lbrack n-k\right\rbrack +\sum_{k=0}^M b_k x\left\lbrack n-k\right\rbrack ,H\left(z\right)=\frac{Y\left(z\right)}{X\left(z\right)}=\frac{\sum_{k=0}^M b_k z^{-k} }{1+\sum_{k=1}^N a_k z^{-k} }
%
%
% (b) Cascade form with second-order sections in transposed direct form I,
% w\left\lbrack n\right\rbrack =-\sum_{k=1}^N a_k w\left\lbrack n-k\right\rbrack +x\left\lbrack n\right\rbrack ,y\left\lbrack n\right\rbrack =\sum_{k=0}^M b_k w\left\lbrack n-k\right\rbrack
%
%
% (c) Cascade form with second-order sections in normal direct form II,
% W\left(z\right)=\frac{1}{1+\sum_{k=1}^N a_k z^{-k} },Y\left(z\right)=\sum_{k=0}^M b_k z^{-k} W\left(z\right)
%
%
% (d) Cascade form with second-order sections in transposed direct form II.
% \begin{array}{l}
% H\left(z\right)=\frac{Y\left(z\right)}{X\left(z\right)}=\frac{b_0 +b_1 z^{-1} +b_2 z^{-2} }{1+a_1 z^{-1} +a_2 z^{-2} },\\
% Y\left(z\right)=z^{-1} V_1 \left(z\right)+b_0 X\left(z\right),V_1 \left(z\right)=z^{-1} V_2 \left(z\right)-a_1 Y\left(z\right)+b_1 X\left(z\right),V_2 \left(z\right)=b_2 X\left(z\right)-a_2 Y\left(z\right)
% \end{array}
%
%
%% P8
% The frequency-sampling form is developed using (9.50) which uses complex arithmetic.
% H\left(z\right)=\frac{1-z^{-N} }{N}\sum_{k=0}^{N-1} \frac{H\left\lbrack k\right\rbrack }{1-z^{-1} e^{j\frac{2\pi }{N}k} },H\left\lbrack k\right\rbrack =H\left(z\right)|_{z=e^{j\frac{2\pi }{N}k} } (9.50)
%
% (a) Develop a MATLAB function [G,sos]=firdf2fs(h) that determines frequency sampling form parameters given in (9.51) and (9.52) given the impulse response in h. The matrix sos (second order system) should contain second-order section coefficients in the form similar to the tf2sos function while G array should contain the respective gains of second-order sections. Incorporate the coefficients for the H[0] and H[N/2] terms in sos and G arrays.
% H\left(z\right)=\frac{1-z^{-N} }{N}\left\lbrace \frac{H\left\lbrack 0\right\rbrack }{1-z^{-1} }+\frac{H\left\lbrack \frac{N}{2}\right\rbrack }{1+z^{-1} }+\sum_{k=1}^K 2|H\left\lbrack k\right\rbrack |H_k \left(z\right)\right\rbrace (9.51)
% H_k \left(z\right)=\frac{\mathrm{cos}\left(\angle H\left\lbrack k\right\rbrack \right)-z^{-1} \mathrm{cos}\left(\angle H\left\lbrack k\right\rbrack -\frac{2\pi k}{N}\right)}{1-2\mathrm{cos}\left(\frac{2\pi k}{N}\right)z^{-1} +z^{-2} } (9.52)

close all; clear;
fprintf('8(a)\n');

open firdf2fs.m;

%
% (b) Verify your function by input h with sampled frequency response (9.53) and compare with the system function (9.54) (see example 9.6 in the textbook)
% H\left\lbrack k\right\rbrack =H\left(e^{j\frac{2\pi }{33}k} \right)=e^{-j\frac{2\pi }{33}k} \times \left\lbrace \begin{array}{cc}
% 1 & k=0,1,2,31,32\\
% 0\ldotp 5 & k=3,30\\
% 0 & \mathrm{otherwise}
% \end{array}\right. (9.53)
% H\left(z\right)=\frac{1-z^{-33} }{33}\left\lbrack \frac{1}{1-z^{-1} }+\frac{-1\ldotp 99+1\ldotp 99z^{-1} }{1-1\ldotp 964z^{-1} +z^{-2} }+\frac{1\ldotp 964-1\ldotp 964z^{-1} }{1-1\ldotp 857z^{-1} +z^{-2} }+\frac{-1\ldotp 96+1\ldotp 96z^{-1} }{1-1\ldotp 683z^{-1} +z^{-2} }\right\rbrack (9.54)

fprintf('8(b)\n');
N = 33;
H = zeros(1, 33);
for k = 1:1:33
    if((k==1)||(k==2)||(k==3))
        H(k) = exp((-1*sqrt(-1)*2*pi*(k-1))/N);
    elseif((k==32)||(k==33))
        H(k) = exp((-1*sqrt(-1)*2*pi*(k-1))/N);
    elseif((k==4)||(k==31))
        H(k) = 0.5*(exp((-1*sqrt(-1)*2*pi*(k-1))/N));
    else
        H(k) = 0;
    end
end
h = ifft(H);
[G, sos] = firdf2fs(h);
sos_after = sos;
for a = 1:1:length(G)
    for b = 1:1:3
        sos_after(a,b) = sos_after(a,b)*G(a);
    end
end
sos_after % sos multiplied by the corresponding G factor

% 由此結果我們可以發現，經由 firdf2fs 處理後的值與 (9.54) 相近，一樣只有前面 4 項有值。