%
% FunSP
% 1º semester 2022/2023
% Week 10
%
% author: AJF
%


% read WAV file
inpfile='soundfile.wav'; NBITS=16;
[x,FS,NBITS]=wavread(inpfile); % or [x,FS]=audioread(inpfile); 
sound(x,FS);

% design filter
order=126 ;
f=[COMPLETEHERE];
m=[COMPLETEHERE] ;
p=[COMPLETEHERE] ;

% uncomment just one of the following lines
h=firpm(order,f,m,p); % length=order+1
% h=fir1(order, COMPLETEHERE, hann(order+1));

% check the filter characteristics
[H W]=freqz(h,1.0, 512);
figure(1);
plot([0:order],h)
xlabel('Time (samples)');
ylabel('Amplitude');
title('Impulse response');

figure(2);
% note the scale change [0, PI] -> [0, Nyquist]
xaxis=W/pi*FS/2; % from 0 up to the Nyquist frequency (Hertz) 
plot(xaxis,20*log10(abs(H)));
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Frequency response magnitude');
pause;

% filter the input signal
y=filter(h,1.0,x) ;

% limit the dynamic range of the output signal and play sound
y=y/(1.2*max(abs(real(y))));
sound(y,FS);

% store the output signal on a WAV file
outfile=('output.wav');
wavwrite(y,FS,NBITS,outfile); % or audiowrite(outfile,y,FS,'BitsPerSample',NBITS);