Digital signal processing: modeling in MATLAB (RUB 367.00).

Youtube

A. I. Solonina S. M. Arbuzov Recommended by the UMO on education in the field of telecommunications as a teaching aid for students of higher educational institutions studying in the field of training for certified specialists 210400 "Telecommunications" St. Petersburg "BHV-Petersburg" 2008

Page 1

UDC 681.3.06 (075.8) BBK 32.973(ya73) S60 Solonina, A.I. S60 Digital signal processing. Modeling in MATLAB / A. I. Solonina, S. M. Arbuzov. - St. Petersburg: BHV-Petersburg, 2008. - 816 p.: ill. - (Tutorial) ISBN 978-5-9775-0259-7 Basic methods and algorithms are considered digital processing signals (DSP) and their computer modeling using the MATLAB system. The main modes of operation of the MATLAB system, matrix calculations, standard numerical methods and plotting are outlined. The specifics of representing DSP signals and systems in MATLAB are discussed in detail, linear discrete systems, synthesis of FIR and IIR filters, adaptive digital filtering, quantization, wavelets and modeling of these DSP objects and processes are described. software MATLAB, as well as a number of graphical programs included in the MATLAB extension package and designed to solve DSP problems using custom GUI without direct access to MATLAB software. For students and university teachers, as well as specialists in digital signal processing Reviewers: M. S. Kupriyanov, Doctor of Technical Sciences, Professor of the St. Petersburg State Electrotechnical University "LETI" V. A. Vargauzin, Ph.D. , associate professor of St. Petersburg State Polytechnic University UDC 681.3.06 (075.8) BBK 32.973 (ya73) Publication preparation group: Editor-in-Chief Deputy. Editor-in-Chief Head. Editorial Editor Computer layout Proofreader Series design Cover design Head. produced by Ekaterina Kondukova Tatyana Lapina Grigory Dobin Nina Sedykh Natalia Smirnova Victoria Piotrovskaya Igor Tsyrulnikov Elena Belyaeva Nikolay Tverskikh License ID No. 02429 dated 07.24.00. Signed for publication on 05/30/08. Format 70×1001/16. Printing of the background. Yes. ïå÷. ë. 65.78. Circulation 1500 ýêç. Order ¹ "Bal-Paratórðr", 194354, Sásíkò-Prátórðr, óë. Sanya, 5A. Sanitary and epidemiological conclusion for the product No. 77.99.60.953.Ä.002108.02.07 dated 02/28/2007. issued Federal service on supervision in the field of consumer rights protection and human well-being. Printed from ready-made transparencies at the State Unitary Enterprise "Printing House "Nauka" 199034, Sashaykt-Paparya, 9 th, 12 ISBN 978-5-9775-0259-7 © Solonina A. I., Arbuzov S. M., 2008 © Design, publishing house " BHV-Petersburg", 2008

Page 2

Contents PREFACE................................................... ........................................................ .............1 LIST OF ABBREVIATIONS IN RUSSIAN LANGUAGE.................................... ...............................3 LIST OF ABBREVIATIONS IN ENGLISH............... ....................................5 PART I. INTRODUCTION TO MATLAB .... ........................................................ ...................7 CHAPTER 1. INTRODUCTION TO THE MATLAB SYSTEM.................................... ...............................9 1.1. Accepted designations......................................................... ....................................10 1.2. Installing and launching MATLAB .................................................... ...........................11 1.3. MATLAB interface................................................... ........................................13 1.4. MATLAB help system................................................................... ..........................16 CHAPTER 2. DIRECT COMPUTATION MODE.................. ................................................24 2.1. Teams........................................................ ........................................................ .........25 2.2. Operators: assignment operator.................................................... ...................26 2.3. Constants........................................................ ........................................................ .......27 2.4. Variables........................................................ ........................................................ ....34 2.5. Functions........................................................ ........................................................ ..........37 2.6. Expressions........................................................ ........................................................ ......43 2.7. Symbols and functions of operations.................................................... ........................44 2.8. Working memory area Workspace: commands who, whos, clear ....................51 2.9. Saving data to disk: save, load commands ........................................................53 2.10. Creating your own folder and saving the path to it...................................53

Page 3

IV Contents CHAPTER 3. MATRIX LABORATORY.................................................... ........................55 3.1. Elements of matrices and access to them.................................................... ................56 3.2. Vector length and matrix size: functions length, size....................................58 3.3. Functions for generating standard matrices.................................................... ...............58 3.4. Formation of vectors and submatrices from a matrix...................................................61 3.5. Concatenation of submatrices and vectors into matrices.................................................63 3.6 . Copying matrices: repmat function .................................................... ...............64 3.7. Element-wise operations with matrices................................................................. ...............65 3.8. Operations with matrices in linear algebra problems....................................66 3.9. Operations with matrices in problems of mathematical statistics: functions max, min, sort, sum, prod, cumsum, diff, mean, std, var, cov, corrcoef ................... ........................................................ ........................................103 CHAPTER 4. TYPES OF ARRAYS.... ........................................................ ...............................111 4.1. Numeric arrays......................................................... ...........................................111 4.2. Non-numeric arrays................................................... ....................................114 4.3. Defining a data type: class function................................................... .......124 CHAPTER 5. GRAPHICS.................................... ........................................................ ..........127 5.1. 2D graphics: commands figure, hold on, hold off; subplot function................................................... ........................................................ ....127 5.2. Design of graphs: grid command, functions title, xlabel, ylabel, gtext, legend, xlim, ylim ..................... ........................................................ ..........129 5.3. Two-dimensional plots: functions plot, loglog, semilogx, semilogy, logspace, fplot................................... ........................................................ ....................130 5.4. Managing graph properties................................................................... ...................135 5.5. Special two-dimensional graphs: functions stem, stairs, polar, compass, bar, pie, hist ................... ........................................................ ........................139 5.6. 3D graphics................................................... ....................................145 5.7. Forming a mesh on the XOY plane: meshgrid function.......................145 5.8. Three-dimensional graphics: functions plot3, mesh, meshc, meshz, surf, surfl, surfc, contour3................................... ........................................................ ...............146 5.9. Managing the properties of three-dimensional graphs: colormap function; shading interp, colorbar commands .................................................... ...............................150 CHAPTER 6. NUMERICAL METHODS FOR SOLVING TYPICAL PROBLEMS............ ...............153 6.1. Operations with polynomials................................................................... ...............................153 6.2. Roots of the equation: fzero function .................................................... ...........................159

Page 4

Contents V 6.3. Approximation and interpolation................................................................. ...................161 6.4. Minimization of functions: functions fminbnd, fminsearch ....................................169 6.5. Numerical integration: functions trapz, cumtrapz, quad, quad1, dblquad ..................................... ........................................................ ...................172 6.6. Numerical integration of ordinary differential equations.................................................................... ........................................................ ..................174 CHAPTER 7. PROGRAMMING MODE.................................... ...................................187 7.1. User programs - M-files................................................... .............188 7.2. Structure of function files: functions nargin, nargout; commands type, global; return statement................................................ ...........................................188 7.3. Structure of script files: commands echo on, echo off....................................192 7.4. Development of programs in MATLAB.................................................... ...................194 7.5. Working with M-files................................................... ...........................................208 PART II. SIMULATION OF DSP USING MATLAB SOFTWARE.................................................... ........................................217 CHAPTER 8. DISCRETE SIGNALS.... ........................................................ ...................219 8.1. Representation of sequences................................................... ..........220 8.2. Random Sequences: functions rand, randn, xcorr, xcov.........244 CHAPTER 9. LINEAR DISCRETE SYSTEMS...................... ............................257 9.1. Modeling linear discrete systems in the time domain....257 9.2. Modeling of linear discrete systems in the z-domain....................................285 9.3. Modeling of linear discrete systems in the frequency domain........302 CHAPTER 10. STRUCTURES OF LINEAR DISCRETE SYSTEMS...................... ....313 10.1. Types of structures of FIR and IIR systems.................................................... .313 10.2. Description of the structures of FIR and IIR systems in the form of dfilt objects................318 10.3. MATLAB Functions for dfilt Objects .................................................... ..........325 10.4. Arrangement of links and scaling in dfilt objects: functions sos, scale .................................................... ........................................................ ..........332 CHAPTER 11. DISCRETE FOURIER TRANSFORM................................... ..........339 11.1. DFT calculation: functions fft, ifft, fftshift; external functions fft_e1, fft_e2 .................................................... ........................................................ ...........................342

Page 5

VI Contents 11.2. Convolution calculation using DFT: fftfilf function; external function iir_iir ................................................... ........................................................ .......356 11.3. Computing convolution with partitioning: fftfilt function .......................363 CHAPTER 12. RANDOM SIGNAL PROCESSING BY LINEAR DISCRETE SYSTEMS......... ........................................................ .......................367 12.1. Formation of random signals with a given probability distribution law.................................................... ....................................367 12.2. Formation of random signals with a given correlation function.................................................... ........................................................ ....................370 12.3. Nonparametric methods of spectral analysis: wvtool, psd functions ............................................................. ........................................................ ....371 12.4. Linear prediction: lpc function .................................................... ............381 12.5. Parametric methods of spectral analysis: functions pcov, arcov, pmcov, armcov, pburg, arburg, pyulear, aryule .....................387 CHAPTER 13. SYNTHESIS OF FIR FILTERS.... ........................................................ ...............391 13.1. Digital filters......................................................... ....................................391 13.2. Properties of FIR filters.................................................... ...............................393 13.3. Specifying requirements for frequency characteristics FIR filters......395 13.4. Synthesis of FIR filters using the window method: functions fir1, kaiserord; external functions check_low, check_high, check_pass, check_stop, plot_fir.......400 13.5. Synthesis of FIR filters using the best uniform (Chebyshev) approximation method: functions firpm, firpmord, firgr; external function MAG_fir................................................... ....................................417 13.6. Digital Hilbert converter.................................................... .......442 13.7. Digital differentiator................................................... ........................449 13.8. The following FIR filter structure: description as a dfilt object.................................................... ........................................................ ...................453 13.9. Analysis of digital filter characteristics: fvtool function......456 CHAPTER 14. SYNTHESIS OF IIR FILTERS............... ........................................................ .....457 14.1. Properties of IIR filters................................................................. ................................457 14.2. Specifying requirements for the frequency characteristics of IIR filters and the synthesis procedure.................................................... ........................................................ ......458 14.3. Synthesis of analog filters: functions butter, cheby1, cheby2, ellip, buttord, cheb1ord, cheb2ord, ellipord, freqs........................ ...............................460 14.4. Synthesis of IIR filters using the impulse response invariance method: functions impinvar, impulse................................................. ...................463

Page 6

Contents VII 14.5. Synthesis of IIR filters using the bilinear Z-transform method: functions bilinear, butter, cheby1, cheby2, ellip, buttord, cheb1ord, cheb2ord, ellipord...................... ........................................................ ........................469 14.6. IIR filter structure: description in the form of a dfilt object....................................480 14.7. Description of requirements for the frequency response of FIR and IIR filters in the form of fdesign objects .................................... ........................................................ ...............481 14.8. Synthesis of FIR and IIR filters in the form of dfilt objects based on fdesign objects .................................... ........................................................ .................486 CHAPTER 15. ADAPTIVE DIGITAL FILTERING................................... ...................493 15.1. Application of adaptation principles in DSP systems...................................493 15.2. Adaptive LMS algorithm: lms and nlms functions for adaptfilt objects.................................................... ........................................................ ...........................499 15.3. Adaptive RLS algorithm: rls function for the adaptfilt object ..................504 CHAPTER 16. QUANTIZATION IN FIXED-POINT DSP SYSTEMS...513 16.1. Quantization effects in the structure of a digital filter...................................515 16.2. Modeling the structure of a digital fixed-point filter.................................................... ........................................................ ........................542 16.3. Analysis of the characteristics of FIR and IIR filters with FT....................................566 16.4. Simulation of quantization in ADC.................................................... ............581 16.5. Calculating the response of FIR and IIR filters with FT: the filter function ......... 601 CHAPTER 17. MULTI-SPEED DSP SYSTEMS..................... ........................618 17.1. Single-shot interpolation systems................................................................... ..........619 17.2. Simulation of single interpolation in MATLAB: functions interp, upfirdn .................................................... ........................................................ 624 17.3. Single-use decimation systems................................................................. ...............629 17.4. Simulation of single decimation in MATLAB: functions decimate, upfirdn .................................................... ........................................635 17.5. Single-shot resampling systems......... ...................................638 17.6. Simulation of single resampling in MATLAB: functions resample, upfirdn................................................... ........................................639 17.7. Description of the polyphase structure of interpolation and decimation systems in the form of mfilt objects.................................................. ................................642 CHAPTER 18. WAVELET SIGNAL PROCESSING......... ........................................651 18.1. Basic concepts of wavelet analysis.................................................... ...............652 18.2. Wavelets in MATLAB: functions wavemngr, waveinfo, wavefun, centfrq ................................................... ........................................................ ...............................655

Page 7

VIII Contents 18.3. Continuous wavelet transform: cwt function ....................................667 18.4. Scaling filters: functions dbwavf, symwavf, coifwavf, biorwavf, rbiowavf...................................... ........................................................ ............671 18.5. Decomposition and recovery filters: functions orthfilt, wfilters, qmf, dwt, iwdt................................... ........................................................ ...................672 18.6. Multilevel wavelet analysis: functions wavedec, waverec, appcoef, detcoef, swt, iswt ......................... ........................................................ ....680 18.7. Wavelet packets: functions wpdec, wpcoef, wprec, wentropy, besttree ........684 CHAPTER 19. INTERACTION WITH EXTERNAL SIGNAL SOURCES........690 19.1. Data formats compatible with signal analysis tools in MATLAB.................................................... ........................................................ ...................690 19.2. Using ready-made signals: wnoise function ....................................692 19.3. Import external files: wavread function...................................................697 19.4. Playing sound: functions sound, soundsc, wavplay.......................700 19.5. Record sound files: functions wavrecord, wavwrite...................................702 PART III. MODELING DSP USING GUI...............................................705 CHAPTER 20. DESIGNING DIGITAL FILTERS USING GUI FDATOOL... ........................................................ ........................................................ ....707 20.1. Synthesis of digital filters................................................................... ...............................708 20.2. Input parameters digital filters........................................................ ..710 20.3. Examples of digital filter synthesis.................................................................... .........715 20.4. Selecting a digital filter structure.............................................................. ..........721 20.5. Analysis of digital filters................................................................... ........................722 20.6. Synthesis of digital Hilbert converters.................................................724 20.7. Synthesis of digital differentiators.................................................... .........726 20.8. Saving digital filters for the duration of a session in the FDATool GUI.........728 20.9. Exporting digital filters as dfilt objects ........................................................729 10.20. Importing digital filters as dfilt objects.................................................732 20.11. Modeling the structure of digital fixed-point filters.................................................... ............................................733 CHAPTER 21. DIGITAL MODELING FILTERING USING GUI SPTOOL .................................................... ........................................................ ...............742 21.1. Synthesis of digital filters................................................................... ...............................743 21.2. Input parameters of digital filters................................................................... ....747

Page 8

Contents IX 21.3. Examples of digital filter synthesis.................................................................... .........750 21.4. Analysis of digital filters................................................................... ........................755 21.5. Importing an input signal................................................................... ...............................756 21.6. Simulation of digital filtering................................................................... .......760 21.7. Time domain signal analysis................................................................... .......761 21.8. Signal analysis in the frequency domain................................................................... ..........763 21.9. Exporting data from GUI SPTool .................................................... ....................766 21.10. Exiting GUI SPTool.................................................... ...................................772 CHAPTER 22. MODELING WAVELET TRANSFORMATIONS USING THE WAVELET TOOLBOX GUI PACKAGE.. ........................................................ 773 22.1. Viewing Wavelets................................................... ....................................774 22.2. One-dimensional discrete wavelet analysis.................................................... .....775 22.3. One-dimensional packet wavelet analysis.................................................... .........782 22.4. Real and complex one-dimensional continuous wavelet analysis.................................................................... ........................................................ ..........785 22.5. Removing noise from a stationary random one-dimensional signal.....787 22.6. Estimation of distribution density......................................................... ................789 22.7. Regression Estimation................................................................ ........................................791 22.8. Selection of wavelet coefficients.................................................... ....................793 REFERENCES.................................... ........................................................ ..................795 INDEX.................................... ........................................................ .........798

Page 9

Preface Modern trends in the field of telecommunications are largely associated with the development of digital equipment and software product, and this radically changes the nature of the work of engineers and scientists - it is increasingly reduced to computer modeling. A feature of digital signal processing (DSP) devices is that the software parts of these devices are created directly in the process of computer modeling, so mastering its modern technologies comes to the fore. Such technologies, of course, include the MATLAB software environment (system), created by The Math Works, Inc. and intended for computer modeling in various fields of science and technology. IN last years The discipline "Digital Signal Processing" and its modifications are included in the general educational standards of Russian universities. It is remarkable that domestic and translated books on DSP theory, its applications and the implementation of DSP algorithms are increasingly being published - “the process has begun.” Often theoretical sections are supported by calculation examples in MATLAB. Books of an applied nature began to appear, which are not fragmentary, but entirely devoted to modeling in MATLAB, for example, the book by R. Gonzalez et al. “Digital image processing in MATLAB”, publishing house “Technosphere”, 2006. Nevertheless, today there is a demand for literature DSP modeling in MATLAB far exceeds the proposal. This is a very broad topic, and this book covers basic DSP methods and algorithms, theoretical basis which are presented in many sources, including textbook by the same authors, "Fundamentals of Digital Signal Processing", publishing house "BHV-Petersburg", 2005. (In the near future, its third re-release is in the works.) Let's immediately say that MATLAB is an immense system, and even in this limited area, the authors in no way claim to provide an exhaustive description MATLAB capabilities. The nuances, subtleties and details, as well as the variety of tools offered to solve problems, can only be understood in practice, using the powerful MATLAB help system. The methodology for teaching computer modeling is special. In fact, it comes down to self-education - independent expansion of knowledge after acquiring

Page 10

2 Preface to initial skills, mastering standard techniques and determining the search vector in the vast MATLAB system, which the authors tried to systematize and describe in application to DSP modeling in MATLAB. The book includes many examples - with their help, “without further ado,” you can quickly master modeling technology. To make it easier for novice users and the self-contained nature of the book, the first part is included, which covers the basics of working in MATLAB. This book can be useful for all engineers and technical workers with an interest in the field of DSP, but it is primarily aimed at undergraduates, graduate students and university teachers and can be recommended, in particular, for the following compulsory disciplines:  Microprocessors and digital signal processing (specialty 210405);  Digital signal processing and signal processors in mobile communication systems (specialty 210402);  Digital processing of audio and video signals (specialty 210312). It is assumed that readers are familiar with the basics of DSP theory and programming in any language high level. The book provides only brief theoretical information on the relevant sections of the DSP. The contents of the book include 22 chapters, which are thematically divided into three parts: 1. Introduction to MATLAB. 2. Modeling of DSP using MATLAB software. 3. DSP modeling using GUI. Authors of parts and chapters:  A. I. Solonina - part I; in part II chapters 8-11, 13, 14, 16, 17; in part III, chapters 20, 21.  S. M. Arbuzov - in part II, chapters 12, 15, 18, 19; in part III - chapter 22. Alla Ivanovna Solonina, prof., Ph.D., and Sergey Mikhailovich Arbuzov, associate professor, Ph.D., teach at the Department of Digital Signal Processing of the St. Petersburg State University of Telecommunications. prof. M.A. Bonch-Bruevich, headed by Prof. Artur Abramovich, Doctor of Technical Sciences. Lanne, Please send all suggestions and comments that will be accepted by the authors with gratitude to the publishing house "BHV-Petersburg" by email: [email protected].

Generating signals in Simulink naturally has its own characteristics. Let's look at them.

Let's take two blocks from the Simulink block library: Sine Wave(from section Sources) And Scope(from section Sinks). By connecting them, we get a simple circuit (Fig. 4).

Fig.4. Circuit for generating and displaying a sinusoidal signal

Then, by double-clicking on the oscilloscope (oscilloscope) block, we activate a window simulating the oscilloscope screen and launch the model (button Start simulation). As a result, we obtain an image of a sinusoid segment (Fig. 5).

Fig.5. Displaying a sine wave segment on the oscilloscope screen

As you can see, generating a harmonic signal in the Simulink environment is even easier than in the MATLAB environment. However, this first impression is very misleading. Indeed, it is also important to be able to control the parameters of the harmonic signal. The fact that the amplitude of the harmonic signal turned out to be equal to unity is simply “lucky” for us. Indeed, by default the amplitude of the generated signal is set to unity. However, we do not yet control the frequency, initial phase and duration of the signal.

Double click on the block Sine Wave– as a result, the parameter settings window will appear (Fig. 6). By clicking the button Help, we get instructions for this block, the essence of which boils down to the fact that an operation is performed in this block

From the above formula and the inscriptions in Fig. 6, the meaning of four variables becomes clear: amplitude, angular frequency, initial phase and constant component. The meaning of the “time” variable remains encrypted for now.

Stopping there important issue, let us note the difference between the concepts of “time” and “model time”. Thus, generating a signal segment with a duration of 1 s (model time) can last significantly longer than a shorter period of time, for example, 0.1 s ( real time). The generation speed depends on the amount of calculations, the speed of the computer, and the chosen “solver”, i.e. modeling algorithm, etc. By the way, the opposite effect is quite possible - for a complex algorithm, the procedure for modeling a signal segment with a duration of 0.1 s can last for several seconds.

The signal can be generated in two types: continuous time-based and discrete sample-based. Accordingly, to simulate the work continuous systems It is recommended to use the continuous type time-based, and for modeling the operation of discrete systems - discrete type sample-based.

If type is set time-based, then the parameter Sample time can take values:

– 0 (default) – the unit operates in continuous mode;

– > 0 – the block operates in discrete mode;

– -1 – the block inherits the same mode as the receiving block.

As stated in Help, operation in continuous mode can lead to large generation errors over large periods of model time.

Fig.6. Block parameters settings window Sine Wave

Operation in discrete mode causes the block to behave as if a block were connected to the output of a continuous oscillator Zero-Order Hold(from the Discrete section). Indeed, by assembling two circuits (Fig. 7) and setting the parameter value in both cases Sample time equal to 0.5 (block settings window Zero-Order Hold shown in Fig. 8), we obtain identical results (Fig. 9).

Fig.7. Inserting a sine wave into a circuit that generates and displays a sine wave

block Zero-Order Hold

Fig.8. Block settings window Zero-Order Hold

Fig.9. Identicality of the result of the circuits shown in Fig. 7

So the block Zero-Order Hold can be interpreted as a “sampler”, i.e. the part of the analog-to-digital converter (ADC) responsible for sampling the signal. Sometimes a block Zero-Order Hold called ADC. In our opinion, this is not correct, since the sampled signal in a “genuine” ADC is also subject to level quantization. In the block Zero-Order Hold, however, no quantization is performed.

A few notes about graphing methods. Beyond the block Scope, the graph can also be constructed using the block X-Y-Graph, to the upper entrance X which you need to submit a sequence of moments in time using a block Clock(clock), and to the lower entrance Y– values of the generated signal (Fig. 10).

Fig. 10. Application for block plotting X-Y-Graph

As a result, pre-configured (in the corresponding settings window, the limit values of the argument and function are set, and the value of the parameter is also indicated Sample time) the plotter will produce the graph shown in Fig. 11 if for the block X-Y-Graph given Sample time=-1(i.e. the sampling period is inherited).

Fig. 11. Displaying a sine wave as a block X-Y-Graph

The graph will be slightly different (Fig. 12) if for the block X-Y-Graph given Sample time=0.5.

Fig. 12. Displaying a sine wave as a block X-Y-Graph, Sample time=0.5

Another way to build graphs. Arrays of samples of time instants and corresponding signal values can be used using the block To Workspace export from the Simulink environment to the MATLAB environment (Fig. 13).

Fig. 13. Application of a block To Workspace

In this case, as practice shows, it is best to set the format array for exported data (Fig. 14).

Fig. 14. Setting the format array for block To Workspace

Further plotting in the MATLAB environment using the plot(x,y) command is not difficult (Fig. 15).

Fig. 15. Plotting data exported

using a block To Workspace

Let us summarize the results obtained.

Signal type time-based when the generation unit operates in continuous time mode, it has the form of a smooth function of time, and in discrete time mode, it has the form of a step signal, such as if the block was connected to the output of the smooth signal generator Zero-Order Hold, which is a sample-storage type sampler.

In other words, by setting the discrete time mode, we avoid the need to use the block Zero-Order Hold.

Now let's generate Simulink a segment of a discrete harmonic signal with the same parameters that were specified in MATLAB: amplitude 1, frequency 100 Hz, sampling frequency 1000 Hz, initial phase π/2, number of samples 20.

We reassemble the circuit from a generator and an oscilloscope. In the generator settings mask window, we specify the required numerical values of the parameters, set the type time-based and assign a value Sample time= 0.001 (Fig. 16).

Fig. 16. Generator settings window

After starting the model, we get on the oscilloscope screen a completely different picture than we expected (Fig. 17).

Fig. 17. Oscilloscope result

The reason is simple - you still need to configure the modeling parameters: set the start and end of the model time (in our case, these are 0 and 0.02 s, respectively), and also select the modeling algorithm (the “solver” type). Figure 18 shows the simulation parameters settings window, which is activated when selecting a menu item Simulation/Simulation parameters. Often these parameters are configured automatically, but knowledge of them is necessary to isolate individual parts of the function.

Fig. 18. Modeling parameters settings window

In addition, let’s configure the oscilloscope parameters by clicking on the button Parameters on the window Scope(Fig. 19 a, b).

Fig. 19. Setting up oscilloscope parameters:

a) setting up the General tab; b) setting up the Data history tab

After starting the model, an image will appear on the oscilloscope screen (Fig. 20).

Fig.20. The result of the oscilloscope after the settings have been made

Since the oscilloscope parameters were set so that two-dimensional array ScopeData argument and function values, using commands

>> y1=ScopeData(:,1);

>> y2=ScopeData(:,2);

>> plot(y1,y2)

you can plot the generated function using MATLAB (Fig. 21).

Fig.21. Graph of the generated function using MATLAB

Comparing Fig. 21 and Fig. 2, we notice only one difference - when modeling in Simulink, 21 points were generated, while in MATLAB 20 points were generated. The reason for the difference is simple: over the model time interval at sampling rate F s located TF s +1 points in time for which the signal will be generated. Obviously, this circumstance can be easily taken into account, achieving complete coincidence of simulation results in the MATLAB and Simulink environments.

Concluding remarks on 1.1.

Sampling of analog signals is the first step towards solving the problem of interfacing analog devices and systems with discrete ones.

Modeling of discrete signals can be done either in MATLAB or in Simulink. Maybe sharing these environments, which increases the flexibility of the toolkit.

There are three ways to generate signals in MATLAB:

in dialog mode (sequence of commands in the command window);

in automatic mode, by creating and running an m-script;

in automatic mode, by creating and calling an m-function.

Modeling signals in Simulink is convenient due to its clarity, but it requires familiar skills in setting the parameters of the blocks from which the model is constructed.

An important feature of simulation in Simulink is the obvious difference between the concepts of “real time” and “model time”. Real time is the time required to carry out calculations (simulations). Model time is the duration of the simulated process.

Model time can be continuous ( time-based) and discrete ( sample-based). Continuous time is recommended to be used when modeling continuous (analog) systems, discrete time - discrete (digital). Monitoring of simulation results in Simulink can be done both by plotting graphs and by printing the values of numerical arrays.

It is more convenient to construct graphs in the Simulink environment. It is more convenient to analyze numerical data by exporting it to the MATLAB workspace.

MINISTRY OF EDUCATION AND SCIENCE OF THE RUSSIAN FEDERATION
GOUVPO "SIBERIAN STATE GEODETIC ACADEMY"

T.V. Dashkovskaya

DIGITAL SIGNAL PROCESSING IN MATLAB ENVIRONMENT
Part 1
Approved by the Academy's Editorial and Publishing Council
as a laboratory workshop for students studying
in the direction 200200 “Optotechnics”, specialty 200203
"Optical-electronic devices and systems"

Novosibirsk
SSGA
2010

UDC 528.7:004
D217
Reviewers: Candidate of Technical Sciences, Associate Professor, NSTU G.A. Syretsky
Doctor of Technical Sciences, Professor, SSGA V.V. Malinin
Dashkovskaya T.V.
D217 Digital signal processing in the matlab environment [Text]: laboratory
workshop at 2 p.m. Part 1 / T.V. Dashkovskaya. - Novosibirsk: SGGA, 2010. - 66 p.

ISBN 978-5-87693-398-0
The laboratory workshop in the disciplines “Image Theory in
optical-electronic systems", "Digital image processing" (4th year,
bachelors), " Computer modelling and design" (3rd year),
“Special chapters of COI” (6th year, masters) is intended for students
students in the direction 200200 “Optotechnics”, specialty 200203
"Optical-electronic devices and systems."
The laboratory workshop consists of two parts. The first part includes
laboratory works related to the presentation of information and elements
programming in the system, modeling and signal processing with
using the Signal Processing ToolBoxes system package, part two
contains laboratory work on digital image processing with
using the system package Image Processing ToolBoxes.
Responsible editor: candidate of technical sciences, associate professor, SSGA E.V.
Gritskevich
Published by decision of the SSGA Editorial and Publishing Council
UDC 528.7:004
ISBN 978-5-87693-399-7 (part 1)
ISBN 978-5-87693-398-0

GOU VPO "Siberian State
Geodetic Academy" (SSGA), 2010

CONTENT
Laboratory work No. 1. General principles operation of the MATLAB system .... 4
Laboratory work No. 2. Programming elements.................................. 15
Laboratory work No. 3. Plotting graphs.................................................... 26
Laboratory work No. 4. Signal generation.................................................... .. 36
Laboratory work No. 5. Representation of linear systems.................................... 48
List of recommended literature for in-depth study of the material
..................................................................................................................... 60

LABORATORY WORK No. 1. GENERAL PRINCIPLES OF WORK
MATLAB SYSTEMS
The purpose of the work is to familiarize yourself with the MatLab system, the rules for creating
numerical arrays and acquiring practical skills in using
system tools for working with them.
1. BRIEF THEORETICAL INFORMATION
The name of the MatLab system comes from the words Matrix Laboratory
(matrix laboratory). The package is focused on processing data sets.
The MatLab system interface includes the following panels:
Command Window, where all calculations and
operations;
Launch Pad, where you can access various
ToolBox modules;
Workspace, where the current set is displayed
variables entered by the user in the command window;
Current Directory, where you can set the current
catalog;
Command History, where the commands you type are stored
by users.
Matrix MatLab system stands out from other systems in that it
Operators and functions have operands in the form of vectors and matrices. Even the operand
in the form of a single number is considered as a matrix of size 1 1. Since
operations with matrices can be either element-wise or matrix-based, then in
elementwise operators add a dot. For example, the symbols dot,
asterisk ‘*’ define element-wise multiplication of arrays, asterisk symbol
‘*’ - matrix multiplication (Table 1). Any command set must
end with a key press . The action performed by the function is
applies to all array elements passed in the input list
arguments.
As an example, we can give the result of calculating the sine:
>> y = sin()
y=
0.3894 0.7174 0.9320 0.9996 0.9093
The semicolon (;) at the end of a statement or function suppresses output
result on screen:
>> y = sin();
reference Information
Get background information possible with the following operators:
Syntax
helpwin
- information about the sections and functions of the MatLab system;
helpdesk
- general information about the MatLab system;

Doc<имя_функции>
- displaying a description of the function in the Help window;
help<имя_функции> - brief information about the function;
type<имя _функции>- output of the text m - function file;
demo
- command to call test cases.
Operators and functions
The complete list of system operators is displayed by the operator:
help ops
In table 1 shows a list of arithmetic operators with their syntax
applications.
Table 1.1 List of arithmetic operators
Function

Name

Plus
uplus
minus
uminus
mtimes
times
mpower
power
mldivide
mrdivide
lvidide
rdivid

Plus
Unary plus
Minus
Unary minus
Matrix multiplication
Element-wise array multiplication
Matrix exponentiation
Element-wise array exponentiation
Reverse (right to left) matrix division
Dividing matrices from left to right
Element-wise division of arrays from right to left
Elementwise division of arrays from left to right

Operator Syntax
+
+
-
-
*
.*
^
.^
\
.\
./

M1+M2
+M
M1-M2
-M
M1*M2
M1.*M2
M1^x
M1.^x
M1\M2
M1/M2
M1.\M2
M1./M2

To enter a comment, use the percent sign - %.
System variables and constants
MatLab objects include a number of system variables and constants,
whose values are set by the system when it boots or automatically
are formed during the calculation process. Listed below are some of these
objects:
Ans - used to write the result if the operator is not used
assignments.
For example:
>>cos()
ans =
1.0000 0.5403 -0.4161 -0.9900 -0.6536 0.2837 0.9602
pi is a number (the ratio of the circumference of a circle to its diameter).
For example:
>> pi
ans =
Angle brackets are used in command syntax for conditional
argument notation.

3.1416
Creating Vectors and Matrices
Let's look at some ways to form numeric arrays. With them
creation you can use:
Square brackets;
Special design j:i:k;
Concatenation (union);
Special matrix functions.
To create a row vector, use square brackets with
listing the elements of a line separated by spaces or commas and a special
construction j:i:k indicating the initial value of the vector - j, step - i and
the final value of the vector - k separated by a colon (if the step value is 1,
it may not be specified).
To create a column vector, the elements of the vector are listed separated by a dot
with a comma in square brackets or the previously obtained string vector is transposed. To perform the transposition operation, a single
quotation mark ("), which is placed after the identifier defining
transposeable structure. For complex matrices, transpose
is complemented by matrix conjugation. Single quote period (.")
used to transpose an array without a conjugation operation for
complex matrices.
To create a matrix, you can use the following input methods
elements in square brackets:
1. On lines separated by semicolons;
2. By columns specified in square brackets;
3. Line by line in interactive mode.
Task 1. Create a row vector, column vector and matrices:
>> x = % creating a row vector
x=
12345
% result
>> x = 1:2:10% 1 - initial vector value, 2 - step, 10 - final
meaning
x=
13579
% result
>> x = 1:10
x=
1 2 3 4 5 6 7 8 9 10% result
>> x = % creating a column vector
x=
1
2
3
4

>> a = 9:2:18% row vector with initial value 9, final value 18,
in steps of 2
a=
9 11 13 15 17
>> a1 = a"
% transpose row vector a to column vector a1
a1 =
9
11
13
15
17
>> x = ;
% creating a matrix using method 1)
>> x = [ ]; % creating a matrix using method 2)
>> x = ;
Executing each of the last three commands will create
matrices:

1 2 3
4 5 6
7 8 9

Concatenation (merging) of arrays
Using the concatenation operation, you can form new arrays from
previously created - vectors, matrices, using these arrays as their
elements. Arrays can be combined horizontally and vertically.
In horizontal concatenation, as an array separator in
square brackets use a comma or space, for example, if B and A -
matrix, then M = [A, B] - horizontal concatenation of matrices A and B. A and B
must have the same number of lines. Horizontal concatenation can
can be applied to any number of matrices within the same brackets.
In vertical concatenation as an array separator
a semicolon is used in square brackets, for example if C and D are
matrices, then M = is the vertical concatenation of matrices C and D. C and D
must have the same number of columns. Vertical concatenation can
can be applied to any number of matrices within the same brackets.
Vertical and horizontal concatenation can be applied in one
operations.
Task 2. Create a matrix using vertical and horizontal
concatenation:
>> a = % creating a column vector of three elements
a=
1
2
3

>> A = % horizontal concatenation of three column vectors a
A=
111
222
333
>> b = ; ] % creates a row vector of three elements
>> B = % vertical concatenation of three row vectors b
B=
148
148
148
>> A1 = ; % creation of matrices A1 and A2 of dimension 2 2
>> A2 = ;
>> A3 = ;
% creates a row vector of four elements
>> H = % horizontal concatenation and vertical
% concatenation
H=
0056
0118
1869
Special Matrix Functions
Below are some functions for creating special matrices
kind.
The zeros function creates a matrix filled with zeros (null matrix):
zeros(m,n) - specifying a matrix of m n zeros;
zeros(n) - task square matrix of n n zeros.
An example of creating a zero matrix from two rows and three columns:
>> x = zeros(2,3)
x=
000
000
The ones function creates a matrix filled with ones (one
matrix):
ones(m,n) - creating a unit matrix of size m n;
ones(n) - specifying a square matrix of n n ones.
Task 3. Create a square identity matrix of dimension 2 2:
>> ones(2)
ans =
11
11
Function
randn
creates
matrix,
completed
normally distributed random numbers.

Randn (m,n) - specifying a matrix of m n normally distributed
random numbers.
randn (n) - specifying a square matrix n n.
Function
rand
creates
matrix,
completed
uniformly distributed random numbers.
rand (m,n) - specifying a matrix of m n uniformly distributed
random numbers.
Task 4: Create a column vector using vertical concatenation
using rand and randn functions:
>> x1 =
x1 =
0.9501
0.2311
-0.4326
-1.6656
Task 5. Create a matrix using horizontal concatenation with
using the ones and zeros functions:
>> x1 =
x1 =
10
10
The repmat() function creates a matrix by copying source array given
number of times vertically and horizontally.
B = repmat(A,M,N) - the function creates a matrix B consisting of M copies
A vertically and N copies of A horizontally, that is, M N copies of the array A
(if A is a number, the function forms a matrix of size M N with the value
elements equal to A).
Task 6. Generate a matrix using a row vector a from
three elements:
>> a = ;
>> A = repmat(a,2,1)
A=
159
159
Task 7. Form a matrix of dimension 2
which is equal to ten:
>> repmat(10,2,3)
ans =
10 10 10
10 10 10

3, all elements

Array indexing
Array elements have two properties: serial number
(index) in the array and the actual value. Numbering of elements in the system
MatLab starts with one. To specify the indices of array elements
parentheses are used (an array indexing error is generated in
if the element index is less than one or greater than the size
array).
Task 8. Set a row vector of four elements and change the third
element to value 8:
>> a = ;
>> a(3) = 8
a=
1289
To indicate an element in a matrix, the indices are the row number and the number
columns separated by commas.
Task 9. Change the value of an element of the random number matrix S,
located in the second row and in the fourth column, at 1:
>> S = rand(4)
S=
0.9501 0.8913 0.8214 0.9218
0.2311 0.7621 0.4447 0.7382
0.6068 0.4565 0.6154 0.1763
0.4860 0.0185 0.7919 0.4057
>> S(2,4) = 1
S=
0.9501 0.8913 0.8214 0.9218
0.2311 0.7621 0.4447 1.0000
0.6068 0.4565 0.6154 0.1763
0.4860 0.0185 0.7919 0.4057
To indicate a block of array elements, use the colon ":" symbol.
Zero in a matrix F of integers with dimensions 4 4 elements,
located in the second and third row and in the first and second columns:
>>F=
F=
1479
5983
9653
7692
>> F(2:3,1:2) = 0
F=
1479
0083
0053

7692
If you need to change the value of an entire column or row, then the numbers
indicating a range of values are not specified and one colon remains.
Task 10. Reset the third and fourth columns from the previous one
example:
>> F(:,3:4) = 0
F=
1400
0000
0000
7600
Empty square brackets remove information from the indexed
structures.
A(m,:) = - removes row m from matrix A.
A(:,n) = - removes column n from matrix A.
Service functions
Below are some of the features you need when working with
arrays:
= size(<идентификатор_массива>) - returns the size of the array,
where M is the number of lines; N - number of columns.
Task 11. Determine the dimension of the unit matrix:
>> s = ones(2,3);
>> = size(s)
M=
2
N=
3
max(<идентификатор_массива>) and min(<идентификатор_массива>) -
computes a string vector containing the maximum or minimum
elements in each column of the matrix.
Task 12. Determine maximum values random matrices
numbers:
>> v = rand(3)
v=
0.9501 0.4860 0.4565
0.2311 0.8913 0.0185
0.6068 0.7621 0.8214
>> ma = max(v)
ma =
0.9501 0.8913 0.8214
% maximum values of each column.
To find an extreme value in a matrix, you need to transform it into
vector. To do this, you can use the reshape (X,M,N) function, where X is

Transformable matrix, M - vertical dimension of the matrix (number of rows),
N - horizontal dimension (number of columns).
Task 13. Convert the matrix from the previous example into a vector row and find extreme values.
>> = size(v);
>>k = reshape(v,1,M*N)
k=
0.9501 0.2311 0.6068 0.4860 0.8913 0.7621 0.4565 0.0185 0.8214
>> max(k)
ans =
0.9501
>>min(k)
ans =
0.0185
There is an easier way to find the extreme value in
array with dimension greater than one - represent this array
one-dimensional, using indexing of all values:
>>max(v(:))
ans =
0.9501
length(<идентификатор_массива>) - determination of the vector length; For
matrix, this is equivalent to executing the functions max(size(X)).
Task 14. Determine the length of a given vector:
>> c = ;
>>length(c)
ans =
7% vector length s
A single quote (") is used to create a string constant,
For example:
>> a = "Enter matrix";
Single quote is also used to perform the operation
transpose (see page 5).

2. TASKS FOR INDEPENDENT SOLUTION
1. Create a row vector: initial element is - , final, step
equals 0.1. Transpose a row into a column.
2. Create three row vectors of 5 elements fi = , where
n = 5 for x = 2, 3, 4. Combine these rows into a matrix A(3 × 5).
3. Create three column vectors from 5 elements of arithmetic
progression. An element of an arithmetic progression is calculated using the formula:
an = an-1 + d,
where an-1 is the previous element; an - subsequent.
Five vector elements are formed starting from specifying the first element
a and c using arithmetic progression step d to specify
following elements:
For the first column vector a = 2; d = 1;
For the second column vector a = 7; d = 2;
For the third column vector a = 10;

The purpose of the laboratory work is to study the methodology for developing programs for complex types of digital signal processing, including a combination of key operations (FFT, correlation, spline approximation and resampling).

Assignment for work

There is a set of experimental data in the form of a numeric array. It is required to design a digital data processing program in the internal MATLAB language that implements an accurate determination of the number of signal periods and frequency in the time domain using several key DSP operations: FFT, correlation, spline approximation and resampling.

Theoretical basis

One of the most important tasks of digital processing of noisy signals is the detection of an informative signal in a data stream distorted by noise and interference, and the determination of its parameters. For this, various methods are used, such as time filtering (accumulation), optimal frequency filtering, direct and inverse Fourier transform, and correlation processing. Each of these operations allows you to perform transformations of the original signal, for example, the transition of a signal from the time domain to the frequency domain or vice versa, and this reduces the noise level in the processed signal. In tasks of detecting and determining the parameters of noisy signals, enhancing the effect of noise suppression and increasing the accuracy of determining signal parameters can be achieved using several digital processing methods in combination. An example would be the task of processing the echo signal of an NMR spectrometer

The result of the FFT of a sampled NMR spectrometer echo signal of a certain frequency is the number of signal periods in the time window. If the sampling frequency or the time discrete interval when measuring a signal is known, then the frequency of the measured signal can be determined by the number of periods in the time window. The accuracy of determining the frequency in the spectrum of the input signal is quite definite and depends on the number of periods p of the signal. If the number of periods is an integer, then the frequency can be found absolutely accurately using the FFT (in the absence of signal noise). If the number of periods is not an integer, then an error in determining the frequency appears. The maximum error value is 1/p. In some practically important cases, for example, when processing echo signals (Fig. 1) of pulsed NMR spectrometers, the number of periods of the analyzed signal in a time window is fundamentally limited to about 10. In this case, the error in determining the frequency using FFT reaches 1/10, i.e. 10%.

The influence of noise in the recorded signal in the time domain can be significantly reduced by repeating the experiment multiple times and synchronously accumulating echo signals (see Fig. 1).

Rice. 1. Initial (A) and accumulated signal (B) in iron borate.

However, as shown in , increasing the number of accumulations makes it possible to improve the signal-to-noise ratio without distorting the shape and reducing the amplitude of the accumulated echo signal only to a certain limit. After this limit, accumulation no longer brings a noticeable improvement in quality. In particular, with a long accumulation time, the echo signal begins to be affected by gradual changes in the parameters of the devices included in the pulsed NMR spectrometer. When limiting the time for analyzing substances, the number of possible signal accumulations should be limited or absent altogether. The problem of processing echo signals in conditions of significant noise also arises when there is a need to study substances of low concentration. Therefore, the task of performing frequency analysis of noisy echo signals is relevant.

Direct use of the FFT for a noisy signal does not allow one to obtain an accurate value for the number of periods in the time window and the frequency of the measured signal in the case where the number of periods is not an integer, when the analyzed signal occupies only part of the time domain, is amplitude modulated and is noisy. In practically important cases of frequency analysis of signals from NMR spectrometers, the shape and initial phase of the echo signals are known. This makes it possible to create reference signals corresponding to the expected echo signal in shape and initial phase, and make their correlation comparison. The correlation coefficient of the echo signal with the reference signal will be equal to one, if the frequencies of the echo signal and the reference signal are equal and the echo signal is not noisy. Therefore, in the absence of noise, the frequency of the echo signal can be found by making a correlation comparison with the reference signals; the frequency of the reference signals can be selected until the condition is met when the correlation coefficient is equal to unity. However, the correlation coefficient decreases both when the frequencies of the echo and reference signals differ, and when the frequencies coincide, but due to the presence of noise. Therefore, it is impossible to determine the frequency of a noisy signal in this way.

Improvement of the algorithm for determining the signal frequency is achieved by combining the positive qualities of the correlation approach and FFT in order to increase the accuracy of determining the frequency of the sampled signal in conditions when the processed signal in the time domain is noisy, the number of periods is not an integer, the signal is amplitude modulated, and the number of signal periods is small.

The idea of the proposed algorithm, described in, is that in a small vicinity of the expected signal frequency or number of periods (the approximate value of the signal frequency can be found using the fast Fourier transform), correlation coefficients are calculated with several reference signals in a certain vicinity of the approximate value, then, using spline interpolation and resampling, a function is constructed that expresses the dependence of the correlation coefficient on the frequency of the standards and the maximum of this function is found; the position of the maximum is used to determine the frequency of the standard signal. The function constructed in this way has the form of a parabola with a clearly defined maximum (see Fig. 2) both in the case of a non-noisy and noisy signal, which makes it possible to determine the frequency of the echo signal more accurately than the FFT allows. In the presence of noise, the shape of the function is preserved, only the absolute value of the maximum decreases.

Rice. Fig. 2. Dependence of the correlation coefficient on frequency in the absence of noise (A) and with a signal-to-noise ratio of 1/3 (B). The exact frequency value is 1010. The fitting curves are constructed using the spaps spline fitting function in MATLAB.

The closer the initial approximation to the true frequency is, the higher the accuracy of determining the signal frequency. Therefore, an iterative calculation is used; at each iteration stage, the refined frequency value obtained at the previous stage is used as an initial approximation. The frequency determined using the FFT is taken as a first approximation.

The number of iterations to calculate the frequency with a given accuracy using the described algorithm depends on how close the initial approximation is to the desired frequency.

Investigate the effectiveness of using the combined use of several digital processing operations and an iterative process to determine the number of periods and frequency of a signal compared to using a fast Fourier transform to solve this problem. Conduct research for various conditions:

For different number of iterations.

With an integer and non-integer number of signal periods less than 10.

At different signal noise levels

Instructions for implementation

1. Use the program given in the appendix as a basis.

2. Make calculations based on the results of the 1st and 2nd iterations separately. As an approximate value of the number of periods at the 1st iteration, take the result obtained using the FFT; for the 2nd iteration, take the value of the number of periods refined as a result of the 1st iteration, the value of the neighborhood near the refined value of the number of periods in which reference signals are created ( parameter procch in the program) decrease by m times, m<=10.

3. When studying the effect of increasing the accuracy of determining the number of periods and frequency of a signal compared to FFT due to additional digital processing operations, measure the number of periods and frequency of the signal in the range of the number of periods from K to K+1 with a step of 0.05 and calculate the maximum, average and root-mean-square errors in determining the number of periods and frequency of the signal in this interval.

4. When studying the effect of increasing the accuracy of determining the number of periods and frequency of a signal in comparison with FFT due to additional digital processing operations when the signal is noisy, measure the number of periods and frequency of the signal in the range of the number of periods from K to K+1 with a step of 0.05 at noise amplitude values from 0 to 1 with a step of 0.1 and calculate the maximum, average and root-mean-square errors in determining the number of periods and frequency of the signal in this interval for each noise level.

5. When studying the dependence of the error in determining the number of periods and frequency of the signal on the number of periods of the signal using the FFT and the combined method, perform steps 2 and 3 at values K = 2, 3, ... 10, and plot the dependence of the maximum, average and root-mean-square errors depending on the number of signal periods.

Assignment for work.

Text of the developed program.

Graphs of the dependence of the maximum, average and root-mean-square errors at different noise levels on the number of signal periods for the FFT and the combined processing method under study.

Application. Basic program text.

%Combined use of key DSP operations

%To improve the accuracy of determining the number of periods and frequency

%"short" signal combination is used

%FFT, cross-correlation, spline approximation, resampling

kt=1024; % number of counts

f=12; %signal frequency

dt=2;% time step during measurement

kp=4.5;%number of signal periods

%1. MODEL SIGNAL GENERATION

for i=1:kt %zeroing the signal array

for i=1:kt %model signal generation

if(i>0)&(i<=kt/2)

y(i)=sin(2*3.14*kp*i/kt)*i*exp(i/k2)/1200;

if(i>kt/2)&(i<(kt))

y(i)=sin(2*3.14*kp*i/kt)*(kt-i)*exp((kt-i)/k2)/1200;

y(i)=y(i)+shum*(2*rand(1)-1);

i=1:kt; %display of the model signal in the time domain

title("Time domain")

xlabel("Sample number")

%2. FUNCTIONAL TRANSFORMATION (FFT)

bpfy=fft(y,kt);%FFT

bpf=bpfy.*conj(bpfy)/kt;%FFT

%f=1000*(0:256)/512;

plot(i(1:257),bpf(1:257));

title("Frequency domain")

xlabel("frequency")

% finding max. meaning FFT functions for array Y

for i=1:kt %search for the number of periods corresponding to the maximum FFT

%3. BENCHMARKING AND CROSS-CORRELATION

procch=0.1;%search area in percentage relative. kp_bpf

shagkor=fr*procch/3;%search step

for iii=fr-fr*procch:shagkor:fr+fr*procch %cycle to create 6 standards in the vicinity of an approximate one

%value of the number of periods determined using FFT.

%Calculation of reference signal arrays

if(i>0)&(i<=kt/2)

x(i)=sin(2*3.14*iii*i/kt)*i*exp(i/k2)/1200;

if(i>kt/2)&(i<(kt))

x(i)=sin(2*3.14*iii*i/kt)*(kt-i)*exp((kt-i)/k2)/1200;

%calculation of average values of model and reference signals

%calculation of standard deviation and coefficient. correlation of model and reference signals

x_sko=x_sko+(x(i)-x_sr)*(x(i)-x_sr);

y_sko=y_sko+(y(i)-y_sr)*(y(i)-y_sr);

kor(k)=kor(k)+(x(i)-x_sr)*(y(i)-y_sr);

kor(k)=kor(k)/(sqrt(x_sko*y_sko));

end % end of the cycle of creating standards and calculating the array of coefficients. corr.

%SPLINE APPROXIMATION AND RESAMPLING

r1=sin(xx); %only for spline fitting testing

yint=interp1(xx,kor,xi,"spline");% spline approximation correlation coefficient

%%apr=csaps(xx,r1);

apr=spaps(xkor,kor,0.000001);%%%%%%%%%%%%%%%%

plot(xkor,r1,"ro")%%%%%%%%%%%%%%%%%%

%FINDING THE REFINITED VALUE OF THE NUMBER OF SIGNAL PERIODS

cmax=max(yint); %finding the maximum coefficient. corr.

for i=1:round((k-1)/0.1+1)

if (yint(i)==cmax)

kp_int=fr-fr*procch+(i-1)*shagkor/10% specified frequency value according to the MAX coefficient function. corr.

A.B.Sergienko. Signal Processing Toolbox - review

Signal processing has always been one of the most important application areas of the MATLAB system. This is primarily evidenced by the fact that Signal Processing Toolbox was one of the first specialized packages - it appeared already in 1988, just four years after the creation of the MATLAB system itself.

To date, the Signal Processing package contains almost two hundred carefully developed specialized functions that allow solving a wide variety of signal analysis and processing problems.

The currently distributed version of MATLAB 6.1 (Release 12.1) contains the Signal Processing package version 5.1. The upcoming release of MATLAB 6.5 (Release 13) will include the Signal Processing package version 6.0.

According to the official documentation of the package, its functions are divided into twenty categories. The list below combines these categories into larger groups. So, according to their purpose, the functions of the Signal Processing package can be divided as follows:

In addition, the package includes three graphical environments:

MATLAB and its extension packages are focused primarily on digital signal processing, so functions related to analog circuit analysis are considered auxiliary. They are primarily intended to be called from other functions that use analog prototypes in the synthesis of digital filters. However, these functions themselves can be very useful. In turn, they can be divided into several groups.

The first group is functions for calculating prototype analog filters, that is, low-pass filters with a cutoff frequency of 1 rad/s. The functions return the zeros, poles, and gains of filters and allow you to calculate prototype Butterworth filters ( buttap), Chebyshev of the first kind ( cheb1ap), Chebyshev of the second kind ( cheb2ap), elliptic (Kauer) ( ellipap) and Bessel ( besselap).

The second group is analog filter conversion functions that allow you to convert a low-pass filter prototype into a low-pass filter with a different cutoff frequency ( lp2lp), into a high pass filter (HPF) with a given cutoff frequency ( lp2hp), into a bandpass filter with a given average frequency and bandwidth ( lp2bp) and into a notch filter with a given average frequency and stopband width ( lp2bs). Functions can accept and return filter descriptions as vectors of coefficients of the numerator and denominator polynomials of the transfer function, or as state space parameters.

The third group is functions for calculating analog filters with given parameters. During the calculation process, they call the functions of the first two groups. The set of initial data required for calculation includes the order of the filter, its type (low-pass filter, high-pass filter, bandpass or notch), frequency or several cutoff frequencies, and also (depending on the prototype) ripple parameters of the amplitude-frequency response (AFC). There are functions for calculating Butterworth filters ( butter), Chebyshev of the first kind ( cheby1), Chebyshev of the second kind ( cheby2), elliptic (Kauer) ( ellip) and Bessel ( beless). All these functions except the function beless, can also be used to calculate discrete filters (see below). A sign of the analog calculation option is the indication of the string "s" as the last input parameter.

The fourth group is functions for determining the required filter order based on specified frequency response parameters (cutoff frequencies of the pass and stop bands, as well as permissible ripples in the pass and stop bands). Each type of filter has its own function for determining the required order: for the Butterworth filter - buttord, for a Chebyshev filter of the first kind - cheb1ord, for a Chebyshev filter of the second kind - cheb2ord, for an elliptic filter - ellipord. Just like the functions of the previous group, these functions allow you to determine the required order for discrete filters (see below). A sign of the analog calculation option is the indication of the string "s" as the last input parameter.

The fifth group is functions for transforming the forms of description of analog linear systems. For analog systems, four such description methods are supported:

The Signal Processing package contains functions that implement mutual transformations of these four forms of representation of analog systems (only the function residue, which works with poles and residues, does not belong to the Signal Processing package, but to the underlying MATLAB library). These features are listed in the following table.

Final form	Transfer function polynomial coefficients	Zeros and poles	Poles and residues	State space
Original form
		tf2zp	residue	tf2ss
Zeros and poles	zp2tf			zp2ss
Poles and residues	residue
State space	ss2tf	ss2zp

residue can convert in both directions. The direction of the transformation is determined by the number of input and output parameters.

The functions listed in the table, except the function residue, can also carry out conversions of discrete systems, since the conversion formulas for analog and discrete systems are the same.

Finally, the sixth group should include the only function freqs, which allows you to calculate or display graphically the amplitude and phase frequency characteristics (AFC and PFC) of an analog linear system. The initial data are the coefficients of the numerator and denominator polynomials of the system transfer function.

As an example of using functions for working with analog systems, let's calculate a fourth-order elliptical low-pass filter with a cutoff frequency of 3 kHz, frequency response ripple in the passband equal to 1 dB, and signal suppression in the stopband equal to 20 dB, and then plot its frequency response and phase response .

Ellip(4, 1, 20, 2*pi*3000, "s"); % filter calculation
f = 0:10:10000; % frequency vector for calculating frequency response and phase response
h = freqs(b, a, 2*pi*f); % complex transmission coefficient
subplot(2, 1, 1)

grid
subplot(2, 1, 2)
plot(f, unwrap(angle(h))*180/pi) % plot of phase response (in degrees)
grid

The function used in the example code unwrap eliminates insignificant jumps in the phase response at 360°.

This category of functions is quite numerous and combines a variety of tools for analyzing discrete linear systems - usually presented in the form of vectors of coefficients of polynomials of the numerator and denominator of the transfer function (in z-regions).

Function freqz is a discrete analogue of the function freqs, it allows you to calculate the complex transmission coefficient or plot the frequency response and phase response of a discrete system. The initial data are the coefficients of the numerator and denominator polynomials of the system transfer function.

Similar to the function freqz function works grpdelay, which allows you to calculate or graphically display the frequency dependence of the group delay introduced by a discrete filter.

Function impz designed to calculate or graphically display the impulse response of a discrete system. The initial data, as for previous functions, are the coefficients of the polynomials of the numerator and denominator of the system transfer function.

Function zplane allows you to display the zeros and poles of the system on the complex plane, additionally depicting a unit circle that limits the permissible area for the location of the poles sustainable discrete system. The initial data can be either the vectors of coefficients of the numerator and denominator polynomials, or directly the vectors of zeros and poles of the system transfer function.

Function filternorm allows you to calculate norm discrete filter. This parameter is used, for example, when choosing scaling factors for individual sections of a filter implemented in sequential (cascade) form in order to reduce rounding errors. Two options are supported: 2-norm and -norm. The 2-norm represents root mean square the frequency response value of the filter (or, which is the same thing, the root of the sum of squared samples of the filter impulse response), and the -norm is maximum frequency response value.

Function fvtool is essentially a graphical environment designed for analyzing and visualizing the characteristics of discrete systems (Filter Visualization Tool). However, unlike other graphical environments in the package, fvtool really is function, since when called it requires the presence of input parameters - the coefficients of the polynomials of the numerator and denominator of the transfer function of the analyzed filter. A significant advantage of this function is the ability to simultaneously view characteristics several filters. The graphical user interface provided by this function is almost the same as in the FDATool filter analysis and synthesis environment. Below is an example of a function call fvtool and the graphic window it creates in the mode of showing the frequency dependence of the group delay introduced by the filter.

b = ;
a = ;
fvtool(b, a)

A large group consists of functions for transforming the forms of descriptions of discrete systems. More forms of representation are supported for discrete systems than for analogue ones:

Coefficients of the numerator and denominator polynomials of the system transfer function.
Zeros, poles and system gain (factorization of the transfer function).
Poles and residues (representation of the system transfer function as a sum of simple fractions).
State space parameters.
Representation in the form of a set of sequentially (cascaded) second-order sections.
Representation in the form of a lattice or lattice-ladder structure.

The functions that convert between different representations are listed in the following table.

Final form	Transfer function polynomial coefficients	Zeros and poles	Poles and residues	State space	Second order sections	Lattice structure
Original form
Transfer function polynomial coefficients		tf2zp	residuez	tf2ss	tf2sos	tf2latc
Zeros and poles	zp2tf			zp2ss	zp2sos
Poles and residues	residuez
State space	ss2tf	ss2zp			ss2sos
Second order sections	sos2tf	sos2zp		sos2ss
Lattice structure	latc2tf

As can be seen from the table, the same function residuez(this is an analogue of the function residue, designed to work with descriptions of discrete systems) can carry out transformations in both directions. The direction of the transformation is determined by the number of input and output parameters.

Two more functions manipulate the coefficient vectors of polynomials, modifying the position of the roots of the polynomial on the complex plane. Function polyscale multiplies all roots of the processed polynomial by a given coefficient, and the function polystab makes all the roots of the polynomial lie inside the unit circle - for this, roots exceeding unity in absolute value are inverted, that is, their modules are replaced by their inverse values, and the sign of their phases is reversed.

The remaining functions in this category are auxiliary. Function freqspace calculates a vector of evenly spaced frequency values, the function freqzplot is intended for plotting graphs of frequency characteristics, and the function unwrap allows you to eliminate insignificant jumps in phase-frequency characteristics by 2p, looking for these jumps between vector elements and shifting the desired vector fragments by 2p n.

The operation of linear discrete filtering is generally described as follows:

Here x(k)- input signal samples, y(k)- output signal samples, a i And b j- constant coefficients. Maximum of the numbers m And n is called filter order.

Previous output samples may not be used in calculations, then all a i = 0 and the filter is called non-recursive or transversal. If previous output samples are used, the filter is called recursive.

Linear discrete filtering refers to technologies for processing arbitrary data, so the corresponding function is filter- does not belong to the Signal Processing package, but is built into the MATLAB core. Function filter2, also part of the MATLAB core library, implements two-dimensional discrete filtering.

Since a constant coefficient filter is a linear time-invariant discrete system, its response to an arbitrary input signal can be represented as discrete convolution input signal with impulse response filter:

Here h(k)- samples of the filter impulse response. The impulse response is the response of the filter to a single sample of a unit value applied to the input.

Of course, calculations using the convolution formula can be implemented in practice only with a finite length of the filter impulse response. This operation is carried out by the function conv; like the implementation of the discrete filtering algorithm, it does not belong to the Signal Processing package, but to the underlying MATLAB library. In reality, implementing the conv function comes down to calling the function filter with the corresponding input parameters. Function conv2 implements two-dimensional discrete convolution. Another function of the MATLAB core library is deconv- implements convolution reversal, allowing the second input vector to be determined from the convolution result and one of the input vectors.

The basic discrete filtering functions mentioned above, as already stated, belong to the base MATLAB library; In fact, the Signal Processing package contains functions that solve more specific filtering problems.

First of all, it should be noted that the function filter makes it possible to access the initial and final internal states of the filter, thereby allowing the organization of block signal processing. Sometimes it becomes necessary to form a vector of the internal state of the filter, knowing a certain number of previous input and output samples. This calculation is performed using the function filtic.

Function fftfilt implements discrete filtering using fast Fourier transform (FFT) in combination with signal division into blocks. Only non-recursive filters can be implemented in this way. The result of the function coincides (up to computational errors) with the results of conventional filtering implemented using the function filter. However, the computational speed of FFT filtering can be significantly higher, especially if the length of the input signal is many times greater than the length of the filter's impulse response (or vice versa).

Function filtfilt allows you to compensate for the phase shift introduced by conventional filtering (in other words, this function implements filtering without introducing a time delay). This is accomplished by bidirectional signal processing. The first pass of filtering is carried out in the usual way, and then the resulting output signal is filtered a second time - from end to beginning. Due to this, phase shifts are compensated, and the resulting filter order is doubled. It should be noted that the resulting filter (equivalent to two filtering passes) does not satisfy the causality condition.

In the practical implementation of high-order recursive filters, they are often represented as second-order sections connected in series. This allows us to reduce computational errors arising from rounding errors and quantization of filter coefficients. Error analysis tools for this type of error are concentrated in the Filter Design package, and the Signal Processing package has a function sosfilt, which allows you to implement discrete data filtering using a filter presented in the form of second-order sections.

Another possible discrete filter structure is a lattice structure. To perform filtering using a filter presented in this form, use the function latcfilt.

Function medfilt1, which implements one-dimensional median filtering, belongs to nonlinear filtering algorithms. The essence of its operation is that a sliding window of a given length is applied to the input signal, the samples within the window are ordered, and the value from the middle of the ordered window (or the half-sum of the two elements closest to the middle if the window has an even length) is returned as an output sample. Median filtering is used, for example, to eliminate impulse noise (clicks) when processing audio signals. Function medfilt2, which implements a two-dimensional version of median filtering, is located in the Image Processing package.

Function sgolayfilt performs discrete filtering using a Savitzky-Golay filter. Its essence lies in the fact that the input signal is divided into blocks of a given size and, within each block, a polynomial approximation of the signal by a polynomial of a given degree is performed according to the criterion of minimum mean square error. If the degree of the polynomials is one less than the size of the blocks, the output signal will be equal to the input signal; with a smaller degree of polynomials, the signal will be smoothed. Savitzky-Golay filters are used to “clean” signals from noise.

Synthesis of a discrete filter means the selection of such sets of coefficients ( a i) And ( b i), at which the characteristics of the resulting filter satisfy the specified requirements. Strictly speaking, the design task also includes choosing a suitable filter structure, taking into account the final accuracy of the calculations. This is especially true when implementing filters “in hardware” - using specialized LSIs or digital signal processors. Effects related to the final precision of calculations can be analyzed using the functions of the Filter Design package; filter synthesis functions do not take these effects into account.

The Signal Processing package has a large number of functions that implement various algorithms for the synthesis of discrete filters. We present the main characteristics of these functions in table form, and then give some additional comments.

Function	Filter type	frequency response	Synthesis method
butter	Recursive	Butterworth	Bilinear z-transform
cheby1	Recursive	Chebyshev of the first kind	Bilinear z-transform
cheby2	Recursive	Chebyshev of the second kind	Bilinear z-transform
ellip	Recursive	Cauera (elliptical)	Bilinear z-transform
bilinear	Recursive		Bilinear z-transform
impinvar	Recursive	Arbitrary analog prototype	Invariant impulse response transformation
yulewalk	Recursive	Piecewise linear	Autoregressive method
invfreqz	Recursive	free	Minimizing the difference between the numerator of the transfer function and the product of its denominator and the desired frequency response
prony	Recursive	Synthesis from a given impulse response	Exponential Prony approximation
fir1	Non-recursive	Multiband
fir2	Non-recursive	Piecewise linear	Inverse Fourier transform using windows
firls	Non-recursive		Minimizing mean square error
fircls	Non-recursive	Piecewise constant
fircls1	Non-recursive	LPF, HPF	Minimizing the root mean square error with a maximum deviation limitation
firrcos	Non-recursive	LPF	Cosine smoothing
intfilt	Non-recursive	LPF	Minimax approximation
remez	Non-recursive	Piecewise linear with transition stripes	Minimax approximation
cremez	Non-recursive (including with nonlinear phase response and complex coefficients)	Piecewise linear with transition stripes	Minimax approximation

Methods for synthesizing discrete filters can be divided into two large groups: with and without the use of an analog prototype. When using a prototype analog filter, it is necessary to somehow represent the analog transfer function defined in the s-domain to a discrete transfer function defined in the z-domain. The Signal Processing package implements two methods of such transformation: the invariant impulse response method and the bilinear z-transform method. Both methods result in recursive discrete filters.

When using the invariant impulse response method, the impulse response of an analog prototype is sampled. The frequency response of the resulting discrete filter is accordingly a periodically repeated frequency response of the analog prototype. For this reason, this method is unsuitable for the synthesis of high-pass filters and, in general, filters whose transmission coefficient does not tend to zero with increasing frequency. The invariant impulse response method is implemented in the Signal Processing package using the function impinvar.

When using the bilinear z-transform method, the characteristics of the analog prototype are distorted only along the frequency axis. In this case, the frequency range of the analog filter (from zero to infinity) is converted to the operating frequency range of the discrete filter (from zero to half the sampling frequency). The frequency axis transformation is described by an arctangent function, so frequencies significantly lower than the sampling rate are converted approximately linearly. This method is implemented using the function bilinear for an arbitrary analog prototype. In addition, there are ready-made functions for calculating low- and high-pass filters, bandpass and notch filters using the bilinear z-transform method using analogue prototypes with frequency response of Butterworth, Chebyshev of the first and second kind, as well as Cauer (elliptic filters). This is according to the function butter, cheby1, cheby2 And ellip. All these functions can also be used to calculate analog filters (see earlier). A sign of a discrete calculation option is the absence of the string “s” in the list of input parameters. There are also functions for determining the required order of these filters based on the specified frequency response parameters (cutoff frequencies of the pass and stop bands, as well as permissible ripple in these bands). This is according to the function buttord, cheb1ord, cheb2ord, ellipord. Just like the filter synthesis functions, these functions allow you to determine the required order for analog filters (see earlier). A sign of a discrete calculation option is the absence of the string “s” in the list of input parameters.

As an example, we synthesize a fourth-order elliptical low-pass filter with the same parameters as the analog filter in one of the previous examples (cutoff frequency 3 kHz, frequency response ripple in the passband 1 dB, and signal suppression in the stopband 20 dB). Let's take the sampling frequency to be 12 kHz. After synthesis, we will construct graphs of the frequency response and phase response of the resulting filter using the function freqz.

Fs = 12000; % sampling frequency
F0 = 3000; % cutoff frequency
= ellip(4, 1, 20, F0/Fs*2); % filter calculation
freqz(b, a, , Fs); % graph output

Synthesis methods that do not use an analog prototype are called direct. They, in turn, can also be divided into two groups: methods for synthesizing recursive and non-recursive filters.

The functions of direct synthesis of non-recursive filters include the following:

Functions that implement filter synthesis by inverse Fourier transform of the desired frequency response and subsequent multiplication of the resulting impulse response by a certain weighting function (window) to attenuate the frequency response ripples that appear due to the Gibbs effect. These are the functions fir1 And fir2. This also includes the function of low-pass filter synthesis with cosine smoothing of the frequency response - firrcos. In addition, the function kaiserord allows, based on the given parameters of the frequency response, to estimate the required filter order during synthesis using the Kaiser window.
Functions that implement minimization of the standard deviation of the frequency response of the resulting filter from the specified one. These are the functions firls, fircls And fircls1. The last two functions solve an optimization problem with a limitation of the maximum deviation of the frequency response from the specified one. This allows you to avoid the appearance of large frequency response emissions near transition strips.
Functions that implement minimax optimization, that is, minimizing the peak deviation of the frequency response of the resulting filter from the specified one. The result is filters with uniform frequency response pulsations. This group includes functions remez(a standard version of the Remez method, implemented in the very first versions of the Signal Processing package) and cremez(an extended version that supports the synthesis of filters with nonlinear phase response and complex coefficients). In addition, the function remezord allows, based on the given parameters of the frequency response, to estimate the required filter order during synthesis using the Remez method.

As an example, we synthesize a non-recursive low-pass filter of the 32nd order using the Remez method with the same cutoff and sampling frequencies as in the previous example (cutoff frequency 3 kHz, sampling frequency 12 kHz). Let's set the start of the stopband to 3.5 kHz. After synthesis, we will construct graphs of the impulse response, as well as the frequency response of the resulting filter (the frequency response of the filter is linear, so it makes no sense to display its graph). We display the frequency response on a linear vertical scale to clearly demonstrate the uniformity of its pulsations.

Fs = 12000; % sampling frequency
F0 = 3000; % cutoff frequency
F1 = 3500; % start of stopband
b = remez(32, , ); % filter calculation
impz(b) % impulse response plot
= freqz(b, 1, , Fs); % complex transmission coefficient
figure
plot(f, abs(h)) % frequency response plot
grid

The functions of direct synthesis of recursive filters include the following:

yulewalk- synthesis of a recursive filter with an arbitrary piecewise linear frequency response using the Yule-Walker method.
invfreqz- this function is intended to solve the problem of identifying systems; it allows you to determine the coefficients of the numerator and denominator of the transfer function of a discrete system from a set of values of this transfer function at different frequencies.

In conclusion, we will list a number of functions that are not included in the groups listed above. Function maxflat is intended for the synthesis of a generalized Butterworth filter (for such filters the number of zeros of the transfer function exceeds the number of its poles). Function intfilt performs the synthesis of filters designed to filter the signal when performing interpolation and decimation. The operation of vector convolution can be represented as a vector-matrix product, and the matrix involved in this product can be calculated using the function convmtx. Finally, the function sgolay performs the synthesis of a Savitzky-Golay smoothing filter. Since, as described above, the Savitzky-Golay filter processes individual blocks of the signal, such a filter is not a stationary system. Therefore the function sgolay returns an entire matrix of time-varying coefficients of an equivalent non-recursive filter.

In addition to the Signal Processing package, a number of discrete filter synthesis functions are available in the Communications and Filter Design packages.

The words "spectral analysis" in the minds of many MATLAB users are strongly associated with the function fft(see further section "Functions for transforming discrete signals"), which performs a discrete Fourier transform (DFT). However, this is just a one-to-one linear transformation, giving performance deterministic signal in the frequency domain. If the analyzed signal is random, for him it only makes sense grade spectral density power, for the calculation of which it is necessary to perform averaging of the available data in one way or another. In addition, in some cases we know some additional information about the signal being analyzed, and it is advisable to take this information into account during spectral analysis.

Methods for spectral analysis of random signals are divided into two large classes - nonparametric and parametric. IN nonparametric(nonparametric) methods use only the information contained in the samples of the analyzed signal. Parametric(parametric) methods assume the presence of some statistical models random signal, and the process of spectral analysis in this case includes determining parameters this model. The term Model-Based Spectrum Analysis (MBSA) is also used.

The Signal Processing package contains functions that implement a variety of spectral analysis methods - both parametric and non-parametric (it must be emphasized once again that by spectral analysis we mean assessment power spectral density of a random process). In addition, there are functions for obtaining other averaged characteristics of random discrete signals.

To determine the spectral characteristics of a discrete random process, the average power spectrum of its fragment limited in length is calculated, and then the length of the fragment tends to infinity:

. (1)

Here x(k) - samples of a random process, T- sampling period. The bar above denotes averaging over the ensemble of realizations.

In addition, this spectrum can be expressed through the correlation function of the random process:

. (2)

This expression is a discrete analogue of the Wiener-Khinchin theorem: the spectrum of a discrete random process is the Fourier transform of its correlation function.

As already mentioned, when using nonparametric methods for calculating the spectrum of a random process, only the information contained in the signal samples is used, without any additional assumptions. The Signal Processing package implements three such methods - periodogram, Welch method and Thomson method.

A periodogram is an estimate of the power spectral density obtained from N counts one implementation random process according to definition (1) (naturally, not by taking a limit, but by averaging a finite number of terms). If a weighting function (window) is used in spectrum calculations, the resulting power spectrum estimate is called modified periodogram(modified periodogram):

Relation (2) is satisfied only for an infinite number of samples used, therefore, for any finite N The periodogram estimate of the power spectral density turns out to be displaced- it turns out that inside the sum (2) the correlation function of the signal is multiplied by a triangular weight function. In addition, it can be shown that the periodogram is not a valid estimate of the power spectral density, since dispersion such an estimate is comparable to the square of its mathematical expectation for any N. As the number of samples used increases, the periodogram values begin to fluctuate faster and faster - its graph becomes more and more jagged.

In the Signal Processing package, the periodogram (including the modified one) is calculated using the function periodogram.

To reduce the jaggedness of the periodogram, it is necessary to apply some kind of averaging. Bartlett proposed dividing the analyzed signal into non-overlapping segments, calculating a periodogram for each segment, and then averaging these periodograms. If the correlation function of the signal during the segment duration decays to negligible values, then the periodograms of individual segments can be considered independent. Welch made two improvements to Bartlett's method: the use of a weighting function and partitioning the signal into overlapping fragments. The use of a weighting function makes it possible to reduce spectrum spreading and reduce the bias in the resulting estimate of the power density spectrum at the cost of a slight deterioration in resolution. Overlapping segments was introduced in order to increase their number and reduce the variance of the estimate.

Calculations using the Welch method (also called the averaged modified periodogram method) are organized as follows: the vector of signal samples is divided into overlapping segments, each segment is multiplied by the weighting function used, modified periodograms are calculated for weighted segments, the periodograms of all segments are averaged .

The Welch method is the most popular periodogram method of spectral analysis. In the Signal Processing package it is implemented using the function pwelch.

Thomson's method implemented by the function pmtm, based on usage prolate spheroidal functions(prolate spheroidal functions). These finite duration functions provide the maximum concentration of energy in a given frequency band. In addition to the spectral assessment itself, the function pmtm can return its confidence interval. It takes some time to calculate prolate spheroidal functions, so when using the function multiple times pmtm You can speed up calculations by calculating the functions necessary for analysis in advance and storing them in the database. To work with such a database (it is a MAT file named dpss.mat) is a family of functions whose names begin with letters dpss (dpss- calculation of functions, dpssload- loading a family of functions from the database, dpsssave- saving a family of functions in the database, dpssdir- output of the database directory, dpssclear- deleting a family of functions from the database).

As an example, we will form an implementation of an exponentially correlated random process and perform its spectral analysis using the three listed methods. The random signal we need is generated by passing normal discrete white noise through a first-order recursive filter:

X0 = randn(1, 1000);
a = 0.9;
X = filter(1, , X0);

We build a periodogram:

periodogram(X, , , 1)

As you can see, the periodogram turns out to be quite jagged. Now let’s estimate the spectrum of the same implementation using Welch’s method:

pwelch(X, , , , 1)

The irregularity of the graph turns out to be much less. Finally, we use Thomson's method:

pmtm(X, , , 1)

On the output of the function pmtm The graph shows the boundaries of the confidence interval along with the power spectrum estimate.

The use of parametric methods implies the presence of some mathematical models analyzed random process. In this case, spectral analysis comes down to solving an optimization problem, that is, searching for such parameters models in which it is closest to the actually observed signal. The Signal Processing package implements a number of varieties of autoregressive analysis and two methods based on the analysis of eigenvalues and vectors of the signal correlation matrix: MUSIC (MUltiple SIgnal Classification) and EV (EigenVectors).

According to autoregressive model the signal is generated by passing discrete white noise through a "purely recursive" filter N-th order. The power spectral density of such a signal is proportional to the square of the modulus of the coefficient of the transfer function of the shaping filter. Thus, this method of spectral analysis comes down to determining the filter coefficients of a given order, estimating the power of exciting white noise, and analytically calculating the power spectral density. To determine the model coefficients, the error is minimized linear prediction signal. Theoretical analysis shows that the optimal coefficients of the model are determined only by the correlation function of the signal.

In practice, we do not know the true correlation function of the signal under study, therefore, to minimize the prediction error, we use assessments CFs obtained by time averaging. A number of autoregressive analysis methods have been developed, differing mainly in the approach to processing edge effects (that is, in the method of involving in the calculations those edge signal samples for which there is no shifted pair when calculating the CF). The Signal Processing package implements the Burg method, the covariance method, the modified covariance method, and the autoregressive Yule-Walker method.

Autoregressive spectrum analysis methods are most suitable for signals that are truly autoregressive processes. In general, these methods give good results when the spectrum of the analyzed signal has clearly defined peaks. In particular, such signals include the sum of several sinusoids with noise.

When using autoregressive methods, it is important to choose the order of the autoregressive model correctly - it should be twice the number of sinusoidal oscillations that are supposed to be contained in the analyzed signal.

Each method of autoregressive analysis in the Signal Processing package corresponds to two functions - the function of calculating model coefficients and the function of spectral analysis itself. The spectrum analysis function calls the model coefficient calculation function and then calculates the spectrum. The function names are summarized in the following table.

Method name	Function for calculating model coefficients	Spectral analysis function
Covariance method	arcov	pcov
Modified covariance method	armcov	pmcov
Berg method	arburg	pburg
Autoregressive Yule-Walker method	aryule	pylear

The exponentially correlated random signal generated in the above example is a first-order autoregressive process, therefore the listed methods of spectral analysis are quite adequate for it. Let's apply Berg's method, setting the order of the autoregressive model equal to one (this is the second parameter of the function pburg):

pburg(X, 1, , 1)

The resulting smooth curve practically coincides with the theoretical spectrum of this random process.

The MUSIC (MUltiple SIgnal Classification) method is intended for spectral analysis of signals that are the sum of several sinusoids (more precisely, in the general case, several complex exponentials) with white noise. The purpose of spectral analysis of such signals, as a rule, is not to calculate the spectrum as such, but to determine the frequencies and levels (amplitudes or powers) of the harmonic components. The MUSIC method is intended precisely for this, therefore the dependence of the signal level on frequency obtained with its help is called pseudospectrum(pseudospectrum).

The method is based on the analysis of eigenvalues and eigenvectors of the signal correlation matrix. When performing an analysis, it is necessary to specify the order of the model, that is, the number of complex exponentials expected to be contained in the signal.

In the Signal Processing package, the MUSIC method is implemented using the function pmusic, and the function rootmusic allows you to obtain estimates of the frequencies and powers of the harmonic components of the signal.

A close relative of MUSIC is the eigenvectors (EV) method. Its only difference is that in the calculation formulas the eigenvectors are multiplied by weighting coefficients inversely proportional to the corresponding eigenvalues. The literature reports that the EV method produces fewer false spectral peaks than MUSIC and generally provides a better representation of the noise spectral shape.

In the Signal Processing package, the EV method is implemented using the function peig, and the function rooteig allows you to obtain estimates of the frequencies and powers of the harmonic components of the signal.

It should be emphasized that pseudospectra are not estimates of the true power density spectrum, but are only spectral pseudo-estimates, allowing one to estimate the frequencies of sinusoidal or narrowband signal components with a resolution slightly higher than that of autoregressive methods.

The discrete Fourier transform, used in all nonparametric spectral estimation methods, implies a periodic continuation of the analyzed signal fragment. In this case, jumps may occur at the junctions of fragments, leading to the appearance of side lobes of a significant level in the spectral region. To weaken this effect, the signal is multiplied by the signal decreasing from the center to the edges before performing the DFT weight function (window). As a result, the magnitude of jumps at the junctions of segments decreases, and the level of unwanted side lobes of the spectrum also becomes smaller - the price for this is a certain broadening of the spectral peaks.

In addition to spectral analysis, weighting functions are used in the synthesis of non-recursive filters by inverse Fourier transform of the desired frequency response. In this case, they allow you to increase the signal suppression in the filter stopband due to some expansion of the passband.

Currently, the Signal Processing package contains about a dozen weighting functions. The spread of some of them is due to computational simplicity, while others are in some sense optimal.

The simplest is a rectangular window implemented by the function rectwin(in package versions up to 5.0 inclusive, this function had the name boxcar). The rectangular window corresponds to the absence of weighing; this function is included in the package only for the formal completeness of the set of weighting functions. The triangular window is implemented by the function triang, the Bartlett window also has a triangular shape (function bartlett), it only differs slightly in the method of calculation.

Several weighting functions are combinations of harmonic components. We list them in increasing order of the number of cosine terms:

Hanna window (function hann), sometimes incorrectly called the Hanning window, is a single cosine term.
Hamming window (function hamming) - one cosine term.
Blackman window (function blackman) - two cosine terms.
Blackman-Harris window (function blackmanharris) - three cosine terms.
Nuttall window (an alternative version of the Blackman-Harris window, function nuttallwin) - three cosine terms.

The remaining windows are described by more complex mathematical relationships. Gaussian window shape (function gausswin) is self-explanatory. Modified Bartlett-Hanna window (function barthannwin) is a linear combination of the Bartlett and Hanna windows. Beauman window (function bohmanwin) is the convolution of two identical cosine pulses. Chebyshev window (function chebwin) has side lobes of a fixed level (specified during calculation) and is calculated by inverse Fourier transform of the frequency response of the window. Kaiser window (function kaiser) also has a parameter that regulates the level of the side lobes and the width of the main lobe; modified Bessel functions are used when calculating this window. Tukey window (function tukeywin) is a rectangle with cosine smoothed edges. At extreme values of the smoothing coefficient, it turns into a rectangular window or Hann window.

Finally, the function window provides a general interface for calling window-specific calculation functions.

Functions belonging to this category calculate various statistical parameters of signals. Functions can be divided into several groups.

The first group relates to the calculation of correlation and covariance functions (it should be recalled here that in domestic and foreign terminology these concepts do not coincide; this review uses the foreign version adopted in MATLAB). Function xcorr allows you to estimate the correlation function of a signal or the cross-correlation function of two signals. A variant of this function designed to work with two-dimensional signals is called xcorr2. Function xcov is intended for estimating the covariance function of a signal or the mutual covariance function of two signals. Functions cov And corrcoef, included in the base MATLAB library, allow you to obtain, respectively, a covariance matrix and a matrix of correlation coefficients by averaging several realizations of random data. Function corrmtx returns a matrix of intermediate data to estimate the correlation matrix of the signal, and can also return this matrix itself.

The next group of functions calculates statistical characteristics in the frequency domain using Welch's nonparametric method (see above). Function csd intended for assessment mutual spectral density two random processes. It can also return a confidence interval for the resulting estimate. Function cohere gives an estimate for the squared modulus mutual coherence functions two random processes. Function tfe allows you to evaluate complex transfer coefficient system for the implementation of its input and output signals.

Finally, the function psdplot Used by all spectral estimation functions to plot the power spectral density. It can also be called explicitly - for example, to display a graph on a linear scale instead of the default logarithmic one, or to show several spectra on one graph.

Parametric modeling and linear prediction functions

Under parametric modeling refers to the selection of a certain mathematical model of a random process and the subsequent selection of the parameters of this model to ensure maximum correspondence between the signal generated by the model and the available real data sample.

One widely used in practice is the autoregressive (AR) model, in which a random signal is generated by passing discrete white noise through a “purely recursive” (that is, not using delayed samples of the input signal) shaping filter. Four features of the Signal Processing package - arburg, arcov, armcov And aryule- are intended to obtain estimates of the coefficients of the shaping filter and the dispersion (power) of the white noise exciting the filter. The calculation methods used by these functions were indicated earlier, in the section “Autoregressive methods”, where autoregressive spectral analysis was discussed.

If we have an estimate complex transmission coefficient systems at different frequencies, it is possible to build a realizable model of the system, the frequency response of which will be as close as possible to the measured one. The realizability of the system here means the representability of its transfer function in the form of a fractional rational function with given orders of the numerator and denominator polynomials. Parametric modeling in this case comes down to finding the optimal coefficients of the numerator and denominator polynomials of the transfer function. This problem is solved by two functions of the Signal Processing package: function invfreqs allows you to build a model of an analog system, and the function invfreqz performs a similar operation when applied to discrete systems.

Another version of the parametric modeling problem involves constructing a model of the system based on an existing estimate of its impulse response. There are two functions in the Signal Processing package for this purpose. Function prony uses the fact that the impulse response of a recursive discrete system, in the absence of multiple poles, is a sum of discrete exponential functions (in the general case complex). The algorithm implemented by this function was originally developed in the 18th century by Baron de Prony for the purpose of fitting the parameters of an exponential analytical model to experimental data. The stability of the resulting system is not guaranteed, but the first n counts ( n- the order of the numerator of the system transfer function) of its impulse response specified during the calculation exactly coincide with the specified ones.

The second function of modeling a system by impulse response is the function stmcb- does not seek to ensure an exact match of the initial fragments of the impulse responses - instead, it minimizes standard deviation the obtained characteristic from the given one, that is, the sum of the squares of the modules of the differences in the samples of the obtained and desired impulse characteristics. The function implements the Steiglitz-McBride iterative method, which reduces to repeatedly solving a system of linear equations with respect to the coefficients of the polynomials of the transfer function of the desired system.

As an example, we obtain a model of a third-order system using the Prony and Steiglitz-McBride methods, setting a triangular impulse response as a sample:

h = ; % impulse response
= prony(h, 3, 3); % Prony method
= stmcb(h, 3, 3); % Steiglitz-McBride method
% graphs of impulse responses of the resulting systems
impz(b1, a1, 30)
title("Prony")
figure
impz(b2, a2, 30)
title("Stmcb")

A comparison of the graphs clearly demonstrates the differences between the two algorithms. When using the Prony method, the first four samples of the resulting impulse response exactly coincide with the specified values, however, further deviations from the specified values increase greatly, and after the end of the specified fragment, a “tail” with a rather high level is observed, since the function prony does not make any assumptions about the required impulse response values outside a given fragment. Function stmcb minimizes quadratic playback error endless impulse response, and at the end of an explicitly specified fragment it is considered equal to zero. As a result, an exact correspondence of the samples to the given and received impulse responses is not observed (with the exception of the first), but the error in reproducing the characteristic is “spread out” more evenly across the samples.

If a discrete random process is not white noise, its samples turn out to be correlated together. This allows, knowing the correlation function of the process, predict the value of its next count. The predicted value is calculated as a linear combination of previous process samples. This is the main idea linear prediction. Linear prediction is used for parametric spectral analysis (see earlier), system identification, analysis of speech signals and compression of information during their transmission.

Models of systems based on linear prediction can be presented in various forms and, accordingly, described using different sets of parameters. A number of functions in the Signal Processing package allow you to convert a model description from one form to another. These features are listed in the following table.

Final form	Autocorrelation sequence	Reflectance coefficients	Prediction coefficients	Arcsine parameters	Logarithmic ratios	Spectral line frequencies
Original form
Autocorrelation sequence		ac2rc, schurrc
Reflectance coefficients
Prediction coefficients
Arcsine parameters
Logarithmic ratios
Spectral line frequencies

In addition, the Signal Processing package has several other functions related to linear prediction. Thus, to calculate the coefficients of a predictive filter, it is necessary to solve a system of linear equations, the matrix of which is the correlation matrix of the input signal. This matrix has a number of properties that can reduce the number of computational operations required to solve a system of linear equations. First, the correlation matrix is self-adjoint(that is, it does not change after applying it to Hermitian conjugation- combinations of transposition with complex conjugation). For a real signal, self-adjointness simply means the symmetry of the matrix about the main diagonal. Secondly, in the case of a stationary random process (and only for such processes can a predictive filter with constant parameters be used), the correlation matrix is Toeplitz matrix- along its diagonals, parallel to the main one, there are identical numbers. Finally, the right side of the system of equations represents the first column of the correlation matrix shifted by one position. Systems of linear equations with matrices having the indicated properties are called systems of Yule-Walker equations, and to solve them a recursive Levinson-Durbin method. This iterative algorithm is implemented by the function Levinson. Function rlevinson solves the inverse problem - allows you to find the vector of samples of the signal correlation function from given linear prediction coefficients.

Function lpc implements the calculation of linear prediction coefficients using the autocorrelation method and is an analogue of the function aryule(See earlier section on parametric spectral analysis). The only difference between the two functions is the MATLAB code used to calculate the correlation matrix estimate. The results they give coincide with accuracy up to computational errors.

Signal generation functions

The Signal Processing package includes a number of functions designed to generate standard waveforms commonly encountered in various signal processing applications.

Generation of non-periodic signals

All functions for generating non-periodic signals receive as parameters a vector of instants of time and additional arguments describing the parameters of the generated pulse. The returned result is a vector of samples of the resulting signal. There are functions for generating signals of the following form:

rectpuls

tripuls

gauspuls

gmonopuls

Generation of periodic signals

Functions belonging to this group receive as parameters a vector of moments in time and additional arguments describing the parameters of the generated pulse. The period of the generated signals is 2p. To generate signals with a different period, it is necessary to scale the time argument passed to the function accordingly. The returned result is a vector of samples of the resulting signal. There are functions for generating periodic signals of the following form:

square

sawtooth

diric

Generation of oscillations with varying frequency

This group includes two functions - chirp And vco. Function chirp generates oscillations, the instantaneous frequency of which varies according to one of three possible laws - linear, quadratic or exponential. The function has more capabilities vco(Voltage Controlled Oscillator - voltage-controlled oscillator), which allows you to generate oscillations with an arbitrary law of change in instantaneous frequency. In essence, this function performs frequency modulation.

Pulse train generation

Function pulstran serves to generate a finite sequence of pulses of the same shape with arbitrarily specified delays and amplitude multipliers. The shape of the pulses can be specified in one of two ways: by the name of the function generating the pulse, or by an already calculated vector of samples.

As an example, consider the use of functions pulstran And sinc to restore an analog signal from its discrete samples according to Kotelnikov’s theorem.

t = -5:0.1:10; % time for analog signal
k = 0:5; % number of samples of a discrete signal
sd = ; % discrete signal
sa = pulstran(t, , "sinc"); % reconstructed analog signal
stem(k, sd) % graph of a discrete signal
hold on
plot(t, sa, "r") % plot of the analog signal
hold off

Discrete signal conversion functions

Perhaps the most famous of the transformations of discrete signals is the discrete Fourier transform (DFT). The corresponding function using the fast Fourier transform (FFT) algorithm in MATLAB is classified as a data processing function and is built-in (functions fft And ifft- one-dimensional option, fft2 And ifft2- two-dimensional version, fftshift And ifftshift- rearrangement of halves of the vector of spectral samples to transfer the zero frequency to the middle of the vector). Actually, the Signal Processing package contains functions that implement more specific transformations.

Like any linear transformation, the DFT can be represented as multiplying the transformation matrix by a column of samples of the converted signal. This transformation matrix for the DFT is calculated by the function dftmtx.

Again, due to the linearity of the DFT, any spectral sample can be represented as the result of processing the original signal by some filter. Representing this filter in recursive form gives Goertzel's algorithm, implemented by the function goertzel. If only a few spectral samples need to be calculated, this algorithm is faster than FFT.

A close relative of the DFT is the discrete cosine transform. When calculating it, instead of the periodic continuation of the signal, which is assumed in DFT, adjacent fragments of the continued signal are mirrored in time. As a result, the signal becomes an even function of time, and its spectrum becomes real. For this reason, in the DFT formula, instead of complex exponentials, only cosines appear, which gave the name to this transformation. Forward and inverse discrete cosine transforms are calculated by the functions dct And idct respectively.

An important method for analyzing discrete number sequences is the z-transform, the result of which is a function of a complex variable z:

In a number of problems it is necessary to calculate z-transformation for points located on the spiral contour: . A computationally efficient calculation of the z-transform along such a contour uses the fast Fourier transform; it is implemented in a function czt.

When implementing some variants of the fast Fourier transform algorithm, to increase efficiency, it is necessary to rearrange the elements of the processed vector into reverse bit order(this means that in the binary representations of the vector's element numbers, the bits are read in reverse order, and then the vector is ordered according to the new element numbers). To perform such a rearrangement, use the function bitrevorder.

When analyzing narrowband signals, it can be useful to think of the signal as an oscillation with amplitude and initial phase varying over time. To obtain such a representation, a complex analytical signal, the real part of which coincides with the original signal, and the imaginary part is determined Hilbert transform from the original signal. In the frequency domain, an analytical signal is characterized by a one-sided spectrum: its spectral function is nonzero only for positive frequencies. Function hilbert calculates the analytical signal in the frequency domain by calculating the forward DFT, zeroing out half the spectrum, and then calculating the inverse DFT.

As an example, let's calculate the analytical signal for a rectangular radio pulse:

s = zeros(256,1);
s(65:192) = cos(pi/2*(0:127)"); % radio pulse counts
sa = hilbert(s); % analytical signal
f = (-128:127)/128; % values of normalized frequencies for graphs
subplot(2,1,1)
plot(f, abs(fftshift(fft(s)))) % real signal spectrum
subplot(2,1,2)
plot(f, abs(fftshift(fft(sa)))) % spectrum of the analytical signal

The lower graph clearly shows the zeroing of the spectrum in the region of negative frequencies and its doubling in the region of positive frequencies.

Cepstral analysis

Cepstral analysis is associated with homomorphic signal processing. Such processing is subject to the generalized principle of superposition: if the input signal of the system is a combination of several signals obtained using a mathematical operation A, then at the output the results of processing individual signals are combined using the operation B. Cepstral analysis has found its application, in particular, in speech processing problems. Signal Processing has several functions related to cepstral analysis.

Function rceps calculates the real cepstrum of the signal, while ignoring the information contained in the phase spectrum. The same function allows you to get minimal-phase reconstruction signal. Zeros z-transformations of the sequence of samples of such a signal lie on the complex plane inside the unit circle, and its cepstrum coincides with the cepstrum of the original signal.

Complex cepstrum calculated using the function cceps, takes into account both amplitude and phase information, so its relationship with the original signal is one-to-one. The inverse transformation, that is, the calculation of the signal using a known complex cepstrum, is performed by the function icceps.

Changing the sample rate

Transformations also include a group of functions that change the sampling frequency of a signal. Most often, the sampling rate needs to be changed by an integer factor. In case of increasing the sampling frequency, this operation is called interpolation, and in case of decrease - thinning(decimation). Any recalculation of the sampling rate requires the sequential execution of several actions. The Signal Processing package contains functions that implement both individual actions and their necessary sequences.

To perform interpolation, the required number of zeros is inserted between signal samples (function upsample), then the resulting signal is passed through a low-pass filter. Together these two stages are implemented by the function interp.

To perform decimation, the signal is passed through a low-pass filter, and then each N th count (function downsample). Together these two stages are implemented by the function decimate.

The combination of interpolation and decimation operations allows you to change the sampling frequency by a factor expressed by an arbitrary rational fraction. To increase the efficiency of calculations, such recalculation is performed by a specialized function resample, and it, in turn, calls the function upfirdn(insertion of zeros, filtering and thinning). The main difference between these two functions is that upfirdn allows you to set the impulse response of the filter used, and when using the function resample the filter is calculated automatically.

Finally, for completely arbitrary recalculation of reference moments, general interpolation functions contained in the basic MATLAB library can be used, such as interp1 And spline.

As an example, let's interpolate the signal used above to demonstrate the functions pulstran And sinc, increasing the sampling rate by four times. Zero samples are added to the edges of the signal, since the function interp requires that the length of the original sequence be at least nine samples.

sd = ; % original signal
t = 0:9; % discrete time for the original signal
sd4 = interp(sd, 4); % interpolation
t4 = (0:39)/4; % discrete time for interpolated signal
stem(t, sd) % graph of the original signal
hold on
plot(t4, sd4, "x") % plot of the interpolated signal
hold off

In the graph shown, the original signal is shown as “stems”, and the interpolated signal is shown as crosses.

This group contains a fairly large number of functions. Many of them are intended primarily for “internal use” - they are called by other functions in the package. However, a number of functions in this category have completely independent value.

Function buffer allows you to represent a vector of signal samples into a matrix of successive frames, and these frames can overlap.

Functions mod And demod perform modulation and demodulation respectively. The following types of modulation are supported: amplitude, amplitude with suppressed carrier, single-sideband, phase, frequency, quadrature, pulse-width, pulse-time.

Functions uencode And udecode implement, respectively, uniform quantization and signal restoration according to the numbers of quantization levels.

Function strips designed to display a signal graph in several lines. This can be useful if you need to look at a long signal in its entirety, and the resolution along the vertical axis is not of great importance.

A family of functions whose names begin with symbols dpss, is intended for calculating discrete prolate spheroidal functions and working with a database designed to store the calculated functions. Discrete prolate spheroidal functions are used in spectral analysis by the Thomson method (see above section “Functions of spectral analysis and statistical signal processing”, subsection “Nonparametric methods”).

Functions cell2sos And sos2cell allow you to store information about the filter, presented in the form of sequentially included sections of the second order, not only in the form of a matrix, but also in the form of an array of cells. These two functions perform representation conversion between these two forms.

Function specgram carries out the calculation of the signal spectrogram, that is, the spectra of successive signal fragments isolated using a sliding window. The calculation results can be returned as a matrix or displayed in color in time-frequency coordinates.

Function cplxpair identifies complex conjugate pairs in vectors of complex numbers.

Function eqtflength allows you to equalize the lengths of two vectors by padding the shorter one with trailing zeros.

Function seqperiod is designed to check vector elements for periodicity and determine their repetition period.

The filter synthesis and analysis program FDATool (Filter Design & Analysis Tool) is a shell for calling functions for the synthesis and analysis of discrete filters. If you have the Filter Design Toolbox package, this program also allows you to analyze the effects associated with the quantization of filter coefficients and make conversions of filter types (low-pass filter to high-pass filter, etc.).

The program allows you not only to calculate filters from scratch, but also to save and load work sessions, and edit previously saved filters. You can also import filter definitions from the MATLAB workspace or from XILINX CORE Generator (*.coe) files. Imported filters can only be analyzed.

During analysis, you can view the following filter characteristics (in the latest versions of the package you can view several graphs simultaneously):

Amplitude-frequency response (AFC).
Phase frequency response (PFC).
Frequency response and phase response simultaneously.
Group delay.
Phase delay.
Impulse response.
Transition characteristic.
The location of zeros and poles on the complex plane.
Filter coefficients.
Information about the filter structure.

The calculated filter can be used as follows:

Save the work session for subsequent loading into the FDATool program for further editing and analysis.
Export to MATLAB workspace.
Export to a text file.
Export to MAT file.
Export to a C header (*.h) file.
Export to SPTool (see below).

The figure below shows the FDATool program window with the results of calculating an elliptical low-pass filter using the bilinear z-transform method. Frequency response and group delay graphs are displayed.

The signal processing program SPTool (Signal Processing Tool) allows you to perform the following operations: import signals from MAT files or the MATLAB workspace; view signal graphs (including several simultaneously); apply various spectral analysis methods to signals and view the resulting graphs; calculate discrete filters using package functions (including by directly editing the location of zeros and poles); pass signals through filters and analyze the resulting output signals.

It should be noted that in terms of analysis and synthesis of filters, the capabilities of the SPTool program are narrower than those of the FDATool program (the only exception is that FDATool does not have the ability to directly edit the location of zeros and poles). However, these limitations are compensated by the ability to export the calculated filter from FDATool to SPTool.

The figures below show the appearance of the SPTool program window when viewing a signal graph, when analyzing a spectrum, and when editing the location of filter zeros and poles.

The WinTool (Window Design and Analysis Tool) program for the synthesis and analysis of weight (window) functions, which appeared in version 6.0 (R13) of the package, is a shell for calling the window calculation functions available in the package. This demonstrates the characteristics of a window (or several windows simultaneously) in the time and frequency domains.

The figure below shows the WinTool window with the results of calculating a 64th order Chebyshev window with a spectral side lobe level of –100 dB. The time domain plot of the window and its spectrum are displayed.

Digital signal processing: modeling in MATLAB (RUB 367.00).

Popular articles

Latest articles

Sections

Pages

Special projects

Contacts