%Filter Designer (FD.m)

choice = input('Choose one of the following or type 999 to quit.\n 1. FIR Notch \n 2. IIR Notch \n 3. All-Pole Resonator \n 4. Pole-Zero Resonator \n 5. Oscillator \n ');

if(choice == 1)

    %FIR notch filter 

    FS = input('Sampling Frequency: ');
    F0 = input('Notch Frequency: ');
    if (F0 > FS/2)
        warnstr =['WARNING! The Sampling frequency is less than twice the notch frequency. Response will be inaccurate.']
    end
    FG = input('Frequency to Set Gain: ');
    if (FG == F0)
        warnstr =['WARNING! User should not set the gain of the notch frequency. Response will be inaccurate.']
    end
    fgstr = ['Gain (dB) at ' num2str(FG) 'Hz: '];
    FGAINdB = input(fgstr);
    FGAIN = 10^(FGAINdB/20);
    
    w0 = 2*pi*F0/FS;
    wg = 2*pi*FG/FS;
    
    GH = abs(1 - 2*cos(w0)*exp(-j*wg) + exp(-j*2*wg));
    G = FGAIN/GH;
    
    h = [1 -2*cos(w0) 1];
    
    freqz(G*h,1,[], FS);
    
    titlestr = ['FIR Notch Filter for ' num2str(F0) 'Hz'];
    title(titlestr)
    
    LCCDE = ['y(n) = ' num2str(G, 4) '[' num2str(h(1), 4) 'x(n) + ' num2str(h(2), 4) 'x(n-1) + ' num2str(h(3), 4) 'x(n-2)]'] 

elseif(choice == 2)

    %IIR Notch

    FS = input('Sampling Frequency: ');
    F0 = input('Notch Frequency: ');
    if (F0 > FS/2)
        warnstr =['WARNING! The Sampling frequency is less than twice the notch frequency. Response will be inaccurate.']
    end
    r = input('r: ');
    if (r >= 1 || r < 0)
        warnstr =['WARNING! r should be a positive value < 1']
    end
    FG = input('Frequency to Set Gain: ');
    if (FG == F0)
        warnstr =['WARNING! User should not set the gain of the notch frequency. Response will be inaccurate.']
    end
    fgstr = ['Gain (dB) at ' num2str(FG) 'Hz: '];
    FGAINdB = input(fgstr);
    FGAIN = 10^(FGAINdB/20);
    
    w0 = 2*pi*F0/FS;
    wg = 2*pi*FG/FS;
    
    bg = abs(1 - 2*cos(w0)*exp(-j*wg) + exp(-j*2*wg));
    ag = abs(1 - 2*cos(w0)*exp(-j*wg) + (r^2)*exp(-j*2*wg));
    GH = bg/ag;
    G = FGAIN/GH;
    
    b = [1 -2*cos(w0) 1];
    a = [1 -2*cos(w0) (r^2)];
    
    freqz(G*b,a,1024, FS);
    titlestr = ['IIR Notch Filter for ' num2str(F0) 'Hz and r = ' num2str(r)];
    title(titlestr)
    
    LCCDE = ['y(n) + ' num2str(b(2), 4) 'y(n-1) + ' num2str((r^2), 4) 'y(n-2) = ' num2str(G, 4) '[x(n) + ' num2str(b(1), 4) 'x(n-1) + x(n-2)]'] 

elseif(choice == 3)

    %All-Pole Resonator

    FS = input('Sampling Frequency: ');
    F0 = input('Resonating Frequency: ');
    if (F0 > FS/2)
        warnstr =['WARNING! The Sampling frequency is less than twice the resonating frequency. Response will be inaccurate.']
    end
    r = input('r: ');
    if (r >= 1 || r < 0)
        warnstr =['WARNING! r should be a positive value < 1']
    end
    FG = input('Frequency to Set Gain: ');
    fgstr = ['Gain (dB) at ' num2str(FG) 'Hz: '];
    FGAINdB = input(fgstr);
    FGAIN = 10^(FGAINdB/20);
    
    wr = 2*pi*F0/FS;
    if (r >= 0.9)
        w0 = wr;
    else
        w0 = acos((2*r*cos(wr))/(1 + r^2));
    end
    
    wg = 2*pi*FG/FS;
    
    ag = abs(1 - 2*r*cos(w0)*exp(-j*wg) + (r^2)*exp(-j*2*wg));
    GH = 1/ag;
    G = FGAIN/GH;
    
    a = [1 -2*r*cos(w0) (r^2)];

    freqz(G,a,[], FS);
    titlestr = ['All-Pole Resonator for ' num2str(F0) 'Hz and r = ' num2str(r)];
    title(titlestr)

    LCCDE = ['y(n) + ' num2str(a(2), 4) 'y(n-1) + ' num2str((r^2), 4) 'y(n-2) = ' num2str(G, 4) 'x(n)'] 
    
elseif(choice == 4)

    %Pole-Zero Resonator

    FS = input('Sampling Frequency: ');
    F0 = input('Resonating Frequency: ');
    if (F0 > FS/2)
        warnstr =['WARNING! The Sampling frequency is less than twice the resonating frequency. Response will be inaccurate.']
    end
    r = input('r: ');
    if (r >= 1 || r < 0)
        warnstr =['WARNING! r should be a positive value < 1']
    end
    FG = input('Frequency to Set Gain: ');
    if (FG == 0)
        warnstr =['WARNING! User should not set the DC (0Hz) gain. Response will be inaccurate.']
    end
    fgstr = ['Gain (dB) at ' num2str(FG) 'Hz: '];
    FGAINdB = input(fgstr);
    FGAIN = 10^(FGAINdB/20);
    
    wr = 2*pi*F0/FS;
    if (r >= 0.9)
        w0 = wr;
    else
        w0 = acos((2*r*cos(wr))/(1 + r^2));
    end
    
    wg = 2*pi*FG/FS;
    
    bg = abs(1 - exp(-j*2*wg));
    ag = abs(1 - 2*r*cos(w0)*exp(-j*wg) + (r^2)*exp(-j*2*wg));
    GH = bg/ag;
    G = FGAIN/GH;
    
    b = [1 0 -1];
    a = [1 -2*r*cos(w0) (r^2)];
    
    freqz(G*b,a,[], FS);
    titlestr = ['Pole-Zero Resonator for ' num2str(F0) 'Hz and r = ' num2str(r)];
    title(titlestr)
    
    LCCDE = ['y(n) + ' num2str(a(2), 4) 'y(n-1) + ' num2str((r^2), 4) 'y(n-2) = ' num2str(G, 4) '[x(n) - x(n-2)]'] 

elseif (choice == 5)
%Oscillator

    FS = input('Sampling Frequency: ');
    F0 = input('Oscillating Frequency: ');
    if (F0 > FS/2)
        warnstr =['WARNING! The Sampling frequency is less than twice the oscillating frequency. Waveform will be inaccurate.']
    end
    
    w0 = 2*pi*F0/FS;
    
    b = [0 sin(w0)];
    a = [1 -2*cos(w0) 1];
   
    LCCDE = ['y(n) + ' num2str(a(2), 4) 'y(n-1) + y(n-2) = ' num2str(b(2), 4) 'x(n-1)'] 
    
    x(1) = 1;
    x(2) = 0;
    y(1) = 0;
    y(2) = sin(w0)*x(1);
    for n = 3:FS;
        x(n) = 0;
        y(n) = 2*cos(w0)*y(n-1) - y(n-2) + sin(w0)*x(n-1);
    end

    t0 = 1:FS;
    t = t0/FS;
    plot(t, y)
    titlestr = [num2str(F0) 'Hz Oscillator'];
    title(titlestr)
    xlabel('Time (sec)')
    ylabel('Amplitude')

elseif (choice == 999)
    %do nothing and quit
else
    
    'Invalid Input.  Please enter a valid integer choice, 1-5 or type 999 to quit.'
    FD()

end
    
if (choice ~= 999)
    clear all
    FD()
end

clear all