LAB MANUAL Of
DIGITAL SIGNAL PROCESSING
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
LIST OF EXPERIMETS 1) To determine frequency response of discrete LTI systems: a) FIR b) IIR 2) To determine and implement n-point DFT algorithms for a given discrete time sequence. 3) To determine circular convolution of two given sequences. 4) To determine and implement the Fast Fourier Transform algorithms 5) To compare the response of different types of window functions used for designing a digital filter. 6) To design an FIR low/high pass filter using windowing techniques 7) To design an IIR low/high pas filter using bilinear technique. 8) To perform interpolation operation on a given signal. 9) To perform decimation operation on a given signal. 10) To design and implement the FIR filter using TMS320C6713 for real time application..
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
EXPERIMENT NO-1 OBJECTIVE: To determine frequency response of discrete LTI systems: a) FIR b) IIR HARDWARE/SOFTWARE REQUIRED: MATLAB THEORY: [h,w] = freqz(b,a,w) returns the frequency response vector h and the corresponding angular frequency vector w for the digital filter whose transfer function is determined by the (real or complex) numerator and denominator polynomials represented in the vectors b and a, respectively. Eg. For general either FIR filter ( )=1+2
b=[1 2 3 2 1]
+3
a=[1 ] eg. For IIR filter H(z) = The vectors h and w are both of length l. The angular frequency vector w has values ranging from 0 to p radians per sample. When you don't specify the integer l, or you specify it as the empty vector [ ], the frequency response is calculated using the default value of 512 samples. freqz(b,a,...) plots the magnitude and unwrapped phase of the frequency response of the filter. The plot is displayed in the current figure window. PROGRAM: % FOR FIR FILTER b=[1 2 3 4]; a=[1 ]; w=-pi:pi; [h,w]=freqz(b,a,w); freqz(b,a,w); DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
% for IIR filter b=[1]; a=[1 2 2 2]; w=-pi:pi; freqz(b,a,w); RESULT/OBSERVATIONS: Freq. response of FIR filter
Magnitude (dB)
50 0 -50 -100 -150 -4
-3
-2 -1 0 1 2 Normalized Frequency ( rad/sample)
3
-3
-2 -1 0 1 2 Normalized Frequency ( rad/sample)
3
Phase (degrees)
2000 1000 0 -1000 -2000 -4
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
for IIR filter
Magnitude (dB)
10 0 -10 -20 -4
-3
-2 -1 0 1 2 Normalized Frequency ( rad/sample)
3
-3
-2 -1 0 1 2 Normalized Frequency ( rad/sample)
3
Phase (degrees)
4000 2000 0 -2000 -4000 -4
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
EXPERIMENT NO-2 OBJECTIVE: To determine and implement n-point DFT algorithms. HARDWARE/SOFTWARE REQUIREMENT: MATLAB THEORY: Spectral analysis is the process of identifying component frequencies in data. For discrete data, the computational basis of spectral analysis is the discrete Fourier transform (DFT). The DFT transforms time- or space-based data into frequency-based data. The DFT of a vector x of length n is another vector X of length n:
Its amplitude and phase are:
PROGRAM/CODE: close all; clear all; N=input('How many point DFT do you want?'); x2=input('Enter the sequence='); n2=length(x2); c= zeros(N); x2=[x2 zeros(1,N-n2)]; for k=1:N for n=1:N w=exp((-2*pi*i*(k-1)*(n-1))/N); %prev.step=>evaluating w-matrix x(n)=w; end c(k,:)=x; end r=[c]*[x2'] %plotting magnitude and angle subplot(211) stem(abs(r)); title('DFT-absolute value'); DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
subplot(212) stem(angle(r)); title('DFT-angle'); RESULT: How many point DFT do you want?4 Enter the sequence=[1 2 3 4] r= 10.0000 -2.0000 + 2.0000i -2.0000 - 0.0000i -2.0000 - 2.0000i DFT-absolute value
10
5
0
1
1.5
2
2.5
3
3.5
4
3
3.5
4
DFT-angle
4 2 0 -2 -4
1
1.5
2
2.5
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
EXPERIMENT NO- 3 OBJECTIVE: To determine circular convolution of two given sequences HARDWRE/SOFTWARE REQUIRED: MATLAB THEORY: Circular convolution is used to convolve two discrete Fourier transform (DFT) sequences. For very long sequences, circular convolution may be faster than linear convolution. The following exercises focus on the circular convolution property of the DFT. An N point circular convolution is defined as follows
Combining N-point vectors according to the circular convolution rule (Equation 1) is easy to visualizewith some experience, and Matlab provides the means to do the visualization. Circular convolutionrequires that all indexing be done in a circular (or periodic) fashion. Thus a shifting operation becomes a rotation. In the exercises below you will develop functions for circular flipping, circular shifting, and then use those functions to implement circular convolution. Method: Steps for cyclic convolution are the same as the usual convolution, except all index calculations are done "mod N" = "on the wheel" Steps for Cyclic Convolution Step1: “Plot f[m] and h[−m]
Subfigure 1.1
Subfigures 1.2
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Step 2: "Spin" h[−m] n times Anti Clock Wise (counter-clockwise) to get h[n-m] (i.e. Simply rotate the sequence, h[n], clockwise by n steps)
Figure 2: Step 2
Step 3: Pointwise multiply the f[m] wheel and the h[n−m] wheel. sum=y[n] Step 4: Repeat for all 0≤n≤N−1 Example 1: Convolve (n = 4)
Subfigure 3.1
Subfigure 3.2
Figure 3: Two discrete-time signals to be convolved.
h[−m] =
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Figure 4 Multiply f[m] and sum to yield: y[0] =3
h[1−m]
Figure 5 Multiply f[m] and sum to yield: y[1] =5
h[2−m]
Figure 6 Multiply f[m] and sum to yield: y[2] =3
h[3−m]
Figure 7 Multiply f[m] and sum to yield: y[3] =1
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
PROGRAM/CODE: x= input ('enter the first sequence'); h= input ('enter the second sequence'); N1= length(x); N2= length(h); N=max(N1,N2);%length of sequence x=[x zeros(1,N-N1)]; %modified first sequence h=[h zeros(1,N-N2)]; %modified second sequence
for n=0:N-1; y(n+1)=0; for i=0:N-1 j=mod(n-i,N); y(n+1)=y(n+1)+x(i+1)*h(j+1); %shifting and adding end end %to display and plot circular convolution n=1:N; disp('output sequence of circular convolution'); disp(y);%to view output in command window pause; stem(n,y);%plotting circular convolution grid minor; xlabel('time index'); ylabel('amplitude'); DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
title('circular convolution sequence of x and h'); RESULT : enter the first sequence[1 2 3 4] enter the second sequence[1 2 3 4] output sequence of circular convolution 26 28 26 20
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
EXPERIMENT NO: 4 OBJECTIVE: To determine and implement the Fast Fourier Transform algorithms HARDWARE/SOFTWRE REQUIRED: MATLAB THEORY: Although the DFT is computable transform, the straightforward implementation is very inefficient, especially when the sequence length N is large. In 1965, Cooley and Tukey showed the a procedure to substantially reduce the amount of computations involved in the DFT. This led to the explosion of applications of the DFT. All these efficient algorithms are collectively known as fast Fourier transform (FFT) algorithms. It can be implemented by using 2 methods: 1) Decimation in Time 2) Decimation in Frequency Let us consider the computation of the N = 2v point DFT by the divide-and conquer approach. We split the N-point data sequence into two N/2-point data sequences f1(n) and f2(n), corresponding to the even-numbered and odd-numbered samples of x(n), respectively, that is,
Thus f1(n) and f2(n) are obtained by decimating x(n) by a factor of 2, and hence the resulting FFT algorithm is called a decimation-in-time algorithm. Now the N-point DFT can be expressed in terms of the DFT's of the decimated sequences as follows:
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
But WN2 = WN/2. With this substitution, the equation can be expressed as
where F1(k) and F2(k) are the N/2-point DFTs of the sequences f1(m) and f2(m), respectively. Since F1(k) and F2(k) are periodic, with period N/2, we have F1(k+N/2) = F1(k) and F2(k+N/2) = F2(k). In addition, the factor WNk+N/2 = -WNk. Hence the equation may be expressed as
We observe that the direct computation of F1(k) requires (N/2)2 complex multiplications. The same applies to the computation of F2(k). Furthermore, there are N/2 additional complex multiplications required to compute WNkF2(k). Hence the computation of X(k) requires 2(N/2)2 + N/2 = N 2/2 + N/2 complex multiplications. This first step results in a reduction of the number of multiplications from N 2 to N 2/2 + N/2, which is about a factor of 2 for N large.
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Figure TC.3.1 First step in the decimation-in-time algorithm.
By computing N/4-point DFTs, we would obtain the N/2-point DFTs F1(k) and F2(k) from the relations
The decimation of the data sequence can be repeated again and again until the resulting sequences are reduced to one-point sequences. For N = 2v, this decimation can be
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
performed v = log2N times. Thus the total number of complex multiplications is reduced to (N/2)log2N. The number of complex additions is Nlog2N. For illustrative purposes, Figure TC.3.2 depicts the computation of N = 8 point DFT. We observe that the computation is performed in tree stages, beginning with the computations of four two-point DFTs, then two four-point DFTs, and finally, one eightpoint DFT. The combination for the smaller DFTs to form the larger DFT is illustrated in Figure TC.3.3 for N = 8.
Figure TC.3.2 Three stages in the computation of an N = 8-point DFT.
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Figure TC.3.3 Eight-point decimation-in-time FFT algorithm.
Figure TC.3.4 Basic butterfly computation in the decimation-in-time FFT algorithm. DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Decimation in Frequency: Decimation-in-frequency algorithm is obtained by using the divide-and-conquer approach. To derive the algorithm, we begin by splitting the DFT formula into two summations, one of which involves the sum over the first N/2 data points and the second sum involves the last N/2 data points. Thus we obtain
Now, let us split (decimate) X(k) into the even- and odd-numbered samples. Thus we obtain
Where we have used the fact that WN2 = WN/2 The computational procedure above can be repeated through decimation of the N/2-point DFTs X(2k) and X(2k+1). The entire process involves v = log2N stages of decimation, where each stage involves N/2 butterflies of the type shown in Figure TC.3.7. Consequently, the computation of the N-point DFT via the decimation-in-frequency FFT requires (N/2)log2N complex multiplications and Nlog2N complex additions, just as in the decimation-in-time algorithm. For illustrative purposes, the eight-point decimation-in-frequency algorithm is given in Figure TC.3.8.
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Figure TC.3.6 First stage of the decimation-in-frequency FFT algorithm.
Figure TC.3.7 Basic butterfly computation in the decimation-in-frequency.
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Figure TC.3.8 N = 8-piont decimation-in-frequency FFT algorithm.
We observe from Figure TC.3.8 that the input data x(n) occurs in natural order, but the output DFT occurs in bit-reversed order. We also note that the computations are performed in place. However, it is possible to reconfigure the decimation-in-frequency algorithm so that the input sequence occurs in bit-reversed order while the output DFT occurs in normal order. Furthermore, if we abandon the requirement that the computations be done in place, it is also possible to have both the input data and the output DFT in normal order.
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Y = fft(X) returns the discrete Fourier transform (DFT) of vector X, computed with a fast Fourier transform (FFT) algorithm. If X is a matrix, fft returns the Fourier transform of each column of the matrix. If X is a multidimensional array, fft operates on the first nonsingleton dimension. Y = fft(X,n) returns the n-point DFT. If the length of X is less than n, X is padded with trailing zeros to length n. If the length of X is greater than n, the sequence X is truncated. When X is a matrix, the lengths of the columns are adjusted in the same manner. PROGRAM/CODE: Fs = 1000;
% Sampling frequency
T = 1/Fs;
% Sample time
L = 1000;
% Length of signal
t = (0:L-1)*T;
% Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid y= 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); plot(Fs*t(1:50),y(1:50)); NFFT = 2^nextpow2(L); % Next power of 2 from length of y Y = fft(y,NFFT)/L; f = Fs/2*linspace(0,1,NFFT/2);
% Plot single-sided amplitude spectrum. plot(f,2*abs(Y(1:NFFT/2))) title('Single-Sided Amplitude Spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)|') WITHOUT FUNCTION function y=mditfft(x) x=input('Enter the Sequences') m=nextpow2(x); n=2^m; if length(x)
x=[x,zeros(1,n-length(x))]; end nxd=bin2dec(fliplr(dec2bin([1:n]-1,m)))+1; y=x(nxd); for mm=1:m le=2^mm;u=1; w=exp(-i*2*pi/le); for j=1:le/2 for k=j:le:n kp=k+le/2; t=y(kp)*u; y(kp)=y(k)-t; y(k)=y(k)+t; end u=u*w; end end RESULT/OBSRERVATION: Single-Sided Amplitude Spectrum of y(t)
1 0.9 0.8 0.7
|Y(f)|
0.6 0.5 0.4 0.3 0.2 0.1 0
0
50
100
150
200 250 300 Frequency (Hz)
350
400
450
500
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
EXPERIMENT NO: 5 OBJECTIVE: To compare the response of different types of window functions used for designing a digital filter. HARDWARE/SOFTWARE REQURIED: MATLAB, PC with windows OS XP/7 THEORY: Window functions The window method is most commonly used method for designing FIR filters. The simplicity of design process makes this method very popular. A window is a finite array consisting of coefficients selected to satisfy the desirable requirements. This chapter provides a few methods for estimating coefficients and basic characteristics of the window itself as well as the result filters designed using these coefficients. The point is to find these coefficients denoted by w[n]. When designing digital FIR filters using window functions it is necessary to specify:
A window function to be used; and The filter order according to the required specifications (selectivity and stopband attenuation).
These two requirements are interrelated. Each function is a kind of compromise between the two following requirements:
The higher the selectivity, i.e. the narrower the transition region; and The higher suppression of undesirable spectrum, i.e. the higher the stopband attenuation.
Table 2-3-1 below contains all window functions mentioned in this chapter and briefly compares their selectivity and stopband attenuation. Window function
Rectangular Triangular (Bartlett) BartlettHanning Hamming Blackman
Normalized length of the main lobe for N=20
Transition region for N=20
0.1π 0.2π
0.041π 0.11π
Minimum stopband attenuation of window function 13 dB 26 dB
Minimum stopband attenuation of designed filter 21 dB 26 dB
0.21π
0.13π
36 dB
39 dB
0.23π 0.32π
0.14π 0.2π
41 dB 58 dB
53 dB 75 dB
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Table 2-3-1. Comparison of window functions Special attention should be paid to the fact that minimum attenuation of window function and that of the filter designed using that function are different in most cases. The difference, i.e. additional attenuation occurs under the process of designing a filter using window functions. This affects the stopband attenuation to become additionally higher, which is very desirable. However, a drawback of this method is that the minimum stopband attenuation is fixed for each function. The exception is the Kaiser window described later in this chapter. The following concepts such as the main lobe, main lobe width, side lobes, transition region, minimum stopband attenuation of window function and minimum stopband attenuation of designed filter are described in more detail in Figure 2-3-1.
Figure 2-3-1. Main lobe, main lobe width, side lobes, transition region As can be seen in the table 2-3-1 above, the stopband attenuation of these windows is not adjustable. It is only possible to affect the transition region by increasing the filter order. For this reason it is preferable to start design process by specifying the appropriate window function on the basis of the stopband attenuation. It is most preferable to specify a window with the least stopband attenuation that satisfies the given requirements. This enables the designed filter to have the narrowest transition region. Example: Design a filter with the following characteristics: 1. Minimum stopband attenuation is 40dB; 2. Transition region between 2KHz and 3KHz; and 3. Sampling frequency 10KHz. Looking at the table 2-3-1 above (last column), it is obvious that required attenuation can be reached using the Hann window. Besides, such attenuation can be obtained using any other DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
window following the Hann window (all the functions are sorted according to their stopband attenuation), thus the higher stopband attenuation would be at a cost of the wider transition region for the same filter order. Since the filter specification also includes transition region, this would further result in increasing the filter order and its complexity as well. After specifying the suitable window function, it is necessary to compute the filter order. The cut-off frequencies of transition regions are normalized first.
The sampling frequency is fs=10KHz; The lower cut-off frequency of transition region is f1=2KHz; and The upper cut-off frequency of transition region is f2=3KHz.
Normalized frequencies are obtained in the following way:
Transition region of the Hann window is 0.12π for a 20th order filter. To satisfy the given specifications, the required transition region is expressed as: fn2 - fn1 = 0.2π Since the required transition region is wider than that of a 20th order filter, the filter order needs to be less than 20. It is found via iteration. Roughly estimated, the initial solution in this case can be a 12th order filter. The required order is somewhat higher and is found after performing several iterations. The resulting filter order is 16. Filter Designer Tool is used for obtaining a first-order filter. The table 2-3-2 provides performed iterations and other important information. Number of iteration Filter order Filter attenuation at 3KHz 1 12 25 dB 2 14 33 dB 3 16 49 dB 4 15 38 dB Table 2-3-2. Calculating filter order After specifying the window function and filter order, it is necessary to compute window coefficients w[n] using expressions for the specified window. DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Rectangular Window The rectangular window is rarely used for its low stopband attenuation. The first lobe (refer to Figure 2-3-2) has attenuation of 13dB and the narrowest transition region, therefore. A filter designed using this window has minimum stopband attenuation of 21 dB. Unlike other window functions being a kind of compromise between requirements for as narrow transition region as possible and as high stopband attenuation as possible, this window is characterized by extreme values. Namely, the minimum transition region is achieved, but at a cost of stopband attenuation. It is easy to find rectangular window coefficients as all coefficients between 0 and N-1 (N-filter order) are equal to 1, which can be expressed in the following way: w[n] = 1; 0 ≤ n ≤ N−1 Note that the rectangular window performs selection of N samples from a sequence of input samples, but it does not perform sample scaling. Figure 2-3-2 illustrates coefficients in the time-domain.
Figure 2-3-2. Rectangular window in the time domain Figure 2-3-3 illustrates the frequency domain of rectangular window
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Figure 2-3-3. Rectangular window in the frequency domain (spectrum) For its less stopband attenuation, the rectangular window is not preferable for digital filter design. Such a less attenuation is a result of cut-off samples within a window (a sequence of sampled frequencies). Up to a zero sample (from which sampling starts), all sampled frequencies are equal to zero. The first sample represents a sudden jump to some value (non-zero sample). Exactly these sudden jumps result in producing relatively sharp high-frequency components which lessen the stopband attenuation. The attenuation gets higher by making cut-off samples less sharp, which results in reducing filter selectivity, i.e. wider transition region. Since initial requirements of a digital filter are predefined and due to less selectivity, it is necessary to increase the filter order to narrow the transition region. Note the fact that the transition region is inversely proportional to the filter order N. The transition region narrows as the filter order increases. Increase in filter order affects the filter to become more complex and need more time for sample processing. This is why it is very important to be careful when specifying the window function and filter order as well. Triangular (Bartlett) window The triangular (Bartlett) window is one among many functions that lessens the effects of final samples. Due to it, the stopband attenuation of this window is higher than that of the rectangular window, whereas the selectivity is less. Namely, filters designed using this window have wider transition region than those designed using the rectangular window. As a result, it is necessary to have a higher order filter in order to keep the same tansition region as that of the filters designed with the rectangular window. This is the price to pay for producing the higher attenuation. This function also represents a kind of compromise between requirements for as narrow transition region as possible and as higher stopband attenuation as possible, where the transition region is considered more important characteristic.
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
One of the advantages of designing filter using the triangular window is the simplicity of computing coefficients. The triangular window coefficients can be expressed as:
Figure 2-3-4 illustrates coefficients in the time-domain.
Figure 2-3-4. Triangular window coefficients in the time-domain Figure 2-3-5. illustrates the coefficients spectrum of the triangular window shown in Figure 2-34.
Figure 2-3-5. Triangular (Bartlett) window in the frequency domain (spectrum)
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
The attenuation of triangular window is low for most digital filter applications, but it is considerably higher than that for rectangular window. In some cases, when high attenuations are not needed, this filter can be used because it provides an easy way of computing coefficients. Hamming Window The Hamming window is one of the most popular and most commonly used windows. A filter designed with the Hamming window has minimum stopband attenuation of 53dB, which is sufficient for most implementations of digital filters. The transition region is somewhat wider than that of the Hann and Bartlett-Hanning windows, whereas the stopband attenuation is considerably higher. Unlike minimum stopband attenuation, the transition region can be changed by changing the filter order. The transition region narrows, whereas the minimum stopband attenuation remains unchanged as the filter order increases. The Hamming window coefficients are expressed as:
The Hamming window belongs to a class of generalized cosine functions which are described in more details at the end of this chapter. Figure 2-3-10 illustrates the Hamming window coefficients in the time domain.
Figure 2-3-10. The Hamming window coefficients in the time domain
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Figure 2-3-11. The Hamming window coefficients in the frequency domain Figure 2-3-11 illustrates the Hamming window in the frequency domain. As seen, the first side lobe is attenuated, so that the minimum stopband attenuation is defined in terms of the second side lobe and amounts to 41dB. All side lobes have almost the same maximum values (about 45dB). Blackman Window The Blackman window is, along with Kaiser, Hamming and Blackman-Harris windows, considered most commonly used and the most popular windows. Relatively high attenuation makes this window very convenient for almost all applications. The minimum stopband attenuation of a filter designed with this window amounts to 75dB. The Blackman window coefficients are expressed as:
Figure 2-3-14 illustrates the Blackman window coefficients in the time domain, whereas Figure 2-3-15. illustrates its coefficients in the frequency domain.
Figure 2-3-14. The Blackman window coefficients in the time domain DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Figure 2-3-15. The Blackman window coefficients in the frequency domain As seen from Figure 2-3-15, the Blackman window frequency domain reminds of the Hann window frequency domain. The difference is in attenuation of the first side lobe which amounts to 51dB, as well as in the main lobe which is somewhat wider. The side lobes, following the first one, cause additional stopband attenuation. Kaiser Window As you know, each of the windows described above is a kind of compromise between requirements for as narrow transition region as possible (greater selectivity) and as higher stopband attenuation as possible. Comparing Hann and Bartlett-Hanning windows, it is obvious that both of them have the same transition region, but the Bartlett-Hanning window has higher attenuation. There is one more thing of concern which says that the minimum stopband attenuation depends on the specified window, whereas an increase in filter order affects the transition region. All this leads us to the conclusion that the windows described here are not optimal. An optimal window is a function that has maximum attenuation according to the given width of the main lobe. The optimal window is also known as Kaiser window. Its coefficients are expressed as:
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
where aa is the minimum stopband attenuation, and Δ ω is the width of (normalized) transition region. The order of band-pass and band-stop filters, obtained from the expression above, should be multiplied by 2. The value of parameter β can be obtained from the table 2-3-4. aa β less than 21 0 between 21 and 50 0.5842(aa - 21)^0.4 + 0.07886(aa - 21) more than 50 0.1102(aa - 8.7) Table 2-3-4. Values of parameter β I0(*) is a modified zero order Bessel function of the first kind. It can be approximated via expression:
In most practical cases it is sufficient to consider the first 20 elements of this order. As can be seen from all mentioned above, in order to design an optimal Kaiser filter it is necessary to know normalized width of transition region as well as minimum desirable stopband attenuation. DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Example: What would a filter design with Kaiser window look like if we take into consideration the requirements given at the beginning of this chapter? The result of a filter designed with the Hann window is a 16th order filter. Assuming that it is required to design a filter with the following characteristics: 1. Minimum stopband attenuation is 40dB (aa); 2. Transition region is between 2KHz and 3KHz (f1, f2); and 3. Sampling frequency is 10KHz (fs). Transition region:
The specified filter order is N=12. It means that the needed filter order is less by 4 than that obtained with the Hann window. It is important to say that all the expressions above are obtained in an empirical way, which means that there will be some exceptions and variations in practice. For example, the expression used to compute the filter order gives accurate results in approximately 98% cases. Otherwise, the resulting filter order should be changed. Fortunately, the changes to be made are slight, and the filter order is increased or decreased by 2 at most.
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
EXPERIMENT NO: 6 OBJECTIVE: To design an FIR low/high pass filter using windowing techniques SOFTWARE/HARDWARE REQUIRED: MATLAB THEORY: FIR filter specification First of all, it is necessary to learn the basic concepts that will be used further in this book. You should be aware that without being familiar with these concepts, it is not possible to understand analyses and synthesis of digital filters. Figure 2-2-1 illustrates a low-pass digital filter specification. The word specification actually refers to the frequency response specification.
Figure 2-2-1a. Low-pass digital filter specification
Figure2-2-1b. Low-pass digital filter specification
ωp – normalized cut-off frequency in the passband; ωs – normalized cut-off frequency in the stopband; δ1 – maximum ripples in the passband; δ2 – minimum attenuation in the stopband [dB]; ap – maximum ripples in the passband; and as – minimum attenuation in the stopband [dB]. DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Frequency normalization can be expressed as follows:
Where:
fs is a sampling frequency; f is a frequency to normalize; and ω is normalized frequency.
Specifications for high-pass, band-pass and band-stop filters are defined almost the same way as those for low-pass filters. Figure 2-2-2 illustrates a high-pass filter specification.
Figure 2-2-2a. High-pass digital filter specification
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Figure 2-2-2b. High-pass digital filter specification Comparing these two figures 2-2-1 and 2-2-2, it is obvious that low-pass and high-pass filters have similar specifications. The same values are defined in both cases with the difference that in the later case the pass-band is substituted by the stop-band and vice versa. Figure 2-2-3 illustrates a band-pass filter specification.
Figure 2-2-3a. Band-pass digital filter specification
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Figure 2-2-3b. Band-pass digital filter specification Figure 2-2-4 illustrates a band-stop digital filter specification.
Figure 2-2-4a. Band-stop digital filter specification
Figure 2-2-4b. Band-stop digital filter specification DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
The rectangular window sequence is given by
In the design of FIR filters using any window technique, the order can be calculated using the formula given by
where δp is the passband ripple, δs is the stopband ripple, fp is the passband frequency, fs is the stopband frequency and Fs is the sampling frequency. BOXCAR Boxcar window: BOXCAR still works but maybe removed in the future. Use RECTWIN instead. RECTWIN Rectangular window: W = RECTWIN(N) returns the N-point rectangular window. FIR1 FIR filter design using the window method. B = FIR1(N,Wn) designs an N'th order lowpass FIR digital filter and returns the filter coefficients in length N+1 vector B. The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0 corresponding to half the sample rate. The filter B is real and has linear phase. The normalized gain of the filter at Wn is -6 dB. B = FIR1(N,Wn,'high') designs an N'th order highpass filter.You can also use B = FIR1(N,Wn,'low') to design a lowpass filter. If Wn is a two-element vector, Wn = [W1 W2], FIR1 returns an order N bandpass filter with passband W1 < W < W2. You can also specify B = FIR1(N,Wn,'bandpass'). If Wn = [W1 W2], B = FIR1(N,Wn,'stop') will design a bandstop filter.
CEIL Round towards plus infinity. CEIL(X) rounds the elements of X to the nearest integers towards infinity.
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
PROGRAM: clc; close all; clear all; format long; rp=input('enter the passband ripple'); rs=input('enter the stopband ripple'); fp=input('enter the passband frequency'); fs=input('enter the stopband frequency'); f=input('enter the sampling frequency'); wp=2*(fp/f); ws=2*(fs/f); num=-20*log10(sqrt(rp*rs))-13; dem=14.6*(fs-fp)/f; n=ceil(num/dem); n1=n+1; if(rem(n,2)~=0) n1=n; n=n-1; end; y=boxcar(n1); %Lowpass filter b=fir1(n,wp,y); [h,o]=freqz(b,1,256); m=20*log10(abs(h)); subplot(2,1,1); plot(o/pi,m); ylabel('gain in db---->'); xlabel('Normalised frequency---->'); title('FIR filter using Rectangular window of LPF ----'); grid on;
%Highpass filter b=fir1(n,wp,'high',y); [h,o]=freqz(b,1,256); m=20*log10(abs(h)); subplot(2,1,2); plot(o/pi,m); ylabel('gain in db---->'); xlabel('Normalised frequency---->'); title('FIR filter using Rectangular window of HPF ----'); DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
title('Magnitude response of high pass filter'); grid on; RESULT: enter the passband ripple
0.04
enter the stopband ripple
0.05
enter the passband frequency 1500 enter the stopband frequency 2000 enter the sampling frequency 9000
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
EXPERIMENT-7 OBJECTIVE: To design an IIR low/high pas filter using bilinear technique. SOFTWARE/HARDWRE REQUIRED: MATLAB THEORY: Bilinear transform is a correction of the backwards difference method.The bilinear transform (also known as Tustin's transformation) is defined as the substitution.It is the most popular method.The bilinear transform produces a digital filter whose frequency response has the same characteristics as the frequency response of the analogue filter (but its impulse response may then be quite different). The bilinear transform:
In some sources you will see the factor (2/T) multiplying the RHS of the bilinear transform; this is an optional scaling, but it cancels and does not affect the final result. The steps of the bilinear transform method are as follows: 1. “Warp” the digital critical (e.g. band-edge or "corner") frequencies ωi , in other words compute the corresponding analogue critical frequencies Ω i = tan(ωi/2). 2. Design an analogue filter which satisfies the resulting filter response specification. 3. Apply the bilinear transform
to the s-domain transfer function of the analogue filter to generate the required z-domain transfer function.
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
PROGRAM: clc; clear all; close all; %Low pass Filter [N,Wn]=buttord(2,4.83,-3-15, 's'); [b,a]=butter(N,Wn,'s'); [num,den]=bilinear(b,a,1); figure; freqz(num,den); %High Pass Filter [N,Wn]=buttord(2,4.83,-3,-15, 's'); [b,a]=butter(N,Wn,'high','s'); [num,den]=bilinear(b,a,1); figure;freqz(num,den);
RESULT:
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
low pass filter
Magnitude (dB)
0 -50 -100 -150
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 0.8 Normalized Frequency ( rad/sample)
0.9
1
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 0.8 Normalized Frequency ( rad/sample)
0.9
1
Phase (degrees)
0 -50 -100 -150 -200
high pass filter
Magnitude (dB)
0 -100 -200 -300 -400
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 0.8 Normalized Frequency ( rad/sample)
0.9
1
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 0.8 Normalized Frequency ( rad/sample)
0.9
1
Phase (degrees)
200 150 100 50 0
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
EXPERIMENT NO-8 OBJECTIVE: To perform interpolation operation on a given signal. SOFTWARE/HARDWRE REQUIRED: MATLAB THEORY: yi = interp1(x,Y,xi) interpolates to find yi, the values of the underlying function Y at the points in the vector or array xi. x must be a vector. Y can be a scalar, a vector, or an array of any dimension, subject to the following conditions:
If Y is a vector, it must have the same length as x. A scalar value for Y is expanded to have the same length as x. xi can be a scalar, a vector, or a multidimensional array, and yi has the same size as xi. If Y is an array that is not a vector, the size of Y must have the form [n,d1,d2,...,dk], where n is the length of x. The interpolation is performed for each d1-by-d2-by-...-dk value in Y. The sizes of xi and yi are related as follows: o If xi is a scalar or vector, size(yi) equals [length(xi), d1, d2, ..., dk]. o If xi is an array of size [m1,m2,...,mj], yi has size [m1,m2,...,mj,d1,d2,...,dk].
yi = interp1(Y,xi) assumes that x = 1:N, where N is the length of Y for vector Y, or size(Y,1) for matrix Y. yi = interp1(x,Y,xi,method) interpolates using alternative methods: 'nearest' Nearest neighbor interpolation 'linear'
Linear interpolation (default)
'spline'
Cubic spline interpolation
'pchip'
Piecewise cubic Hermite interpolation
'cubic'
(Same as 'pchip')
'v5cubic' Cubic interpolation used in MATLAB® 5. This method does not extrapolate. Also, if x is not equally spaced, 'spline' is used/ For the 'nearest', 'linear', and 'v5cubic' methods, interp1(x,Y,xi,method) returns NaN for any element of xi that is outside the interval spanned by x. For all other methods, interp1 performs extrapolation for out of range values. yi = interp1(x,Y,xi,method,'extrap') uses the specified method to perform extrapolation for out of range values. yi = interp1(x,Y,xi,method,extrapval) returns the scalar extrapval for out of range values. NaN and 0 are often used for extrapval. PROGRAM:
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
x = 0:10; y = sin(x); xi = 0:.25:10; yi = interp1(x,y,xi); plot(x,y,'o',xi,yi) RESULT: 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1
0
1
2
3
4
5
6
7
8
9
10
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
EXPERIMENT NO.-9 OBJECTIVE: To perform decimation operation on a given signal. SOFTWARE/HARDWARE REQUIRED: MATLAB THEORY: Decimation reduces the original sampling rate for a sequence to a lower rate, the opposite of interpolation. The decimation process filters the input data with a lowpass filter and then resamples the resulting smoothed signal at a lower rate. y = decimate(x,r) reduces the sample rate of x by a factor r. The decimated vector y is r times shorter in length than the input vector x. By default, decimate employs an eighth-order lowpass Chebyshev Type I filter with a cutoff frequency of 0.8*(Fs/2)/r. It filters the input sequence in both the forward and reverse directions to remove all phase distortion, effectively doubling the filter order. y = decimate(x,r,n) uses an order n Chebyshev filter. Orders above 13 are not recommended because of numerical instability. In this case, a warning is displayed. PROGRAM/CODE: t = 0:.00025:1;
% Time vector
x = sin(2*pi*30*t) + sin(2*pi*60*t); y = decimate(x,4); %View the original and decimated signals: stem(x(1:120)), axis([0 120 -2 2]) % Original signal title('Original Signal') figure stem(y(1:30))
% Decimated signal
title('Decimated Signal') RESULT:
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Original Signal
2 1.5 1 0.5 0 -0.5 -1 -1.5 -2
0
20
40
60
80
100
120
20
25
30
Decimated Signal
2 1.5 1 0.5 0 -0.5 -1 -1.5 -2
0
5
10
15
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
EXPERIMENT NO.-10 OBJECTIVE: To design and implement the FIR filter using TMS320C6713 for real time application.. SOFTWARE/HARDWARE REQUIRED: MATLAB Finite Impulse Response Filter DESIGNING AN FIR FILTER : Following are the steps to design linear phase FIR filters Using Windowing Method. I.
Clearly specify the filter specifications. Eg: Order = 30; Sampling Rate = 8000 samples/sec Cut off Freq. = 400 Hz.
II.
Compute the cut-off frequency Wc Eg: Wc = 2*pie* fc / Fs = 2*pie* 400/8000 = 0.1*pie
III.
Compute the desired Impulse Response h d (n) using particular Window Eg: b_rect1=fir1(order, Wc , 'high',boxcar(31));
IV.
Convolve input sequence with truncated Impulse Response x (n)*h (n)
USING MATLAB TO DETERMINE FILTER COEFFICIENTS : Using FIR1 Function on Matlab B = FIR1(N,Wn) designs an N'th order lowpass FIR digital filter and returns the filter coefficients in length N+1 vector B. The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0 corresponding to half the sample rate. The filter B is real and has linear phase, i.e., even symmetric coefficients obeying B(k) =B(N+2-k), k = 1,2,...,N+1. If Wn is a two-element vector, Wn = [W1 W2], FIR1 returns an order N bandpass filter with passband W1 < W < W2. B = FIR1(N,Wn,'high') designs a highpass filter. B= FIR1(N,Wn,'stop') is a bandstop filter if Wn = [W1 W2]. If Wn is a multi-element vector, Wn = [W1 W2 W3 W4 W5 ... WN], FIR1 returns an order N multiband filter with bands 0 < W < W1, W1 < W < W2, ..., WN < W < 1. B = FIR1(N,Wn,'DC-1') makes the first band a passband. B = FIR1(N,Wn,'DC-0') makes the first band a stopband. DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
For filters with a passband near Fs/2, e.g., highpass and bandstop filters, N must be even. By default FIR1 uses a Hamming window. Other available windows, including Boxcar, Hanning, Bartlett, Blackman, Kaiser and Chebwin can be specified with an optional trailing argument. For example, B = FIR1(N,Wn,kaiser(N+1,4)) uses a Kaiser window with beta=4. B = FIR1(N,Wn,'high',chebwin(N+1,R)) uses a Chebyshev window. By default, the filter is scaled so the center of the first pass band has magnitude exactly one after windowing. Use a trailing 'noscale' argument to prevent this scaling, e.g. B = FIR1(N,Wn,'noscale'), B = FIR1(N,Wn,'high','noscale'), B = FIR1(N,Wn,wind,'noscale'). Matlab Program to generate ‘FIR Filter-Low Pass’ Coefficients using FIR1 % FIR Low pass filters using rectangular, triangular and kaiser windows % sampling rate - 8000 order = 30; cf=[500/4000,1000/4000,1500/4000]; cf--> contains set of cut-off frequencies[Wc ] % cutoff frequency - 500 b_rect1=fir1(order,cf(1),boxcar(31)); b_tri1=fir1(order,cf(1),bartlett(31)); b_kai1=fir1(order,cf(1),kaiser(31,8));
Rectangular Triangular Kaisar [Where 8-->Beta Co-efficient]
% cutoff frequency - 1000 b_rect2=fir1(order,cf(2),boxcar(31)); b_tri2=fir1(order,cf(2),bartlett(31)); b_kai2=fir1(order,cf(2),kaiser(31,8)); % cutoff frequency - 1500 b_rect3=fir1(order,cf(3),boxcar(31)); b_tri3=fir1(order,cf(3),bartlett(31)); b_kai3=fir1(order,cf(3),kaiser(31,8)); fid=fopen('FIR_lowpass_rectangular.txt','wt'); fprintf(fid,'\t\t\t\t\t\t%s\n','Cutoff -400Hz'); fprintf(fid,'\nfloat b_rect1[31]={'); fprintf(fid,'%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,\n',b_rect1); fseek(fid,-1,0); fprintf(fid,'};'); fprintf(fid,'\n\n\n\n'); fprintf(fid,'\t\t\t\t\t\t%s\n','Cutoff -800Hz');
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
fprintf(fid,'\nfloat b_rect2[31]={'); fprintf(fid,'%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,\n',b_rect2); fseek(fid,-1,0); fprintf(fid,'};'); fprintf(fid,'\n\n\n\n'); fprintf(fid,'\t\t\t\t\t\t%s\n','Cutoff -1200Hz'); fprintf(fid,'\nfloat b_rect3[31]={'); fprintf(fid,'%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,\n',b_rect3); fseek(fid,-1,0); fprintf(fid,'};'); fclose(fid); winopen('FIR_highpass_rectangular.txt');
T.1 : Matlab generated Coefficients for FIR Low Pass Kaiser filter: Cutoff -500Hz float b_kai1[31]={-0.000019,-0.000170,-0.000609,-0.001451,-0.002593,-0.003511,0.003150,0.000000,0.007551,0.020655,0.039383,0.062306,0.086494,0.108031,0.122944, 0.128279,0.122944,0.108031,0.086494,0.062306,0.039383,0.020655,0.007551,0.000000, -0.003150,-0.003511,-0.002593,-0.001451,-0.000609,-0.000170,-0.000019};
IMPLEMENTATION OF AN FIR FILTER : Cutoff -1000HzTO IMPLEMENT : ALGORITHM float b_kai2[31]={-0.000035,-0.000234,-0.000454,0.000000,0.001933,0.004838,0.005671, We need to realize an advance FIR filter by implementing its difference equation as per the specifications. -0.000000,-0.013596,-0.028462,-0.029370,0.000000,0.064504,0.148863,0.221349,0.249983, A direct form I implementation approach is taken. (The filter coefficients are taken as ai as generated by 0.221349,0.148863,0.064504,0.000000,-0.029370,-0.028462,-0.013596,-0.000000,0.005671, the Matlab program.) 0.004838,0.001933,0.000000,-0.000454,-0.000234, -0.000035};
Cutoff -1500Hz float b_kai3[31]={-0.000046,-0.000166,0.000246,0.001414,0.001046,-0.003421,-0.007410, 0.000000,0.017764,0.020126,-0.015895,-0.060710,-0.034909,0.105263,0.289209,0.374978, 0.289209,0.105263,-0.034909,-0.060710,-0.015895,0.020126,0.017764,0.000000,-0.007410, -0.003421,0.001046,0.001414,0.000246,-0.000166, DEPARTMENT OF ELECTRONICS AND-0.000046}; COMMUNICATION ENGINEERING
T.2 :Matlab generated Coefficients for FIR Low Pass Rectangular filter Cutoff -500Hz float b_rect1[31]={-0.008982,-0.017782,-0.025020,-0.029339,-0.029569,-0.024895, -0.014970,0.000000,0.019247,0.041491,0.065053,0.088016,0.108421,0.124473,0.134729, 0.138255,0.134729,0.124473,0.108421,0.088016,0.065053,0.041491,0.019247,0.000000, -0.014970,-0.024895,-0.029569,-0.029339,-0.025020,-0.017782,-0.008982};
Cutoff -1000Hz float b_rect2[31]={-0.015752,-0.023869,-0.018176,0.000000,0.021481,0.033416,0.026254,-0.000000,0.033755,-0.055693,-0.047257,0.000000,0.078762,0.167080,0.236286,0.262448, 0.236286,0.167080,0.078762,0.000000,-0.047257,-0.055693,-0.033755,-0.000000,0.026254, FLOWCHART FOR FIR : 0.033416,0.021481,0.000000,-0.018176,-0.023869,-0.015752};
Cutoff -1500Hz
T.3 Matlab generated Coefficients for FIR Low Pass Triangular filter float: b_rect2[31]={-0.020203,-0.016567,0.009656,0.027335,0.011411,-0.023194,-0.033672, Cutoff -500Hz 0.000000,0.043293,0.038657,-0.025105,-0.082004,-0.041842,0.115971,0.303048,0.386435, float b_tri1[31]={0.000000,-0.001185,-0.003336,-0.005868,-0.007885,-0.008298,-0.005988, 0.303048,0.115971,-0.041842,-0.082004,-0.025105,0.038657,0.043293,0.000000,-0.033672, 0.000000,0.010265,0.024895,0.043368,0.064545,0.086737,0.107877,0.125747,0.138255, -0.023194,0.011411,0.027335,0.009656,-0.016567,-0.020203}; 0.125747,0.107877,0.086737,0.064545,0.043368,0.024895,0.010265,0.000000,-0.005988, -0.008298,-0.007885,-0.005868,-0.003336,-0.001185,0.000000};
Cutoff -1000Hz float b_tri2[31]={0.000000,-0.001591,-0.002423,0.000000,0.005728,0.011139,0.010502, -0.000000,-0.018003,-0.033416,-0.031505,0.000000,0.063010,0.144802,0.220534,0.262448, 0.220534,0.144802,0.063010,0.000000,-0.031505,-0.033416,-0.018003,-0.000000,0.010502,
MATLAB Program to generate ‘FIR Filter-High Pass’ Coefficients using FIR1 0.011139,0.005728,0.000000,-0.002423,-0.001591,0.000000};
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING Cutoff -1500Hz float b_tri3[31]={0.000000,-0.001104,0.001287,0.005467,0.003043,-0.007731,-0.013469,
% FIR High pass filters using rectangular, triangular and kaiser windows % sampling rate - 8000 order = 30; cf=[400/4000,800/4000,1200/4000];
;cf--> contains set of cut-off frequencies[Wc]
% cutoff frequency - 400 b_rect1=fir1(order,cf(1),'high',boxcar(31)); b_tri1=fir1(order,cf(1),'high',bartlett(31)); b_kai1=fir1(order,cf(1),'high',kaiser(31,8)); Where Kaiser(31,8)--> '8'defines the value of 'beta'. % cutoff frequency - 800 b_rect2=fir1(order,cf(2),'high',boxcar(31)); b_tri2=fir1(order,cf(2),'high',bartlett(31)); b_kai2=fir1(order,cf(2),'high',kaiser(31,8)); % cutoff frequency - 1200 b_rect3=fir1(order,cf(3),'high',boxcar(31)); b_tri3=fir1(order,cf(3),'high',bartlett(31)); b_kai3=fir1(order,cf(3),'high',kaiser(31,8)); fid=fopen('FIR_highpass_rectangular.txt','wt'); fprintf(fid,'\t\t\t\t\t\t%s\n','Cutoff -400Hz'); fprintf(fid,'\nfloat b_rect1[31]={'); fprintf(fid,'%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,\n',b_rect1); fseek(fid,-1,0); fprintf(fid,'};'); fprintf(fid,'\n\n\n\n'); fprintf(fid,'\t\t\t\t\t\t%s\n','Cutoff -800Hz'); fprintf(fid,'\nfloat b_rect2[31]={'); fprintf(fid,'%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,\n',b_rect2); fseek(fid,-1,0); fprintf(fid,'};'); fprintf(fid,'\n\n\n\n'); fprintf(fid,'\t\t\t\t\t\t%s\n','Cutoff -1200Hz'); fprintf(fid,'\nfloat b_rect3[31]={'); fprintf(fid,'%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,\n',b_rect3); fseek(fid,-1,0); fprintf(fid,'};'); fclose(fid); winopen('FIR_highpass_rectangular.txt');
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
T.1 : MATLAB generated Coefficients for FIR High Pass Kaiser filter: Cutoff -400Hz float b_kai1[31]={0.000050,0.000223,0.000520,0.000831,0.000845,-0.000000,-0.002478, -0.007437,0.015556,-0.027071,-0.041538,-0.057742,-0.073805,-0.087505,-0.096739, 0.899998,-0.096739,0.087505,-0.073805,-0.057742,-0.041538,-0.027071,-0.015556, -0.007437,-0.002478,0.000000,0.000845,0.000831,0.000520,0.000223,0.000050}; Cutoff -800Hz float b_kai2[31]={0.000000,-0.000138,-0.000611,-0.001345,-0.001607,-0.000000,0.004714, 0.012033,0.018287,0.016731,0.000000,-0.035687,-0.086763,-0.141588,-0.184011,0.800005,0.184011,-0.141588,-0.086763,-0.035687,0.000000,0.016731,0.018287,0.012033,0.004714, ALGORITHM TO IMPLEMENT : 0.000000,-0.001607,-0.001345,-0.000611,-0.000138,0.000000};
IMPLEMENTATION OF AN FIR FILTER :
Cutoff -1200Hz We need to realize an advance FIR filter by implementing its difference equation as per the specifications. Afloat directb_kai3[31]={-0.000050,-0.000138,0.000198,0.001345,0.002212,-0.000000,-0.006489,-0.012033,form I implementation approach is taken. (The filter coefficients are taken as ai as generated by 0.005942,0.016731,0.041539,0.035687,-0.028191,-0.141589,-0.253270,0.700008,-0.253270,the Matlab program.) 0.141589,-0.028191,0.035687,0.041539,0.016731,-0.005942,-0.012033,-0.006489, -0.000000,0.002212,0.001345,0.000198,-0.000138,-0.000050};
T.2 :MATLAB generated Coefficients for FIR High Pass Rectangular filter Cutoff -400Hz float b_rect1[31]={0.021665,0.022076,0.020224,0.015918,0.009129,-0.000000,-0.011158,-0.023877,0.037558,-0.051511,-0.064994,-0.077266,-0.087636,-0.095507,-.100422,0.918834,-0.100422,0.095507,-0.087636,-0.077266,-0.064994,-0.051511,-0.037558,-0.023877,-0.011158,0.000000,0.009129,0.015918,0.020224,0.022076,0.021665}; Cutoff -800Hz float b_rect2[31]={0.000000,-0.013457,-0.023448,-0.025402,-0.017127,-0.000000,0.020933, 0.038103, 0.043547, 0.031399,0.000000,-0.047098,-0.101609,-0.152414,-0.188394,0.805541,-0.188394,.152414, -0.101609,-0.047098,0.000000,0.031399,0.043547,0.038103,0.020933,-0.000000,-0.017127,0.025402,-0.023448,-0.013457,0.000000}; Cutoff -1200Hz float b_rect3[31]={-0.020798,-0.013098,0.007416,0.024725,0.022944,-0.000000,-0.028043,-0.037087,0.013772,0.030562,0.062393,0.045842,-0.032134,-0.148349,-0.252386,0.686050,-0.252386,0.148349,-0.032134,0.045842,0.062393,0.030562,-0.013772,-0.037087,-0.028043,0.000000,0.022944,0.024725,0.007416,-0.013098,-0.020798}; DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
T.3 : MATLAB generated Coefficients for FIR High Pass Triangular filter Cutoff -400Hz float b_tri1[31]={0.000000,0.001445,0.002648,0.003127,0.002391,-0.000000,-0.004383,-0.010943,0.019672,-0.030353,-0.042554,-0.055647,-0.068853,-0.081290,-0.092048,0.902380,-0.092048,0.081290,-0.068853,-0.055647,-0.042554,-0.030353,-0.019672,-0.010943,-0.004383,-0.000000, 0.002391, 0.003127,0.002648,0.001445,0.000000}; Cutoff -800Hz float b_tri2[31]={0.000000,-0.000897,-0.003126,-0.005080,-0.004567,-0.000000,0.008373, 0.017782, 0.023225, 0.018839, 0.000000, -0.034539, -0.081287, -0.132092, -0.175834, 0.805541,-0.175834, 0.132092, -0.081287, -0.034539, 0.000000, 0.018839, 0.023225, 0.017782, 0.008373,-0.000000,0.004567, -0.005080, -0.003126, -0.000897, 0.000000}; Cutoff -1200Hz float b_tri3[31]={0.000000,-0.000901,0.001021,0.005105,0.006317,-0.000000,-0.011581,-0.017868,0.007583,0.018931,0.042944,0.034707,-0.026541,-0.132736,-0.243196,0.708287,-0.243196,0.132736,-0.026541,0.034707,0.042944,0.018931,-0.007583,-0.017868,-0.011581,-0.000000, 0.006317, 0.005105, 0.001021, -0.000901, 0.000000};
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
C PROGRAM TO IMPLEMENT FIR FILTER:
fir.c #include "filtercfg.h" #include "dsk6713.h" #include "dsk6713_aic23.h" float filter_Coeff[] ={0.000000,-0.001591,-0.002423,0.000000,0.005728, 0.011139,0.010502,-0.000000,-0.018003,-0.033416,-0.031505,0.000000, 0.063010,0.144802,0.220534,0.262448,0.220534,0.144802,0.063010,0.000000, -0.031505,-0.033416,-0.018003,-0.000000,0.010502,0.011139,0.005728, 0.000000,-0.002423,-0.001591,0.000000 }; static short in_buffer[100]; DSK6713_AIC23_Config config = {\ 0x0017, /* 0 DSK6713_AIC23_LEFTINVOL Leftline input channel volume */\ 0x0017, /* 1 DSK6713_AIC23_RIGHTINVOL Right line input channel volume*/\ 0x00d8, /* 2 DSK6713_AIC23_LEFTHPVOL Left channel headphone volume */\ 0x00d8, /* 3 DSK6713_AIC23_RIGHTHPVOL Right channel headphone volume */\ 0x0011, /* 4 DSK6713_AIC23_ANAPATH Analog audio path control */\ 0x0000, /* 5 DSK6713_AIC23_DIGPATH Digital audio path control */\ 0x0000, /* 6 DSK6713_AIC23_POWERDOWN Power down control */\ 0x0043, /* 7 DSK6713_AIC23_DIGIF Digital audio interface format */\ 0x0081, /* 8 DSK6713_AIC23_SAMPLERATE Sample rate control */\ 0x0001 /* 9 DSK6713_AIC23_DIGACT Digital interface activation */ \ }; /* * main() - Main code routine, initializes BSL and generates tone */ void main() { DSK6713_AIC23_CodecHandle hCodec; Uint32 l_input, r_input,l_output, r_output; /* Initialize the board support library, must be called first */ DSK6713_init(); /* Start the codec */ hCodec = DSK6713_AIC23_openCodec(0, &config); DSK6713_AIC23_setFreq(hCodec, 1);
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
while(1) { /* Read a sample to the left channel */ while (!DSK6713_AIC23_read(hCodec, &l_input)); /* Read a sample to the right channel */ while (!DSK6713_AIC23_read(hCodec, &r_input)); l_output=(Int16)FIR_FILTER(&filter_Coeff ,l_input); r_output=l_output; /* Send a sample to the left channel */ while (!DSK6713_AIC23_write(hCodec, l_output)); /* Send a sample to the right channel */ while (!DSK6713_AIC23_write(hCodec, r_output)); } /* Close the codec */ DSK6713_AIC23_closeCodec(hCodec); } signed int FIR_FILTER(float * h, signed int x) { int i=0; signed long output=0; in_buffer[0] = x; /* new input at buffer[0] */ for(i=29;i>0;i--) in_buffer[i] = in_buffer[i-1]; /* shuffle the buffer */ for(i=0;i<31;i++) output = output + h[i] * in_buffer[i]; return(output); }
PROCEDURE :
Switch on the DSP board. Open the Code Composer Studio. Create a new project Project New (File Name. pjt , Eg: FIR.pjt) Initialize on board codec. Note: “Kindly refer the Topic Configuration of 6713 Codec using BSL” DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Add the given above ‘C’ source file to the current project (remove codec.c source file from the project if you have already added). Connect the speaker jack to the input of the CRO. Build the program. Load the generated object file(*.out) on to Target board. Run the program Observe the waveform that appears on the CRO screen. Vary the frequency on function generator to see the response of filter.
MATLAB GENERATED FREQUENCY RESPONSE High Pass FIR filter(Fc= 800Hz).
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
Low Pass FIR filter (Fc=1000Hz)
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING