MATLAB Script for Filter Coefficients Parameterization, Papers of Digital Signal Processing

Download MATLAB Script for Filter Coefficients Parameterization and more Papers Digital Signal Processing in PDF only on Docsity! ECE 491 6/5/2007 DSP & Quantization 1 Tanim Taher ECE 491-153 Fall 05 Undergraduate Research: Digital Signal Processing & Quantization Effects Done By: Tanim Taher SID# 10370800 Date: December 19, 2005 ECE 491 6/5/2007 DSP & Quantization 2 Tanim Taher ACKNOWLEDGEMENTS I would like to express my gratitude towards Dr. Erdal Oruklu and Dr. Jafar Saniiee for their strong support in making this project successful and a good learning experience for me. A special thanks to Dr. Saniiee for realizing my potential and suggesting that I should take up this research course. Thank you professors so much! Tanim Taher ECE 491 6/5/2007 DSP & Quantization 5 Tanim Taher 2. INTRODUCTION Digital Signal Processing is distinguished from the rest of the digital computer field by the unique type of data it uses: signals. In most cases, these signals originate as sensory data from the real world, for example, sound signals being picked up by a microphone. Digital filters are a very important part of DSP. In fact, their extraordinary performance is one of the key reasons that DSP has become so popular. Filters have two uses: signal separation and signal restoration: 1) Signal Separation is the separation of signals that have been combined. Signal separation is needed when a signal has been contaminated with interference, noise, or other signals undesired. For example, when talking in a telephone’s speakerphone, in addition to the person’s voice at this end of the line, the microphone picks up the sound being produced by the loudspeaker. This would normally cause the person at the other end of the line to hear echoes of his or her own voice. Thus a filter is used that separates the sounds and eliminates the echo. 2) Signal restoration is used when a signal has been distorted in some way. For example, an audio recording made with poor equipment may be filtered to better represent the sound as it actually occurred. Filtering can be done with analog filters, but digital filters achieve much better results. This is because fine tuning analog filters can be difficult, as the analog components have tolerance issues with their stated property values and the actual ones. On the other hand, digital filters can be designed with precise characteristics. However, digital filters have their own limitations. The major limit is due to a finite word length. This limits the degree of precision with which the signal can be represented. Digital filtering is carried out by arithmetic operations (processing) that involve ‘filter ECE 491 6/5/2007 DSP & Quantization 6 Tanim Taher coefficients’. The finite word length also means that the results of the processing are not 100% accurate, but only an approximation depending on the word length and representation format. There are two main kinds of digital filters – infinite impulse response (IIR) filters, and finite impulse response (FIR) filters. IIR filters have feedback, while FIR filters do not. Any digital filter can be represented as a polynomial transfer function. The coefficients of the polynomial represent the digital filter’s coefficients mentioned before. The transfer function is usually represented as follows: H(z) = (b0 + b1z-1 + b2z-2 + … bMz-M) / (1 + a1z-1 + a2z-2 + … aNz-N) Here z-1 is called a unit delay operator. It represents a digital component whouse output value is the same as the input value of the previous bus cycle. An FIR filter has the denominator (a) coefficients as all zero, while IIR filters have at least one non-zero ‘a’ coefficients. The values of the ‘a’ and ‘b’ coefficients are determined based on the filtering characteristics that are required. Now these coefficients must also be represented in the computer in a digital format with finite precision, and hence the response of the filter can vary depending on how precisely these numbers are represented in the digital signal processor. In summary, quantization of the signal data, the coefficients and the finite precision possible in arithmetic processing affect the performance of digital filters. The purpose of my research was to look at the quantization effects more closely. An important objective of my work was also to learn more about the field of digital signal processing and digital signal processors. Hence a large portion of this report has been dedicated to summarizing what I have learnt. ECE 491 6/5/2007 DSP & Quantization 7 Tanim Taher 3. TYPES OF PROCESSING There are two main types of processing in computer science – real-time processing and off-line processing. Off-line processing is where data has been obtained before, but the processing is done later. Real-time processing is where the data is processed immediately after input. Off-line processing is used in situations where the results of the operations on the data collected are not immediately required. However, in real-time processing, the results need to be obtained immediately. For example, in speech to text computer systems, the spoken voice (input signal) must be processed in real-time to generate the text (result) when preparing a document in a word processing software. Time is a very crucial factor in real-time systems. Since the results are required immediately, digital filters in real-time applications need to be very fast. Hence, sometimes shortcuts may need to be made during processing, in order to produce results faster. If a computer system is being used in an application as a real time digital filter, this may mean that single precision word format is used instead of double precision, as single precision processing is faster. On the other-hand, time is not so crucial a factor in off-line processing. Hence, if we use the same computer in an off-line application that needs to use a digital filter, we can now use double precision arithmetic. Thus we can use higher precision in the digital filter which makes quantization effects minimal. Thus depending on the type of processing required – real-time or off-line, the quantization effects on digital filters can vary. Off course we can use the same amount of higher precision in real-time digital filters, but this may make the design of the digital processor more complex, and hence more expensive. ECE 491 6/5/2007 DSP & Quantization 10 Tanim Taher 6. QUANTIZATION ISSUES No digital system in the world can give us infinite precision while representing all numbers. There is always a size limit for a digital word. Rounding is the term used to describe the process of approximating a number in the digital domain. This is done by rounding the least significant bit such that the value of the digital number is closest to the actual value of the number. We are loosing some degree of precision by this approximation, and generating what is called a ‘roundoff error’. On average, this error amount is equal to a quarter of the least significant bit value. During analog to digital sampling, the roundoff error can be quite significant if the amount of noise inherent in the original signal is close to the average error amount. Another type of error is truncation error which occurs when the results of arithmetic operations like division are truncated after the word length has been reached. In recursive filtering operations, the quantization errors simply keep adding up each time with the digitally filtered data. Furthermore, roundoff error can cause stability issues when it comes to IIR filters. Rounding, as it causes loss in precision, changes the filter coefficients such that they no longer are the ideal values. If the pole (denominator in the transfer function) coefficients are close to 1, and rounding produces a coefficient equal to 1 or more, then the digital filter is unstable. There are three binary number formats that are in common use. The type of format used changes the impact of quantization effects in digital filters. The following are the formats: 1) Single Precision floating point – this uses a 32 bit word. The format for storing the number is standard number format (e.g +1.234 X 103 is a standard format number). 2) Double Precision floating point – 64 bit, standard number format. 3) Fixed Point – varying word lengths are in use (common 16, 32, 64 bits), where the decimal point is fixed and does not vary. ECE 491 6/5/2007 DSP & Quantization 11 Tanim Taher 7. SIGNAL FLOW GRAPHS Previously, the transfer function for any digital filer was shown. The transfer function can be used to draw the direct form system diagrams. The system diagram simply represents the flow of data within a digital filter, and shows the operations carried out at each stage, like multiplication by the coefficients, and summation of the results. The direct form diagrams are very common, and I was introduced to them while taking the course, Digital Signal Processing (ECE 437). Since my focus in this report is on material that I learnt afterwards, we will not get into the details of the direct form diagrams. However, we shall look into signal flow graphs, which are another way of representing the digital filter pictorially. Signal flow graphs use nodes to represent the adders and as signal sources. If an arrow is pointing towards a node, it means that that node is adding that signal, or performing a process on the signal (like multiplication). If the signal arrow is pointing away, then it simply means that the result of the node’s operation is taken out in that data path. Multiplication coefficients are represented simply as numbers on the data flow lines, and the unit delay as z-1on the line. The following diagram makes this clear: ECE 491 6/5/2007 DSP & Quantization 12 Tanim Taher Based on this principle, an N order FIR filter would look like the following if shown in the form of a signal flow graph: Reference: http://www.onesmartclick.com/engineering/digital-signal-processing.html ECE 491 6/5/2007 DSP & Quantization 15 Tanim Taher Example: if A2(z) = 1 + 0.2z-1 + 0.3z-2 then B2(z) = 0.3 + 0.2z-1 + 1z-2 so α2(2) = β 2(0) = 0.3 And α2(1) = β 2(1) = 0.2 And α2(0) = β 2(2) = 1 First Step: Km = αm(m) Since A2(z) = 1 + 0.2z-1 + 0.3z-2 Km = K2 = αm(m) = α2(2) = 0.3 Second: αm-1(n) = [αm(n) - Km βm(n)] / (1 – K2m) α1(0) = α2(0) – K2 β2(0)] / (1 – K22) = [1 – (.3 X .3)] / [1 - .32] = 1 α1(1) = α2(1) – K2 β2(1)] / (1 – K22) = [.2 – (.3 X .2)] / [1 - .32] = .1538 So K1 = α1(1) = 0.1538 All the Am(n), for n = 0 to m, need to be determined so that the Bm(n) can be found for all required filter orders. Bm(z) is also required later on to find the V coefficients. Now that Bm(z) is known for all m orders, we can find the V coefficients. To do this, we first create an FIR filter by copying the numerator of the transfer function of the original IIR. Let us call this filter Cm(z). Then: Let cm(n) be the numerator coefficients. It is the response of an ‘m’ order FIR filter at time nT. i.e. hm(n)=cm(n) First Step: vm = cm(m) Second: Cm-1(n) = Cm(n) - vm βm(n) Thus both the V and K coefficients can be found in this way. ECE 491 6/5/2007 DSP & Quantization 16 Tanim Taher Now the question is why use lattice filters. It is indeed quite a lot of work to determine the K and V coefficients initially. However, during processing, a bigger problem occurs as is evident from the system diagram shown earlier. The problem is that a large number of extra multiplication and summation operations need to be carried out in lattice filters in comparison to their direct form counterparts. So there needs to be a good reason for choosing lattice filters. The answer is that the output of a lattice filter is less sensitive to the parameterization of the filter coefficients. This is specially relevant to stability, as otherwise slight changes in the coefficients during parameterization can make the entire filter system unstable. The following two diagrams explain this concept. Source of diagrams: homes.esat.kuleuven.ac.be/~moonen/ This figure shows all the possible pole locations for any IIR filter where the coefficient numbers are represented in 5 bits. Note due to areas of concentration, generally, slight changes in coefficients equal a significant change in the filter’s characteristics. This figure shows all the possible pole locations for any IIR lattice filter where the K coefficient numbers are represented in 5 bits. Note due to uniform concentration, slight changes in coefficients here do not affect much the filter’s characteristics. -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 ECE 491 6/5/2007 DSP & Quantization 17 Tanim Taher 10a. RESULTS: Simulation Program Overview In an effort to explore the quantization effects in digital filters and to obtain a first-hand experience in programming digital filters in a computer, I developed a Digital Filter Simulator and Analysis Tool. I used Matlab to develop the software, and used Matlab’s Graphics User Interface Design Environment (GUIDE) to make the user interface. Most of my research time was spent on this software development, and in my opinion the investment was worthwhile. The software lets the user easily change the quantization parameters, and at the click of a button change the structure of a digital filter. The output graphs are very useful in analyzing the performance of the filters in response to quantization parameter changes. The software has to GUIDE forms/menus. The first form is called the “Waveform Generator” and it is used to generate sample data to be fed into the digital filter. The advantage of this form is that virtually any periodic waveform can be used as input provided the fourier series coefficients are known. This is also very helpful to the electrical engineering student for two reasons: a) He/she can observe firsthand how the fourier series can be used to generate any waveform of choice, and hence understand the power of Fourier series. b) He/she can also observe the effect of aliasing which occurs significantly if the sampling rate is below the Nyquist frequency. The other form is the “Digital Filter” form. The user inputs the coefficients of a digital filter (FIR or IIR) from the transfer function coefficients. The user also chooses the filter structure to be used. The quantization parameters like word length, fixed-point/floating point can be changed at will. Then at the click of a button, the software simulates the digital filter. There are two outputs that can be produced by the program: a) Data output: Using the data-set of the waveform generated by the “Waveform Generator”, an output data-set is created and displayed by digitally filtering the ECE 491 6/5/2007 DSP & Quantization 20 Tanim Taher The entry fields on the right side let the user determine the length of the data-set by specifying the start and stop times for the input waveform. The “Total Window Points” field is used to zoom in or zoom out, and it refers to the number of data points to be plotted on screen. This screenshot is that of the “Digital Filter” form. ECE 491 6/5/2007 DSP & Quantization 21 Tanim Taher The quantization options lets the user to select the word type (fixed point, single or double floating points), and the word length for the fixed point type. The output shown is that produced by the digital filter with input being the waveform shown in the “Waveform Generator” shown before. The digital filter is specified by the ‘a’ and ‘b’ coefficients that the user can input. The filter structure is chosen by selecting any bullet from the “Implementation Type” box area. The most important thing to note is that as the quantization options are changed, the actual filter coefficients used by the digital filter changes, and these are displayed in the two list boxes at the bottom of the form (A Coefficients, B coefficients). The following screenshot shows the magnitude phase response which is the other output. ECE 491 6/5/2007 DSP & Quantization 22 Tanim Taher Note, the aliasing effect can be observed here also, by lowering the ‘Fs(Hz)’ entry-field value. This is shown in the next screenshot: The third output waveform is for the phase response: ECE 491 6/5/2007 DSP & Quantization 25 Tanim Taher 10c. RESULTS: Software Algorithm Description 1) The Waveform Generator uses the following algorithm to generate the input waveform: Time Period, T = 1 / Sampling Frequency a. For sample n, input x(n) = 0 b. T0 = Start time c. Start with harmonic number h = 1 d. For harmonic number h, x(n) = [h’th fourier Coeff] * Sin (h * fundamental frequency * n * T + T0) + x(n) e. Increment h by 1 f. If h <= highest harmonic number, repeat steps c to f. g. Increment n h. If n not last sample, then repeat steps c to h. The last sample count is determined from the start and stop times, and the time period T. 2) Once the quantization options are chosen, the program creates a fixed point data object type is created. Then the program uses this object type for all arithmetic, if the word type is fixed point. Otherwise, it uses single or double floating point objects (as specified by the user) for all calculations. Hence during filtering, the simulated digital filter is fully exposed to quantization effects. 3) Direct Form 1 Filter Implementation: The software uses two separate loops to perform the two separate summation operations of the IIR direct form 1 implementation equation shown below: ECE 491 6/5/2007 DSP & Quantization 26 Tanim Taher 4) Direct Form 2 Filter Implementation: The software uses two separate loops to perform this filter function. Each of the loops perform one of the summation functions as specified by the IIR direct form 2 implementation equations shown below: 5) Generating Lattice Filter Coefficients: The software uses several loops to generate the lattice coefficients recursively using the equations given before in the Lattice Filters sections. While I was able to successfully generate the lattice coefficients, I faced constant difficulty recreating the filter model for the lattice structure, and never was able to successfully implement it. Hence, I used the lattice coefficient vectors generated by my algorithm, with Matlab’s built in lattice filter function (keyword latcfilt()) for the filtering. 6) In order to generate the magnitude response, the software uses complex arithmetic. This algorithm uses the quantized coefficients, and replaces ‘z’ within the filter’s transfer function (given before) by z = e-jwt where, w is the sweep frequency being used. Then the value of H(z) , i.e. H(e-jwt) is calculated. This complex value’s magnitude represents the magnitude response, and its phase represents the phase response. ECE 491 6/5/2007 DSP & Quantization 27 Tanim Taher 10d. RESULTS: Software Applicability This was briefly mentioned before. Basically, the software is a powerful learning tool. I certainly learned a lot as the programmer. By simulating a lot of filters using the software, I was also able to learn and see the quantization effects, and to classify how different filter structures hold up to quantization effects. For an electrical engineering student, the software can help for the following: Learning tools for digital signal processing – aliasing, fourier series, filter structures, stability, magnitude and phase response, etc. Students can use and modify this code to learn about: MATLAB Simulation Filter structure Students can test Quantization effects First Hand. By examining code, then can appreciate how Quantization errors build up. I would like to allow the source to be used as freeware, that can be modified by students and users. The condition is that the comments with my ID information are not modified. ECE 491 6/5/2007 DSP & Quantization 30 Tanim Taher 10kHz LP, Order 4 Chebyshev IIR, Double, Floating Point, Direct Form1 8 Bit, Fixed Point, DF1 ECE 491 6/5/2007 DSP & Quantization 31 Tanim Taher 10kHz HP,Order 13 Chebyshev IIR, 6 Bit, Fixed Point, Lattice 8 Bit, Fixed Point, Lattice ECE 491 6/5/2007 DSP & Quantization 32 Tanim Taher 11. INTERPRETATION: A lot was learned about quantization effects from the simulations. Listed below are the interpretation of the simulation results: Lattice Structure is most stable. Direct Form 1 and Direct Form 2 can give issue with stability. Direct Form 1 and Direct Form2 sensitive to slight parameterization changes. Cascade and Parallel structures (done in MATLAB) more stable than Direct Form 1 and Direct Form 2. This is due to the fact that the lower order structural units used in the cascade and parallel forms have lesser feedback continuation (only one level of feedback for second order units). The total number of coefficients in the IIR and FIR filter is the main factor that determines the filter performance. That is the order is a very important factor affecting filter performance. As long as the word length is adequate to hold the maximum and minimum input and output values, word format has little effect on performance. However, for input and output values with large ranges, the floating point data formats were decidedly better than the fixed point ones. The fixed point outputs often saturated. Even though lattice structure is more stable, the extra processing required (more multiplication and addition) cannot be really justified. This is because today, having large word lengths an inexpensive feat, and with large words, the direct form filters give performance close to that of lattice filters when parameterization effects are considered. ECE 491 6/5/2007 DSP & Quantization 35 Tanim Taher 13. APPENDIX: Program Listing The following is the listing for the “Waveform Generator” form. %Programmer Name: Tanim Taher %e-mail: tanim_02@yahoo.com %Done as coursework for ECE 491, undergraduate research work %Fall 2005 function varargout = wavegen(varargin) % WAVEGEN M-file for wavegen.fig % WAVEGEN, by itself, creates a new WAVEGEN or raises the existing % singleton*. % % H = WAVEGEN returns the handle to a new WAVEGEN or the handle to % the existing singleton*. % % WAVEGEN('CALLBACK',hObject,eventData,handles,...) calls the local % function named CALLBACK in WAVEGEN.M with the given input arguments. % % WAVEGEN('Property','Value',...) creates a new WAVEGEN or raises the % existing singleton*. Starting from the left, property value pairs are % applied to the GUI before wavegen_OpeningFunction gets called. An % unrecognized property name or invalid value makes property application % stop. All inputs are passed to wavegen_OpeningFcn via varargin. % % *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one % instance to run (singleton)". % % See also: GUIDE, GUIDATA, GUIHANDLES % Copyright 2002-2003 The MathWorks, Inc. % Edit the above text to modify the response to help wavegen % Last Modified by GUIDE v2.5 07-Dec-2005 16:00:42 % Begin initialization code - DO NOT EDIT gui_Singleton = 1; gui_State = struct('gui_Name', mfilename, ... 'gui_Singleton', gui_Singleton, ... 'gui_OpeningFcn', @wavegen_OpeningFcn, ... 'gui_OutputFcn', @wavegen_OutputFcn, ... 'gui_LayoutFcn', [] , ... 'gui_Callback', []); if nargin && ischar(varargin{1}) ECE 491 6/5/2007 DSP & Quantization 36 Tanim Taher gui_State.gui_Callback = str2func(varargin{1}); end if nargout [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); else gui_mainfcn(gui_State, varargin{:}); end % End initialization code - DO NOT EDIT % --- Executes just before wavegen is made visible. function wavegen_OpeningFcn(hObject, eventdata, handles, varargin) global h_max ampl h_no Amp %intialize harmonic values h_max = 1; ampl = 1; Amp(1) = 1 h_no = 1; % Choose default command line output for wavegen handles.output = hObject; % Update handles structure guidata(hObject, handles); % UIWAIT makes wavegen wait for user response (see UIRESUME) % uiwait(handles.figure1); % --- Outputs from this function are returned to the command line. function varargout = wavegen_OutputFcn(hObject, eventdata, handles) % Get default command line output from handles structure varargout{1} = handles.output; function edit4_Callback(hObject, eventdata, handles) % --- Executes during object creation, after setting all properties. function edit4_CreateFcn(hObject, eventdata, handles) % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end function edit5_Callback(hObject, eventdata, handles) global h_max ampl h_no Amp h_no = str2double(get(handles.edit5,'String')); % Get entered harmonic no value ECE 491 6/5/2007 DSP & Quantization 37 Tanim Taher if h_no <= h_max set(handles.edit6,'String',Amp(h_no)); else set(handles.edit6,'String',0); end % --- Executes during object creation, after setting all properties. function edit5_CreateFcn(hObject, eventdata, handles) % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end function edit6_Callback(hObject, eventdata, handles) global h_max ampl h_no Amp %redundant % --- Executes during object creation, after setting all properties. function edit6_CreateFcn(hObject, eventdata, handles) % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes on button press in pushbutton3. function pushbutton3_Callback(hObject, eventdata, handles) global h_max ampl h_no Amp ampl = str2double(get(handles.edit6,'String')); % Get entered amplitude value h_no = str2double(get(handles.edit5,'String')); % Get entered harmonic no value Amp(h_no)=ampl; set(handles.edit5,'string',h_no+1); set(handles.edit6,'string',0); if h_no > h_max h_max = h_no set(handles.maxH,'string',h_max) end; %Find Nyquist frequency and set it as default value min_samp = (2 * h_max * str2double(get(handles.edit4,'String'))) / (2 * pi); w = str2double(get(handles.edit7,'String')); if w < min_samp set(handles.edit7,'string',min_samp); end ECE 491 6/5/2007 DSP & Quantization 40 Tanim Taher if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes on button press in clearH. function clearH_Callback(hObject, eventdata, handles) global h_max ampl h_no Amp global st_time end_time pts tot_pts xs dpts T% for slider pts = 0; %Make all the fourier coefficients = 0, clear the plot tot_pts = 0; clear dpts; clear xs; h_max = 1; clear Amp; h_no = 1; set(handles.edit5,'string',1); set(handles.edit6,'string',1); set(handles.maxH,'string',1); Amp(1) = 1; plot(0,0); THE FOLLOWING IS THE CODE FOR THE DIGITAL FILTER SIMULATOR %Programmer Name: Tanim Taher %e-mail: tanim_02@yahoo.com %Done as coursework for ECE 491, undergraduate research work %Fall 2005 function varargout = Filter_Coefficients_Parameterization(varargin) % FILTER_COEFFICIENTS_PARAMETERIZATION M-file for Filter_Coefficients_Parameterization.fig % FILTER_COEFFICIENTS_PARAMETERIZATION, by itself, creates a new FILTER_COEFFICIENTS_PARAMETERIZATION or raises the existing % singleton*. % % H = FILTER_COEFFICIENTS_PARAMETERIZATION returns the handle to a new FILTER_COEFFICIENTS_PARAMETERIZATION or the handle to % the existing singleton*. % % FILTER_COEFFICIENTS_PARAMETERIZATION('CALLBACK',hObject,eventData,handl es,...) calls the local % function named CALLBACK in FILTER_COEFFICIENTS_PARAMETERIZATION.M with the given input arguments. % % FILTER_COEFFICIENTS_PARAMETERIZATION('Property','Value',...) creates a new FILTER_COEFFICIENTS_PARAMETERIZATION or raises the % existing singleton*. Starting from the left, property value pairs are % applied to the GUI before Filter_Coefficients_Parameterization_OpeningFunction gets called. An % unrecognized property name or invalid value makes property application ECE 491 6/5/2007 DSP & Quantization 41 Tanim Taher % stop. All inputs are passed to Filter_Coefficients_Parameterization_OpeningFcn via varargin. % % *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one % instance to run (singleton)". % % See also: GUIDE, GUIDATA, GUIHANDLES % Copyright 2002-2003 The MathWorks, Inc. % Edit the above text to modify the response to help Filter_Coefficients_Parameterization % Last Modified by GUIDE v2.5 14-Dec-2005 19:24:45 % Begin initialization code - DO NOT EDIT gui_Singleton = 1; gui_State = struct('gui_Name', mfilename, ... 'gui_Singleton', gui_Singleton, ... 'gui_OpeningFcn', @Filter_Coefficients_Parameterization_OpeningFcn, ... 'gui_OutputFcn', @Filter_Coefficients_Parameterization_OutputFcn, ... 'gui_LayoutFcn', [] , ... 'gui_Callback', []); if nargin && ischar(varargin{1}) gui_State.gui_Callback = str2func(varargin{1}); end if nargout [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); else gui_mainfcn(gui_State, varargin{:}); end % End initialization code - DO NOT EDIT % --- Executes just before Filter_Coefficients_Parameterization is made visible. function update_listbox(handles) %This functions displays the A and B coefficients global a_no a_max b_no b_max c_A c_B a_string b_string oCA oCB oCK oCV clear a_string b_string a2_string b2_string k_string v_string b2_string(1)={' '}; a2_string(1) = {'A0 = 1'}; for n=1:a_max %Show the a coeff tmp = {strcat('A',num2str(n),' = ',num2str(c_A(n)))}; a_string(n) = tmp; end for n=1:b_max %Show the b coeffs tmp = {strcat('B',num2str(n-1),' = ',num2str(c_B(n)))}; b_string(n) = tmp; ECE 491 6/5/2007 DSP & Quantization 42 Tanim Taher end if get(handles.Latt, 'Value') == 0 set(handles.paramAK, 'Title','A Coefficients'); set(handles.paramBV, 'Title','B Coefficients'); for n=1:numel(oCA) tmp = {strcat('A',num2str(n),' = ',num2str(oCA(n)))}; a2_string(n+1) = tmp; end for n=1:numel(oCB) tmp = {strcat('B',num2str(n-1),' = ',num2str(oCB(n)))}; b2_string(n) = tmp; end set(handles.listAK,'String',a2_string); set(handles.listBV,'String',b2_string); else %Show the K and V coeffs if filter type is lattice set(handles.paramAK, 'Title','K Coefficients'); set(handles.paramBV, 'Title','V Coefficients'); for n=1:numel(oCK) tmp = {strcat('K',num2str(n),' = ',num2str(oCK(n)))}; k_string(n) = tmp; end for n=1:numel(oCV) tmp = {strcat('V',num2str(n-1),' = ',num2str(oCV(n)))}; v_string(n) = tmp; end set(handles.listAK,'String',k_string); set(handles.listBV,'String',v_string); end set(handles.a_list,'String',a_string); set(handles.b_list,'String',b_string); function Filter_Coefficients_Parameterization_OpeningFcn(hObject, eventdata, handles, varargin) % This function has no output args, see OutputFcn. % hObject handle to figure % varargin command line arguments to Filter_Coefficients_Parameterization (see VARARGIN) global a_no a_max b_no b_max c_A c_B PT WL FL oCA oCB %Initialize oCA=[]; oCB=[]; PT = 1; b_no = 1; b_max = 2; a_no = 1; a_max = 2; WL = 8; FL = 4; %c_B = [0.028 0.053 0.071 0.053 0.028] % c_A = [-2.026 2.148 -1.159 0.279] c_B = [0.5 .8125] %Initialize the A and B coefficients to any desired value c_A = [0.5 .3125] pushbutton4_Callback(0, 0, handles) ECE 491 6/5/2007 DSP & Quantization 45 Tanim Taher % hObject handle to listbox1 (see GCBO) if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes during object creation, after setting all properties. function edit3_CreateFcn(hObject, eventdata, handles) % hObject handle to edit3 (see GCBO) if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes during object creation, after setting all properties. function edit4_CreateFcn(hObject, eventdata, handles) % hObject handle to edit4 (see GCBO) if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes during object creation, after setting all properties. function listbox2_CreateFcn(hObject, eventdata, handles) % hObject handle to listbox2 (see GCBO) if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end function wordlength_Callback(hObject, eventdata, handles) % hObject handle to wordlength (see GCBO) global WL FL WL = str2num(get(hObject,'String')) FL = str2num(get(handles.frac,'String')) if WL <= FL WL = FL + 2; set(hObject,'String',WL) end % --- Executes during object creation, after setting all properties. function wordlength_CreateFcn(hObject, eventdata, handles) % hObject handle to wordlength (see GCBO) if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end ECE 491 6/5/2007 DSP & Quantization 46 Tanim Taher function frac_Callback(hObject, eventdata, handles) % hObject handle to frac (see GCBO) global WL FL WL = str2num(get(handles.wordlength,'String')) FL = str2num(get(handles.frac,'String')) if WL <= FL WL = FL + 2; set(handles.wordlength,'String',WL) end % --- Executes during object creation, after setting all properties. function frac_CreateFcn(hObject, eventdata, handles) % hObject handle to frac (see GCBO) if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes on button press in togglebutton1. function togglebutton1_Callback(hObject, eventdata, handles) % Hint: get(hObject,'Value') returns toggle state of togglebutton1 global PT %Change the word type used for processing PT = PT + 1 % Sets the Float or fixed Point data types if PT > 3 PT = 1 end if PT == 1 set(handles.togglebutton1,'String','Fixed Point'); elseif PT == 2 set(handles.togglebutton1,'String','Single, Floating Point'); elseif PT == 3 set(handles.togglebutton1,'String','Double, Floating Point'); end if PT == 1 set(handles.frac,'Enable','on'); set(handles.wordlength,'Enable','on'); else set(handles.frac,'Enable','off'); set(handles.wordlength,'Enable','off'); end function f_start_Callback(hObject, eventdata, handles) % hObject handle to f_start (see GCBO) % --- Executes during object creation, after setting all properties. function f_start_CreateFcn(hObject, eventdata, handles) % hObject handle to f_start (see GCBO) if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) ECE 491 6/5/2007 DSP & Quantization 47 Tanim Taher set(hObject,'BackgroundColor','white'); end function f_end_Callback(hObject, eventdata, handles) % hObject handle to f_end (see GCBO) % --- Executes during object creation, after setting all properties. function f_end_CreateFcn(hObject, eventdata, handles) % hObject handle to f_end (see GCBO) if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes on button press in pushbutton6. function pushbutton6_Callback(hObject, eventdata, handles) sweep_freq(handles,'magnitude'); % hObject handle to pushbutton6 (see GCBO) % --- Executes on button press in pushbutton7. function pushbutton7_Callback(hObject, eventdata, handles) sweep_freq(handles,'phase'); % hObject handle to pushbutton7 (see GCBO) %Calculates the lattice coefficients function calc_latt(handles) global oCK oCV a_max b_max c_A c_B K_C V_C Km Am K_C = []; V_C = []; oCK = []; oCV = []; % X_C are the ones used for calculations, oCX are used for printing Km=[]; Am=[]; Bm=[]; %My Code for determining K and V Am(a_max,1) = 1; for i=1:a_max Am(a_max,i+1) = c_A(i); end for m=a_max:-1:1 %From m'th order to 1st order Am(m,1) = 1; %Am(0)=1 Km(m) = Am(m,m+1); %Km=Am(m) if m>1 %If there is any lower order coefficients to find for n=1:(m-1) %for Am-1(1) to Am-1(m-1) Am(m-1,n+1)=(Am(m,n+1)-Km(m)*Am(m,m-n+1))/(1- (Km(m)^2)); %The Equation for finding lower order A coeff's end end end %This finds all values of Bm(k): ECE 491 6/5/2007 DSP & Quantization 50 Tanim Taher for i = 1:b_max B_C(i) = double(c_B(i)); oCB(i) = B_C(i); end for i = 1:numel(oCK) K_C(i) = oCK(i); end for i = 1:numel(oCV) V_C(i) = oCV(i); end end update_listbox(handles) % hObject handle to pushbutton4 (see GCBO) % --- Executes on button press in del_coeff. function del_coeff_Callback(hObject, eventdata, handles) clear c_A c_B global a_no b_no a_max b_max c_A c_B b_no = 1; b_max = 1; a_no = 1; a_max = 1; c_A(1) = 0; c_B(1) = 0; set(handles.ac_no,'String',1) %set(handles.ac,'String',0) set(handles.bc_no,'String',0) %set(handles.bc,'String',0) update_listbox(handles); % plot(xs(1:tot_pts),dpts(1:tot_pts)) function ac_no_Callback(hObject, eventdata, handles) function ac_Callback(hObject, eventdata, handles) %Store a coefficient values global a_no a_max b_no b_max c_A c_B if a_no > 0 c_A(a_no) = str2double(get(hObject,'String')); a_no = a_no + 1; set(handles.ac_no,'String',a_no); if a_no <= a_max set(handles.ac,'String',c_A(a_no)); else set(handles.ac,'String',0); a_max = a_no; c_A(a_no) = 0; end end update_listbox(handles) % hObject handle to ac (see GCBO) ECE 491 6/5/2007 DSP & Quantization 51 Tanim Taher function bc_no_Callback(hObject, eventdata, handles) % hObject handle to bc_no (see GCBO) function bc_Callback(hObject, eventdata, handles) %Store B coefficient values from guide global a_no a_max b_no b_max c_A c_B if b_no > 0 c_B(b_no) = str2double(get(hObject,'String')); b_no = b_no + 1; if b_no < b_max set(handles.bc_no,'String',b_no-1); set(handles.bc,'String',c_B(b_no)); else set(handles.bc,'String',0); b_max = b_no; c_B(b_no) = 0; end end update_listbox(handles) % --- Executes on selection change in a_list. function a_list_Callback(hObject, eventdata, handles) global a_max a_no c_A; a_no = get(hObject,'Value'); set(handles.ac_no,'String',a_no); set(handles.ac,'String',c_A(a_no)); % --- Executes on selection change in b_list. function b_list_Callback(hObject, eventdata, handles) global b_max b_no c_B; b_no = get(hObject,'Value'); set(handles.bc_no,'String',b_no-1); set(handles.bc,'String',c_B(b_no)); function param_obj(handles) global PT FL if strcmp(get(handles.togglebutton1,'String'),'Fixed Point') PT = 1; else PT = 2; end function find_lattice_op(handles) %Failed function that tried to implement lattice structure %Calculate the lattice output global oCK oCV PT FL T1 dpts tot_pts F otpt xs oCK oCV [oCK,oCV]=tf2latc(c_B,[1, c_A]); %Fm(n)=Fo(n), G(1,n)=Go(1,n) ECE 491 6/5/2007 DSP & Quantization 52 Tanim Taher if PT == 1 otpt = fi([],T1,F); inpt = fi([],T1,F); elseif PT==2 otpt = single([]); inpt = single([]); w = single([]) elseif PT==3 otpt = double([]) inpt = double([]) w = double([]) end for i = 1:tot_pts if PT == 1 xn = fi(dpts(i),T1,F); elseif PT==2 xn = single(dpts(i)); elseif PT==3 xn = double(dpts(i)); end Fm(numel(oCK)+1)=xn; for jj=1:numel(oCK) if i < 2 G(jj,2)=0; else G(jj,2)=G(jj,1); end end for j=(numel(oCK)+1):-1:2 Fm(j-1) = Fm(j) - oCK(j-1) * G(j-1,2); G(j,1)=oCK(j-1) * Fm(j-1) + G(j-1,2); end G(1,1)=Fm(1); for n=1:numel(oCV) otpt(i)=0; otpt(i) = otpt(i) + G(n,1) * oCV(n); end end otpt = latc plot(xs(1:tot_pts), otpt(1:tot_pts)); function find_FP_output(handles) %Calculate the output waveform global a_no a_max b_no b_max c_A c_B PT WL FL A_C B_C F T1 global pts tot_pts xs dpts T % for slider otpt = []; inpt = []; % Set the data types if PT == 1 %Quantize the input values otpt = fi([],T1,F); inpt = fi([],T1,F); w = fi([],T1,F); elseif PT==2 otpt = single([]) inpt = single([]) ECE 491 6/5/2007 DSP & Quantization 55 Tanim Taher [k,v] = tf2latc(B_C,tmpA); [otpt,tmp_g] = latcfilt(k,v,inpt2); otpt = filter(B_C,tmpA, inpt2); for i = 1: tot_pts if PT == 1 otpt2(i) = fi(otpt(i),T1,F); elseif PT == 2 otpt2(i) = single(otpt(i)); else otpt2(i) = double(otpt(i)); end otpt2(i)=double(otpt(i)); end plot(xs(1:tot_pts), otpt2(1:tot_pts)) else find_lattice_op(handles); end end otpt = []; inpt = []; w = []; tmp_g = []; % --- Executes on button press in show_Output. function show_Output_Callback(hObject, eventdata, handles) pushbutton4_Callback(0, 0, handles) find_FP_output(handles) function sweep_freq(handles, plot_type) %This calculates the frequency response pushbutton4_Callback(0, 0, handles) global a_no a_max b_no b_max oCA oCB PT WL FL A_C B_C F T1 c_A = oCA; c_B = oCB; %only in this function due to later change T = 1 / str2num(get(handles.samplerate,'String')); freq = str2num(get(handles.f_start,'String')); fend = str2num(get(handles.f_end,'String')); f_step = (fend - freq) / 60; fn = 0; for fn = 1:60 frequ(fn) = freq; comp_sum_n = complex(0,0); w = 2 * pi * freq; %Vary the frequency for bk = 1:b_max %Calculate the complex value of numerator of transfer function with z=e^jw comp_sum_n = comp_sum_n + c_B(bk) * exp(complex(0,-(bk- 1) * w * T)); ECE 491 6/5/2007 DSP & Quantization 56 Tanim Taher end comp_sum_d = complex(1,0); for ak = 1:a_max %Calculate the complex value of denominator of transfer function with z=e^jw comp_sum_d = comp_sum_d + c_A(ak) * exp(complex(0,-ak * w * T)); end mag(fn) = abs(comp_sum_n / comp_sum_d); %Find the magnitude phs(fn) = angle(comp_sum_n / comp_sum_d) * (180 / pi); %Find the phase freq = freq + f_step; end %Plot the frequency response if strcmp(plot_type, 'magnitude') && get(handles.dBplot,'Value') == 0 plot(frequ(1:fn),mag(1:fn),'--.b'); elseif strcmp(plot_type, 'magnitude') && get(handles.dBplot,'Value') == 1 plot(frequ(1:fn),20 * log(mag(1:fn)),'--.b'); else plot(frequ(1:fn),phs(1:fn),'--*r'); end % --- Executes when user attempts to close figure1. function figure1_CloseRequestFcn(hObject, eventdata, handles) delete(hObject) % Hint: delete(hObject) closes the figure % --- Executes on button press in dBplot. function dBplot_Callback(hObject, eventdata, handles) % Hint: get(hObject,'Value') returns toggle state of dBplot function samplerate_Callback(hObject, eventdata, handles) % hObject handle to samplerate (see GCBO) % --- Executes during object creation, after setting all properties. function samplerate_CreateFcn(hObject, eventdata, handles) % hObject handle to samplerate (see GCBO) if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes during object creation, after setting all properties. function axes3_CreateFcn(hObject, eventdata, handles) % hObject handle to axes3 (see GCBO) % Hint: place code in OpeningFcn to populate axes3 % --- Executes on selection change in listAK. function listAK_Callback(hObject, eventdata, handles) ECE 491 6/5/2007 DSP & Quantization 57 Tanim Taher % hObject handle to listAK (see GCBO) % --- Executes during object creation, after setting all properties. function listAK_CreateFcn(hObject, eventdata, handles) % hObject handle to listAK (see GCBO) if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes on selection change in listBV. function listBV_Callback(hObject, eventdata, handles) % hObject handle to listBV (see GCBO) % --- Executes during object creation, after setting all properties. function listBV_CreateFcn(hObject, eventdata, handles) % hObject handle to listBV (see GCBO) if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes on button press in poles. function poles_Callback(hObject, eventdata, handles) %This plots the poles of the system pushbutton4_Callback(0, 0, handles) global oCA oCB a_max A = [1 oCA]'; % random polynomial poles = roots(A); % system poles for i = 1:(a_max) if abs(poles(i))>=1 %For unstable pole use red color plot(poles(i),'*r'); else plot(poles,'+b'); %For stable pole use blue color end hold on; end axis ([-5 5 -2 2]); grid on; hold off; % --- Executes on button press in latttype. function latttype_Callback(hObject, eventdata, handles)