% M09-bpfilt.diary % Simple bandpass filter example % 2007-10-18 % w_c = 0.4 pi, bandwidth B = 0.1 pi wc = 0.4*pi; B = 0.1*pi; % hence... beta = cos(wc) %beta = % 0.3090 % quadratic in alpha for B roots([cos(B) -2 cos(B)]) %ans = % 1.3764 % 0.7265 % |alpha| must be < 1 for stable system alpha = ans(2) %alpha = % 0.7265 % Hence numerator and denominator coefficients b = (1-alpha)/2*[1 0 -1]; a = [1, -beta*(1+alpha), alpha]; % Plot the response freqz(b,a) % Adjust the axes for more detail in passband axis([0 1 -20 0]) % Peak at 0.4pi, -3 dB bandwidth about 0.1pi % Apply the filter to a signal [x,sr]=wavread('mpgr1_sx419.wav'); soundsc(x,sr); % Apply the filter y = filter(b,a,x); soundsc(y,sr); % Bandpassed sound! % Visualize it: take spectrum of first 1 sec of sound subplot(311) X = fft(x,16384); plot([0:8192]/8192, 20*log10(abs(X(1:8193)))); % plot in dB, 1st half % Make all the dB axes the same; 80 dB (-60..20) is plenty of range axis([0 1 -60 20]) grid % compare to spectrum of output subplot(312) Y = fft(y,16384); plot([0:8192]/8192, 20*log10(abs(Y(1:8193)))); axis([0 1 -60 20]) grid % .. and compare to spectrum of filter [H,W] = freqz(b,a); subplot(313) plot(W/pi, 20*log10(abs(H))); axis([0 1 -60 20]) grid % Output is clearly input spectrum 'shifted' by magnitude response % (multiplication becomes an additive shift with a db (log) axis) % .. or try spectrograms to see whole sound subplot(211) specgram(x,512,sr); subplot(212) specgram(y,512,sr);