% multiband script

T = 10;         % Observation window
fc = ???;       % Cutoff frequency
fsamp = ???;    % Sampling frequency

% Remember that 0 and fsamp/2 are implicit band edges
% and are never included in ftrue.
ftrue =    [ ??? ];
Hideal =   [ ??? ];
deltavec = [ ??? ];

figure(1)
plotfilterspecs(ftrue,Hideal,deltavec,fsamp)

% Get Parks-McClellan Filter
[twotau,fnorm,Hout,w] = firpmord(ftrue,Hideal,deltavec,fsamp);
fprintf('Filter order = %g\n',twotau)
h = firpm(twotau,fnorm,Hout,w);

% Plot actual frequency response of filter on top of specs
[H,f] = dtftfft(h,0,-8000);
k = find(f>=0);
hold on
plot(fsamp*f(k),abs(H(k)),'k')
hold off

figure(2) % Plot transform of square wave
N = floor(T*fsamp);
t = [0:N]/fsamp;
x = sqwave(t,5);
X = dtftfft(x,0,-8000)/fsamp;
subplot(2,1,1); plot(fsamp*f(k),abs(X(k)),'m'); grid on

% Filter square wave and plot transform of output
y = filter(h,1,x);
Y = dtftfft(y,0,-8000)/fsamp;
subplot(2,1,2); plot(fsamp*f(k),abs(Y(k))); grid on

figure(3) % Plot input and output waveforms and ideal output
subplot(2,1,1); plot(t,x); grid on
z = ??? % Compute true sum of harmonics 1 and 5 of the square wave
subplot(2,1,2); plot(t,y,t,z,'r'); grid on