EECE 301 Signals & Systems Prof. Mark FowlerEECE 301 Signals & SystemsProf. Mark FowlerDiscussion #12a• DT Filter ApplicationThese notes explore the use of DT filters to remove an interference (in this case a single tone) from an audio signal. Imagine thatyou are in the your home recording studio and have just recordedwhat you feel is a “perfect take” of a guitar solo for a song you are recording, but you discover that someone had turned on some nearby electronic device that caused electromagnetic radiation that was picked up somewhere in the audio electronics and was recorded on top of the guitar solo. Rather than try to recreatethis “perfect take” you decide that maybe you can design a DT filter to remove it. We will explore two different cases: i. a high-pitched tone that lies above the significant portion of the guitar signal’s spectrum, and ii. a mid-pitched tone that lies in the middle of the guitar signal’s spectrum.Amp & AA FilterADCStore in MemoryInterferenceDT Filter(Implemented in S/W & run by CPU)Store in MemoryDACAmpReal Set-UpI. Signal Access and Exploration1. Use MATLAB’s wavread command to load the guitar1.wav file2. Listen to the guitar signal using MATLAB’s sound command.3. Plot the first second or so of the signal in the time domain to see what the signal looks like.4. Look at the guitar signal in the frequency domain by computing and plotting (in dB) the DFT of various 16384-pt blocks of the guitar signal. Verify that the significant portion of the guitar signal’s spectrum lies below 5 kHz.%%%%%% I. Signal Access and Exploration %%%%%%%%%%[x,Fs]=wavread('guitar1.wav');x=x.'; % convert into row vectorsound(x,Fs);t=(0:49999)*(1/Fs);plot(t,x(1:50000))X1=fftshift(fft(x(20000+(1:16384)),65536));X2=fftshift(fft(x(40000+(1:16384)),65536));X3=fftshift(fft(x(60000+(1:16384)),65536));X4=fftshift(fft(x(80000+(1:16384)),65536));f=(-32768:32767)*Fs/65536;Figure;subplot(2,2,1); plot(f/1e3,20*log10(abs(X1)));subplot(2,2,2); plot(f/1e3,20*log10(abs(X2)));subplot(2,2,3); plot(f/1e3,20*log10(abs(X3)));subplot(2,2,4); plot(f/1e3,20*log10(abs(X4)));II. Adding A High Frequency Interference1. Create a sinusoid whose frequency is 10kHz that is sampled at the same rate as the guitar signal and has the same length. The amplitude of this sinusoid should be 1.2. Add this signal to the guitar signal to create the simulated recorded signal that has the interference (call this signal x_10 to indicate that it has an interference at 10 kHz).3. Listen to the guitar signal with interference using MATLAB’s sound command.4. Plot the first second or so of the signal with interference and the signal without interference.5. Compute the DFTs of the signal that has interference. Verify that the interference is outside the significant portion of the guitar spectrum.%%%%%% II. Adding A High Frequency Interference %%%%%%%omega=2*pi*(10000/Fs); % convert 10 kHz into DT frequencyN=length(x);n=0:(N-1);x_10=x+cos(omega*n);sound(x_10,Fs);figuret=(0:49999)*(1/Fs);plot(t,x_10(1:50000),'r',t,x(1:50000))X_10_1=fftshift(fft(x_10(20000+(1:16384)),65536));figure;subplot(2,1,1); plot(f/1e3,20*log10(abs(X_10_1)));subplot(2,1,2); plot(f/1e3,20*log10(abs(X1)));III. Filter DesignMATLAB contains some easy to use routines for designing FIR filters – FIR (finite-impulse response) filters don’t use any output feedback – therefore they don’t really have any poles and they will always be stable. They are the most widely used type of DT filter in practice.A simple FIR filter: ]2[]1[][][313131−+−+= nxnxnxnyA more general FIR filter: N = “Order of Filter”∑=−=Niiinxbny0][][Such filters are quite easy to design using software-based tools. We’ll use the MATLAB FIR design routines called remezord.m and remez.mThe command remezord will give an estimate of the FIR filter order needed to achieve given specifications. The routine remez.m will then give the required design.Ω|H(Ω)|1 + δp1 – δpδsΩpΩssppHHδ≤Ω≤δ+≤Ω≤δ−)(01)(1dB )(log20n Attenuatio StopbanddB )1(log20 Ripple Passband1010spδ−=δ+=Lowpass Filter SpecificationπHere is how we state the filter specifications:Ω|H(Ω)|1 + δp1 – δpδsΩpΩsHighpass Filter SpecificationπΩ|H(Ω)|1 + δp1 – δpδs1Ωp2Ωs1Bandpass Filter SpecificationπΩs2Ωp1δs2Ω|H(Ω)|1 + δp11 – δp1δsΩp2Ωs1Bandstop Filter SpecificationπΩs2Ωp11 + δp21 – δp21. Use the “remezord” and “remez” commands to design lowpass filter to get:− 60 dB of attenuationin the stopband for the undesired signal − 1 dB of passband ripple− passband edge at 7kHz− stopband edge at 9 kHz. Look at DFTs to see why 60 dB of attenuation is a reasonable choice.Use the MATLAB variable b for the vector that holds the FIR filter coefficients2. Plot the filter’s impulse response.− For an FIR filter it is easy to show that the impulse response is nothing more than the bicoefficients in its difference equation:3. Compute and plot the filter’s frequency response.4. Make a pole-zero plot for the filter’s transfer function:∑=−=Niiinxbny0][][NNNNNNNzbzbzbzbzbzbzbbzH++++=++++=−−−−−LL2211022110)(III. Lowpass Filter Design Steps%%%%%% III. Lowpass Filter Design %%%%%%%%%% Lowpass Filter Design Specifications:% · Passband cutoff frequency = 7 kHz% · Stopband cutoff frequency = 9 kHz% · Sampling Frequency = 44.1 kHz (frequencies of interest 0 to 22.05 kHz)% · At least 60 dB of stopband attenuation% · No more than 1 dB passband ripplerp=1; rs=60; % specify passband ripple & stopband attenuation in dBf_spec=[7000 9000]; % specify passband and stopband edges in HzAA=[1 0]; %%% specfies that you want a lowpass filterdev=[(10^(rp/20)-1)/(10^(rp/20)+1) 10^(-rs/20)]; % parm. needed by design routineFs=44.1e3;[N,fo,ao,w]=remezord(f_spec,AA,dev,Fs) % estimates filter order and gives other parmsb=remez(N,fo,ao,w); % Computes the designed filter coefficients in vector b[H,ff]=freqz(b,1,1024,Fs); % Compute the frequency responsefigure; stem(0:N,b) % Plots filter's impulse responsefigure; subplot(2,1,1); plot(ff,20*log10(abs(H))) % Plot magnitude in dBsubplot(2,1,2); plot(ff,unwrap(angle(H))) % Plot unwrapped angle in radiansfigurezeros=roots(b); % compute roots of transfer function (i.e. compute zeros)plot(real(zeros), imag(zeros),'o') % Plot zeros in complex z planehold onplot(0,0,'x')theta=linspace(-pi,pi,10000);plot(cos(theta),sin(theta)) % Put