%實驗4:用窗函數法設計FIR數字濾波器
clc;
N=input('輸入窗函數長度(默認值15,選擇范圍[10-64]):\nN=','s');

N=str2num(N);
if N<1
    N=15;
else
    if N<10 | N>64
        N=15;
    end
end
w=input('輸入截止頻率(默認值pi/4,選擇范圍[0.5-pi]):\nwc=','s');

if w<1
    w=pi/4;
else
    if w<0.5 | w>pi
        w=pi/4;
    end
end
b=1;
close all;

i=0;
while(b)
    n=[0:(N-1)];

%產生理想低通濾波器單位脈沖響應hd(n)
alpha=(N-1)/2;
n=[0:(N-1)];
Q=n-alpha+eps;
hd=sin(w*Q)./(pi*Q);

    k=input('選擇窗口類型(默認1):\n[1]矩形窗\n[2]Hamming窗\n[3]Hanning窗\n[4]Blackman\n選擇:','s');
    switch(k)
        case{'1','2','3','4'}
            k=str2num(k);
        otherwise
            k=1;
    end
    switch(k)
        case 1
            B=boxcar(N);
            string='Boxcar';
        case 2
            B=hamming(N);
            string='Hamming';
        case 3
            B=hanning(N);
            string='Hanning';
        case 4
            B=blackman(N);
            string='Blackman';
    end
    string=[string,' N=',num2str(N)];
    
    h=hd.*(B)';
    [H,m]=freqz(h,[1],1024,'whole');
    mag=abs(H);
    db=20*log10((mag+eps)/max(mag));
    pha=angle(H);
    i=i+1;
    
    figure(i);
    subplot(2,2,1);
    n=0:N-1;
    stem(n,h,'.');
    axis([0 N-1 -0.1 0.3]);
    hold on;
    n=0:N-1;
    x=zeros(N);
    plot(n,x,'-');
    xlabel('n');
    ylabel('h(n)');
    title('實際低通濾波器的h(n)');
    text((0.3*N),0.27,string);
    hold off;
    
    subplot(2,2,2);
    plot(m/pi,db);
    axis([0 1 -100 0]);
    xlabel('w/pi');
    ylabel('20lg|H(e^j^w)|');
    title('衰減特性/dB');
    grid;
    
    subplot(2,2,3);
    plot(m,pha);
    hold on;
    n=0:7;
    x=zeros(8);
    plot(n,x,'-');
    axis([0 3.15 -4 4]);
    xlabel('w');
    ylabel('相位/rad');
    title('相頻特性');
    
    subplot(2,2,4);
    plot(m,mag);
    axis([0 3.15 0 1.5]);
    xlabel('w');
    ylabel('|H(e^j^w)|');
    title('幅頻特性');
    text(0.9,1.2,string);
    
    b=input('繼續實驗嗎?(默認1)\n[1]繼續\n[0]退出\n選擇:','s');
    switch(b)
        case{'0','1'}
            b=str2num(b);
        otherwise
            b=1;
    end
    if b==1
        clc;
        N=input('輸入窗函數長度(默認值15,選擇范圍[10-64]):\nN=','s');
        N=str2num(N);
        if length(N)<1
            N=15;
        else
            if N<10 | N>64
                N=15;
            end
        end
        w=input('輸入截止頻率(默認值pi/4,選擇范圍[0.5-pi]):\nwc=','s');
        w=str2num(w);
        if length(w)<1
            w=pi/4;
        else
            if w<0.5 | w>pi
                w=pi/4;
            end
        end
    end
end