% filtdemo.m
%
% The program includes various low pass, bandpass
% and notch filters.

% Consider 256 Hz carrier wave
% sampled at 8192 Hz, i.e. about 32 samples
% per cycle.
% Select parameters by editing this file.

% =======================================
% MATLAB HELP
% to learn more about Matlab, type 'intro'
% (without quotes) at the command prompt
% also type 'demo' for more complete intro
% for help on any function, type 'help functionname'
% for help on syntax:
% type 'help punct' 'help ops'  'help arith'
% for web-based help, type 'helpdesk'
% keyword search, type 'lookfor keyword'
% to retain variables in workspace on error
% type 'dbstop if error' before running 


% =======================================
% set universal parameters
% =======================================
fs=8192; % sampling rate, suggest 8192 (power of 2)
T=0.03125; %wavefile time in sec, suggest 1/32
fc=256; % carrier frequency, suggest 256 (middle C)


% =======================================
% input signal x
% =======================================

% want to write sampled cos (2 pi fc t)
% but in sampled systems t = n*ts = n/fs
% where n is an integer index and fs=1/ts is the 
% sampling rate 
% fs*T is the total number of samples in the file
% (1:fs*T)/(fs*T) is vector of fs*T numbers
% ranging from 0+ to 1
% index n is represented by vector of increasing
% integers (1:fs*T) = [ 1 2 3 ... fs*T]
% for more info, type 'help colon'

% cosine wave
% 2 pi fc t is written as below
twopi_fc_t=(1:fs*T)*2*pi*fc/fs; 
a=1;
phi=0;
x = a * cos(twopi_fc_t + phi);

% square wave
% x(x>0)=1;
% x(x<0)=-1;

% impulse signal
% x = [1 zeros(1,fs*T-1)];

% =======================================
% FILTERS H(z)
% =======================================
% Here we choose filter coefficient vectors for
% four example filters:
% two LPFs, one BPF and one notch filter.
% There are two vectors for each filter:
% b (numerator) and a (denominator) of H(z)
% then use Matlab filter function
% with input vector x to get output
% y = filter(b,a,x) 
% =======================================
% filter gains and bandwidths may need
% adjustment
% check these out to be sure they are correct

% low pass filter, 1 pole 1 zero, cutoff fmax
% gain set to be unity at zero frequency (z=1)
% H(z)=gain*(z+1)/(z-p)
% pole p on real axis, zero at -1
% magnitude pmag=|p| near but < 1
% choose pole pmag for about 3dB down at fmax
fmax = fs/32;
pmag=1-(1/2)*(2*pi*fmax/fs);
gain=(1-pmag)/2; % to get unity gain at z=1
bfmax=gain*[1 1];% numerator (zeros)
afmax=[1 -pmag];% denominator (poles)

% low pass filter, 2 pole 2 zero, cutoff fmax
% gain set to be unity at zero frequency (z=1)
% H(z)=gain*(z+1)(z+1)/(z-p)(z-p)
%     =gain*(z^2+2*z+1)/(z^2-2*p*z+p^2)
% double pole p on real axis, double zero at -1
% magnitude pmag=|p| near but < 1
% choose pole pmag for about 3dB down at fmax
% 2-sided 3dB BW(Hz) = (1-pmag)/(pi*sqrt(pmag))
pmag2=1-(2*pi*fmax/fs);
gain2=(1-2*pmag2+pmag2^2)/4; % to get unity gain at z=1
bfmax2=gain2*[1 2 1];% numerator (zeros)
afmax2=[1 -2*pmag2 pmag2^2];% denominator (poles)

% 2-pole 2-zero bandpass filter at fc
% zeros at +1, -1, poles at p=|p|e^+-j(theta)
% where theta = 2pi fc/fs
% for peak response at fc
% H(z) = (z+1)(z-1)/(z-p)(z-p*)
%      = (z^2-1)/(z^2-2*|p|*cos(theta)*z+|p|^2)
theta3=2*pi*(fc)/fs;
p3=0.95; % sets bandwidth
% 2-sided 3dB BW(Hz) = (1-p1)/(pi*sqrt(p1))
%gain= % should scale for unity output at peak
afc=[1 -2*p3*cos(theta3) p3*p3]; % denominator (poles)
bfc=[1 0 -1]; % numerator (zeros)
  
% 2-pole 2-zero notch filter at fc
% zeros at z1=e^+-j(theta),
% poles at p=|p|e^+-j(theta)
% where theta =  2*pi fc/fs
% for peak response at fc
% H(z) = (z-z1)(z-z1*)/(z-p)(z-p*)
% complete yourself

% plot frequency response of one of these filters
% two plots: linear scale and log scale
% first choose values of a and b
a=afmax;b=bfmax;
zz=exp(j*(0:fs/2-1)*2*pi/fs); % z values along unit circle
%H=(b(1)*zz.^2+b(2)*zz+b(3))./(a(1)*zz.^2+a(2)*zz+a(3));
H=(b(1)*zz+b(2))./(a(1)*zz+a(2));
H(H==0)=10^(-6); % H(1)=0 is really zero, want to avoid log(0)
figure(1);
subplot(2,1,1);
plot((1:fs/2)-1,abs(H));
axis([0 4*fc 0 max(abs(H))]);
ylabel('amplitude of H');
subplot(2,1,2);
plot((1:fs/2),20*log10(abs(H/max(H))));
axis([0 fs/2 -40 0]);
title('|H(z)| 1-pole, 1-zero LPF');
xlabel('frequency (Hz)');
ylabel('20*log|H|');

% ===============================================
% output y = filter(b,a,x)
% ===============================================


y = filter(b,a,x);

% plot input x and output y
% observe amplitude and phase shift of y compared to x

figure(2);
plot(x,'b-'); hold on;
plot(y,'r-');
hold off;

