A Simplified Matlab Function for Power Spectral Density

In an earlier post, I showed how to compute power spectral density (PSD) of a discrete-time signal using the Matlab function pwelch.  Pwelch is a useful function because it gives the correct output, and it has the option to average multiple Discrete Fourier Transforms (DFTs).  However, a typical function call has five arguments, and it can be hard to remember how to set them all and how they default.

In this post, I create a simplified PSD function by putting a wrapper on pwelch that sets some parameters and converts the output units from W/Hz to dBW/bin.  The function is named psd_simple.m, and its code is listed in the appendix.

Fractional Delay FIR Filters

Consider the following Finite Impulse Response (FIR) coefficients:

b = [b₀ b₁ b₂ b₁ b₀]

These coefficients form a 5-tap symmetrical FIR filter having constant group delay [1,2] over 0 to f_s/2 of:

D = (ntaps - 1)/2 = 2      samples

For a symmetrical filter with an odd number of taps, the group delay is always an integer number of samples, while for one with an even number of taps, the group delay is always an integer + 0.5 samples.  Can we design a filter with arbitrary delay, say 9.3 samples?  The answer is yes -- It is possible to design a non-symmetrical FIR filter with arbitrary group delay which is approximately constant over a wide band, with approximately flat magnitude response [3,4].  Let the desired group delay be:

D = (ntaps - 1)/2 + u 

   = D₀ + u     samples,            (1)

where we call u the fractional delay and -0.5 <= u <= 0.5.  D₀ is the fixed portion of the total delay; it is determined by ntaps.  The appendix lists a simple Matlab function frac_delay_fir.m to compute FIR coefficients for a given value of u and ntaps.  The function provides coefficients with approximately flat delay and frequency responses over a frequency range approaching 0 to f_s/2.

In this post, we'll present a couple of examples using the function, then discuss the theory behind it.  Finally, we'll look at an example of a fractional delay lowpass FIR filter with arbitrary cut-off frequency.

Digital Filtering in the Frequency Domain

Time domain digital filtering, whether implemented using FIR or IIR techniques, has been very well documented in literature and been thoroughly used in many base band processor designs. However, with the advent of software defined radios as well as CPU support in more recent baseband processors, it has become possible and often desirable to filter signals in software rather than digital hardware. Whereas, time domain digital filtering can certainly be implemented in software as well, it becomes highly inefficient as the number of filter taps grows. Frequency domain filtering, using FFT and IFFT operations, is significantly more efficient and surprisingly easy to understand. This document introduces the reader to frequency domain filtering both in theory and in practice via a MatLab script.

The DFT of Finite-Length Time-Reversed Sequences

Recently I've been reading papers on underwater acoustic communications systems and this caused me to investigate the frequency-domain effects of time-reversal of time-domain sequences. I created this article because there is so little coverage of this topic in the literature of DSP.

Model Signal Impairments at Complex Baseband

In this article, we develop complex-baseband models for several signal impairments: interfering carrier, multipath, phase noise, and Gaussian noise. To provide concrete examples, we'll apply the impairments to a QAM system. The impairment models are Matlab functions that each use at most seven lines of code. Although our example system is QAM, the models can be used for any complex-baseband signal.

Update To: A Wide-Notch Comb Filter

This article presents alternatives to the wide-notch comb filter described in Reference [1].

A Wide-Notch Comb Filter

This article describes a linear-phase comb filter having wider stopband notches than a traditional comb filter.

The Risk In Using Frequency Domain Curves To Evaluate Digital Integrator Performance

This article shows the danger in evaluating the performance of a digital integration network based solely on its frequency response curve. If you plan on implementing a digital integrator in your signal processing work I recommend you continue reading this article.

Reduced-Delay IIR Filters

This document describes a straightforward method to significantly reduce the number of necessary multiplies per input sample of traditional IIR lowpass and highpass digital filters.

Physical Principles Involved in Transistor Action

A Quadrature Signals Tutorial: Complex, But Not Complicated

Quadrature signals are based on the notion of complex numbers and perhaps no other topic causes more heartache for newcomers to DSP than these numbers and their strange terminology of j operator, complex, imaginary, real, and orthogonal. If you're a little unsure of the physical meaning of complex numbers and the j = √-1 operator, don't feel bad because you're in good company. Why even Karl Gauss, one the world's greatest mathematicians, called the j operator the "shadow of shadows". Here we'll shine some light on that shadow so you'll never have to call the Quadrature Signal Psychic Hotline for help. Quadrature signal processing is used in many fields of science and engineering, and quadrature signals are necessary to describe the processing and implementation that takes place in modern digital communications systems. In this tutorial we'll review the fundamentals of complex numbers and get comfortable with how they're used to represent quadrature signals. Next we examine the notion of negative frequency as it relates to quadrature signal algebraic notation, and learn to speak the language of quadrature processing. In addition, we'll use three-dimensional time and frequency-domain plots to give some physical meaning to quadrature signals. This tutorial concludes with a brief look at how a quadrature signal can be generated by means of quadrature-sampling.

Introduction of C Programming for DSP Applications

Appendix C of the book : Real-Time Digital Signal Processing: Implementations, Application and Experiments with the TMS320C55X

Digital Envelope Detection: The Good, the Bad, and the Ugly

Recently I've been thinking about the process of envelope detection. Tutorial information on this topic is readily available but that information is spread out over a number of DSP textbooks and many Internet web sites. The purpose of this blog is to summarize various digital envelope detection methods in one place. Here I focus of envelope detection as it is applied to an amplitude-fluctuating sinusoidal signal where the positive-amplitude fluctuations (the sinusoid's envelope) contain some sort of information. Let's begin by looking at the simplest envelope detection method.

Use Matlab Function pwelch to Find Power Spectral Density - or Do It Yourself

In this article, I'll present some examples to show how to use pwelch. You can also "do it yourself", i.e. compute spectra using the Matlab fft or other fft function. As examples, the appendix provides two demonstration mfiles; one computes the spectrum without DFT averaging, and the other computes the spectrum with DFT averaging.

Understanding the 'Phasing Method' of Single Sideband Demodulation

There are four ways to demodulate a transmitted single sideband (SSB) signal. Those four methods are: synchronous detection, phasing method, Weaver method, and filtering method. Here we review synchronous detection in preparation for explaining, in detail, how the phasing method works. This blog contains lots of preliminary information, so if you're already familiar with SSB signals you might want to scroll down to the 'SSB DEMODULATION BY SYNCHRONOUS DETECTION' section.

Design IIR Butterworth Filters Using 12 Lines of Code

While there are plenty of canned functions to design Butterworth IIR filters [1], it's instructive and not that complicated to design them from scratch. You can do it in 12 lines of Matlab code.

Fractional Delay FIR Filters

Consider the following Finite Impulse Response (FIR) coefficients:

b = [b₀ b₁ b₂ b₁ b₀]

These coefficients form a 5-tap symmetrical FIR filter having constant group delay [1,2] over 0 to f_s/2 of:

D = (ntaps - 1)/2 = 2      samples

For a symmetrical filter with an odd number of taps, the group delay is always an integer number of samples, while for one with an even number of taps, the group delay is always an integer + 0.5 samples.  Can we design a filter with arbitrary delay, say 9.3 samples?  The answer is yes -- It is possible to design a non-symmetrical FIR filter with arbitrary group delay which is approximately constant over a wide band, with approximately flat magnitude response [3,4].  Let the desired group delay be:

D = (ntaps - 1)/2 + u 

   = D₀ + u     samples,            (1)

where we call u the fractional delay and -0.5 <= u <= 0.5.  D₀ is the fixed portion of the total delay; it is determined by ntaps.  The appendix lists a simple Matlab function frac_delay_fir.m to compute FIR coefficients for a given value of u and ntaps.  The function provides coefficients with approximately flat delay and frequency responses over a frequency range approaching 0 to f_s/2.

In this post, we'll present a couple of examples using the function, then discuss the theory behind it.  Finally, we'll look at an example of a fractional delay lowpass FIR filter with arbitrary cut-off frequency.

Python For Audio Signal Processing

This paper discusses the use of Python for developing audio signal processing applications. Overviews of Python language, NumPy, SciPy and Matplotlib are given, which together form a powerful platform for scientific computing. We then show how SciPy was used to create two audio programming libraries, and describe ways that Python can be integrated with the SndObj library and Pure Data, two existing environments for music composition and signal processing.

Method to Calculate the Inverse of a Complex Matrix using Real Matrix Inversion

This paper describes a simple method to calculate the invers of a complex matrix. The key element of the method is to use a matrix inversion, which is available and optimised for real numbers. Some actual libraries used for digital signal processing only provide highly optimised methods to calculate the inverse of a real matrix, whereas no solution for complex matrices are available, like in [1]. The presented algorithm is very easy to implement, while still much more efficient than for example the method presented in [2]. [1] Visual DSP++ 4.0 C/C++ Compiler and Library Manual for TigerSHARC Processors; Analog Devices; 2005. [2] W. Press, S.A. Teukolsky, W.T. Vetterling, B.R. Flannery; Numerical Recipes in C++, The art of scientific computing, Second Edition; p52 : “Complex Systems of Equations”;Cambridge University Press 2002.

Design IIR Filters Using Cascaded Biquads

This article shows how to implement a Butterworth IIR lowpass filter as a cascade of second-order IIR filters, or biquads. We'll derive how to calculate the coefficients of the biquads and do some examples using a Matlab function biquad_synth provided in the Appendix. Although we'll be designing Butterworth filters, the approach applies to any all-pole lowpass filter (Chebyshev, Bessel, etc). As we'll see, the cascaded-biquad design is less sensitive to coefficient quantization than a single high-order IIR, particularly for lower cut-off frequencies.