Multirate DSP, part 2: Noninteger sampling factors (original) (raw)

Part 1 explains how to upsample and downsample by an integer factor.

Part 3 looks at oversampling in analog-to-digital converters. It applies these principles to a sigma-delta ADC and revisits the CD player case study.

Part 4 shows how to undersample bandpass signals.

12.1.3 Changing Sampling Rate by a Non-Integer Factor L/M
With an understanding of the downsampling and upsampling processes, we now study the sampling rate conversion by a non-integer factor of L /M . This can be viewed as two sampling conversion processes. In step 1, we perform the upsampling process by a factor of integer L following application of an interpolation filter H 1 (z ); in step 2, we continue filtering the output from the interpolation filter via an anti-aliasing filter H 2 (z ), and finally operate downsampling. The entire process is illustrated in Figure 12-8.

The Opportunity for Liquid Sensing

By Infineon Technologies 06.05.2025

Nine Compelling Reasons Why Menta eFPGA Is Essential for Achieving True Crypto Agility in Your ASIC or SoC 

By Menta 06.02.2025

Ex-NASA Engineer Now Guiding Manufacturers to Greater Efficiency

By MRPeasy 06.02.2025


Figure 12-8. Block diagram for sampling rate conversion.

Since the interpolation and anti-aliasing filters are in a cascaded form and operate at the same rate, we can select one of them. We choose the one with the lower stop frequency edge and choose the most demanding requirement for passband gain and stopband attenuation for the filter design. A lot of computational saving can be achieved by using one lowpass filter. We illustrate the procedure via the following simulation. Let us generate the signal x (n ) by:

x (n ) = 5 sin ((2π × 1000_n_ )/8000) + cos ((2π × 2500_n_ )/8000),

with a sampling rate of f s = 8,000 Hz and frequencies of 1 kHz and 2.5 kHz. Now we resample x (n ) to 3,000 Hz by a non-integer factor of 0.375, that is,

L/M = 0.375 = 3/8.

Upsampling is at a factor of L = 3 and the upsampled sequence is filtered by an FIR lowpass filter designed with the filter length N = 53 and a cutoff frequency of 3,250 Hz at the sampling rate of 3 × 8000 = 24,000 Hz. The spectrum for the upsampled sequence and the spectrum after application of the interpolation filter are plotted in Figure 12-9a.


Figure 12-9A. (Top) Spectrum after upsampling and (bottom) spectrum after interpolation filtering.

The sequence from step 1 can be filtered via another FIR lowpass filter designed with the filter length N = 159 and a cutoff frequency of 1,250 Hz, followed by downsampling by a factor of M = 8. The spectrum after the antialiasing filter and the spectrum for the final output y (m ) are plotted in Figure 12-9b.


Figure 12-9B. (Top) Spectrum after anti-aliasing filtering and (bottom) spectrum after downsampling.

Note that the anti-aliasing filter removes the frequency component of 2.5 kHz to avoid aliasing. This is because after downsampling, the Nyquist limit is 1.5 kHz. As we discussed previously, we can select one filter for implementation. We choose the FIR lowpass with N = 159 and a cutoff frequency of 1,250 Hz because its bandwidth is smaller than that of the interpolation filter. The MATLAB implementation is listed in Program 12-3.

Therefore, three steps are required to accomplish the process:

  1. Upsampling by a factor of L = 3
  2. Filtering the upsampled sequence by an FIR lowpass filter designed with the filter length N = 159 and a cutoff frequency of 1,250 Hz at the sampling rate of 3 × 8000 = 24,000 Hz
  3. Downsampling by a factor of M = 8.

(Click to enlarge)
Program 12-3. MATLAB program for changing sampling rate with a noninteger factor.

Example 12.3.
Given a sampling conversion DSP system with the following specifications, determine the filter length and cutoff frequency for the combined antialiasing filter H (z ), and window types, respectively, if the window design method is used:


Figure 12-10A. Sampling conversion in Example 12.3.

Solution: The filter frequency specifications and corresponding block diagram are developed in Figure 12-10b.


Figure 12-10B. Filter frequency specifications for Example 12.3.

Specifications for the interpolation filter H 1 (z ):

Specifications for the anti-aliasing filter H 2 (z ):

Combined specifications H (z ):

We use the FIR filter with the Hamming window. Since

Δ_f_ = (f stop – f pass )/f sL = (3000 – 2500)/18000 = 0.0278,

the length of the filter and the cutoff frequency can be determined by

N = 3.3/Δ_f_ = 3.3/0.0278 = 118.8.

We choose N = 119, and

f c = (f pass + f stop )/2 = (3000 + 2500)/2 = 2750 Hz.

12.1.4 Application: CD Audio Player
In this application example, we will discuss principles of the upsampling and interpolation-filter processes used in the CD audio system to help the reconstruction filter design.

Each raw digital sample recorded on the CD audio system contains 16 bits and is sampled at the rate of 44.1 kHz. Figure 12-11 describes a portion of one channel of the CD player in terms of a simplified block diagram.


Figure 12-11. Sample rate conversion in the CD audio player system.

Let us consider the situation without upsampling and application of a digital interpolation filter. We know that the audio signal has a bandwidth of 22.05 kHz, that is, the Nyquist frequency, and digital-to-analog conversion (DAC) produces the sample-and-hold signals that contain the desired audio band and images thereof. To achieve the audio band signal, we need to apply a reconstruction filter (also called a smooth filter or anti-image filter) to remove all image frequencies beyond the Nyquist frequency of 22.05 kHz. Due to the requirement of the sharp transition band, a higher-order analog filter design becomes a requirement.

The design of the higher-order analog filter is complex and expensive to implement. In order to relieve such design constraints, as shown in Figure 12-11, we can add the upsampling process before DAC, followed by application of the digital interpolation filter (assume L = 4). Certainly, the interpolation filter design must satisfy the specifications studied in the previous section on increasing the sampling rate by an integer factor. Again, after digital interpolation, the audio band is kept the same, while the sampling frequency is increased fourfold (L = 4), that is, 44.1 × 4 = 176.4 kHz.

Since the audio band of 22.05 kHz is now relatively low compared with the new folding frequency (176.4/2 = 88.2 kHz), the use of a simple first-order or second-order analog anti-image filter may be sufficient. Let us look at the following simulation.

A test audio signal with a frequency of 16 kHz and a sampling rate of 44.1 kHz is generated using the formula

x (n ) = sin ((2π × 16000_n_ )/441000).

If we use an upsampling factor of 4, then the bandwidth would increase to 88.2 kHz. Based on the audio frequency of 16 kHz, the original Nyquist frequency of 22.05 kHz, and the new sampling rate of 176.4 kHz, we can determine the filter length as

Δ_f_ = (22.05 – 16)/176.4 = 0.0343.

Using the Hamming window for FIR filter design leads to

N = 3.3/Δ_f_ = 96.2.

We choose N = 97. The cutoff frequency therefore is

f c = (16 + 22.05)/2 = 19.025 kHz.

The spectrum of the interpolated audio test signal is shown in Figure 12-12, where the top plot illustrates that after upsampling, the audio test signal has the frequency of 16 kHz, along with image frequencies coming from 44.1 kHz – 16 = 28.1 kHz, 44.1 + 16 = 60.1 kHz, 88.2 – 16 = 72.2 kHz, and so on. The bottom graph describes the spectrum after the interpolation filter. From lowpass FIR filtering, the interpolated audio signal with a frequency of 16 kHz is observed.


Figure 12-12. (Top) The spectrum after upsampling and (bottom) the spectrum after applying the interpolation filter.

Let us examine the corresponding process in time domain, as shown in Figure 12-13. The upper left plot shows the original samples. The upper right plot describes the upsampled signals. The lower left plot shows the signals after the upsampling process and digital interpolation filtering. Finally, the lower right plot shows the sample-and-hold signals after DAC. Clearly, we can easily design a reconstruction filter to smooth the sample-and- hold signals and obtain the original audio test signal. The advantage of reducing hardware is illustrated. The MATLAB implementation can be seen in Program 12-4.


Figure 12-13. Plots of signals at each stage according to the block diagram in Figure 12-11.

(Click to enlarge)
Program 12-4. MATLAB program for CD player example.

12.1.5 Multistage Decimation
The multistage approach for downsampling rate conversion can be used to dramatically reduce the anti-aliasing filter length. Figure 12-14 describes a two-stage decimator.


Figure 12-14. Multistage decimation.

As shown in Figure 12-14, a total decimation factor is M = M 1 ×M 2 . Here, even though we develop a procedure for a two-stage case, a similar principle can be applied to general multistage cases. Using the two-stage decimation in Figure 12-15, the final Nyquist limit is f s /2_M_ after final downsampling. So our useful information bandwidth should stop at the frequency edge of f s /2_M_ .


Figure 12-15. Stopband frequency edge for the anti-aliasing filter at stage 1 for two-stage decimation.

Next, we need to determine the stop frequency edge for the anti-aliasing lowpass filter at stage 1 before the first decimation process begins. This stop frequency edge is actually the lower frequency edge of the first image replica centered at the sampling frequency of f s /M 1 after the stage 1 decimation. This lower frequency edge of the first image replica is then determined by

f s /M 1 – f s /(2_M_ ),

After downsampling, we expect that some frequency components from f s /M 1 – f s /2_M_ to f s /M1 to be folded over to the frequency band between f s /2_M_ and f s /2_M_ 1 . However, these aliased frequency components do not affect the final useful band between 0 Hz to f s /2_M_ and will be removed by the anti-aliasing filter(s) in the future stage(s). As illustrated in Figure 12-15, any frequency components beyond the edge f s /M 1 – f s /2_M_ can fold over into the final useful information band to create aliasing distortion. Therefore, we can use this frequency as the lower stop frequency edge of the anti-aliasing filter to prevent the aliasing distortion at the final stage. The upper stopband edge (Nyquist limit) for the anti-image filter at stage 1 is clearly f s /2, since the filter operates at f s samples per second. So the stopband frequency range at stage 1 is from f s /M 1 – f s /2_M_ to f s /2. The aliasing distortion, introduced into the frequency band from f s /2_M_ to f s /2_M_ 1 , will be filtered out after future decimation stage(s).

Similarly, for stage 2, the lower frequency edge of the first image developed after stage 2 downsampling is

f s /(M 1 M2 )– f s /(2_M_ ) = f s /(2_M_ ).

As is evident in our two-stage scheme, the stopband frequency range for the second anti-aliasing filter at stage 2 should be from f s /(2_M_ ) to f s /(2_M_ 1 ).

We summarize specifications for the two-stage decimation as follows:

Filter requirements for stage 1:

Filter requirements for stage 2:

Example 12.4 illustrates the two-stage decimator design.

Example 12.4.
Determine the anti-aliasing FIR filter lengths and cutoff frequencies for the two-stage decimator with the following specifications and block diagram:


Figure 12-16A. Multistage decimation in Example 12.4.

Solution:

M = 240 kHz/8 kHz = 30 = 10 × 3

We choose M 1 = 10 and M 2 = 3; there could be other choices. Figure 12-16b shows the block diagram and filter frequency specifications.


Figure 12-16B. Filter frequency specifications for Example 12.4.

Note that the lower stopband edge can be determined as

f stop = f s /M 1 – f s /(2 × M ) = 240000/10 – 240000/(2 × 30) = 20000 Hz

Δ_f_ = (f stop – f pass )/f s = (20000 – 3400)/240000 = 0.06917.

The length of the filter and the cutoff frequency can be determined by

N = 3.3/Δ_f_ = 47.7.

We choose N = 49, and

f c = (f pass + f stop )/2 = (20000 + 3400)/2 = 11700 Hz.

Filter specification for H 2 (z ):

Note that

Δ_f_ = (f stop – f pass )/f sM1 = (4000 – 3400)/24000 = 0.025.

The length of the filter and the cutoff frequency can be determined by

N = 3.3/Δ_f_ = 132.

We choose N = 133, and

f c = (f pass + f stop )/2 = (4000 + 3400)/2 = 3700 Hz.

The reader can verify for the case by using only one stage with a decimation factor of M = 30. Using the Hamming window for the FIR filter, the resulting number of taps is 1,321, and the cutoff frequency is 3,700 Hz. Thus, such a filter requires a huge number of computations and causes a large delay during implementation compared with the two-stage case.

The multistage scheme is very helpful for sampling rate conversion between audio systems. For example, to convert the CD audio at the sampling rate of 44.1 kHz to the MP3 or digital audio type system (professional audio rate), in which the sampling rate of 48 kHz is used, the conversion factor L /M = 48/44.1 = 160/147 is required. Using the single-stage scheme may cause impractical FIR filter sizes for interpolation and downsampling. However, since L /M = 160/147 = (4/3)(8/7)(5/7), we may design an efficient three-stage system, in which stages 1, 2, and 3 use the conversion factors L /M = 8/7,L /M = 5/7, and L /M = 4/3, respectively.

12.2 Polyphase Filter Structure and Implementation
Due to the nature of the decimation and interpolation processes, polyphase filter structures can be developed to efficiently implement the decimation and interpolation filters (using fewer numbers of multiplications and additions). As we will explain, these filters are all-pass filters with different phase shifts (Proakis and Manolakis, 1996), thus we call them polyphase filters.

Here, we skip their derivations and illustrate implementations of decimation and interpolation using simple examples. Consider the interpolation process shown in Figure 12-17, where L = 2.


Figure 12-17. Upsampling by a factor of 2 and a four-tap interpolation filter.

We assume that the FIR interpolation filter has four taps, shown as:

H (z ) = h (0) + h (1)z -1 + h (2)z -2 + h (3)z -3

and the filter output is

y (m ) = h (0)w (m ) + h (1)w (m – 1) + h (2)w (m – 2) + h (3)w (m – 3).

For the purpose of comparison, the direct interpolation process shown in Figure 12-17 is summarized in Table 12-1, where w (m ) is the upsampled signal and y (m ) the interpolated output. Processing each input sample x (n ) requires applying the difference equation twice to obtain y (0) and y (1). Hence, for this example, we need eight multiplications and six additions.


Table 12-1. Results of the direct interpolation process in Figure 12-17. (8 multiplications and 6 additions for processing each input sample x(n)).

The output results in Table 12-1 can be easily obtained by using the polyphase filters shown in Figure 12-18.


Figure 12-18. Polyphase filter implementation for the interpolation in Figure 12-17. (4 multiplications and 4 additions for processing each input sample x(n)).

In general, there are L polyphase filters. Having the designed interpolation filter H (z ) of the N taps, we can determine each bank of filter coefficients as follows:

ρk (n ) = h (k + nL ) for k = 0, 1, … , L – 1 and n = 0, 1, … , (N /L) – 1. (12.12)

For our example L = 2 and N = 4, we have L – 1 = 1, and N /L – 1 = 1, respectively. Hence, there are two filter banks, ρ0 (z ) and ρ1 (z ), each having a length of 2, as illustrated in Figure 12-18. When k = 0 and n = 1, the upper limit of time index required for h (k + nL ) is k + nL = 0 + 1 – 2 = 2. When k = 1 and n = 1, the upper limit of the time index for h (k + nL ) is 3. Hence, the first filter ρ0 (z ) has the coefficients h (0) and h (2). Similarly, the second filter ρ1 (z ) has coefficients h (1) and h (3). In fact, the filter coefficients of ρ0 (z) are a decimated version of h(n) starting at k = 0, while the filter coefficients of ρ1 (z) are a decimated version of h (n ) starting at k = 1, and so on.

As shown in Figure 12-18, we can reduce the computational complexity from eight multiplications and six additions down to four multiplications and four additions for processing each input sample x (n ). Generally, the computation can be reduced by a factor of L as compared with the direct process.

The commutative model for the polyphase interpolation filter is given in Figure 12-19.


Figure 12-19. Commutative model for the polyphase interpolation filter.

Example 12.5.
Verify y (1) in Table 12-1 using the polyphase filter implementation in Figures 12.18 and 12.19, respectively.

Solution:

Applying the polyphase interpolation filter as shown in Figure 12-18 leads to

w 0 (n ) = h (0)x (n ) + h (2)x (n – 1)

w 1 (n ) = h (1)x (n ) + h (3)x (n – 1);

when n = 0,

w 0 (0) = h (0)x (0)
w 1 (0) = h (1)x (0).

After interpolation, we have

y 0 (m ):w 0 (0) 0 …

and

y 1 (m ): 0 w 1 (0) 0 …

Note: there is a unit delay for the second filter bank. Hence

m = 0, y 0 (0) = h (0)x (0), y 1 (0) = 0

m = 1, y 0 (1) = 0, y 1 (1) = h (1)x (0).

Combining two channels, we finally get

m = 0, y (0) = y 0 (0) + y 1 (0) = h (0)x (0),

m = 1, y (1) = y 0 (1) + y 1 (1) = h (1)x (0).

Therefore, y (1) matches that in the direct interpolation process given in Table 12-1.

Applying the polyphase interpolation filter using the commutative model in Figure 12-19, we have

y 0 (n ) = h (0)x (n ) + h (2)x (n – 1)

y 1 (n ) = h (1)x (n ) + h (3)x (n – 1);

when n = 0,

m = 0, y (0) = y 0 (0) = h (0)x (0) + h (2)x ( – 1) = h (0)x (0)

m = 1, y (1) = y 1 (0) = h (1)x (0) + h (3)x ( – 1) = h (1)x (0).

Clearly, y (1) = h (1)x (0) matches the y (1) result in Table 12-1.

Next, we will explain the properties of polyphase filters (all-pass gain and possible different phases). Each polyphase filter ρk (n ) operating at the original sampling rate f s (assuming 8 kHz) is a downsampled version of the interpolation filter h (n ) operating at the upsampling rate Lf s (32 kHz, assuming an interpolation factor of L = 4). Considering that the designed interpolation FIR filter coefficients h (n ) are the impulse response sequence having a flat frequency spectrum up to a bandwidth of f s = 2 (assume a bandwidth of 4 kHz with a perfect flat frequency magnitude response, theoretically) at the sampling rate of Lf s (32 kHz), we then downsample h (n ) to obtain polyphase filters by a factor of L = 4 and operate them at a sampling rate of f s (8 kHz).

The Nyquist frequency after downsampling should be (Lf s = 2) = L = f s = 2 (4 kHz); at the same time, each downsampled sequence ρk (n ) operating at f s (8 kHz) has a flat spectrum up to f s = 2 (4 kHz) due to the f s = 2 (4 kHz) bandlimited sequence of h (n ) at the sampling rate of f s (32 kHz). Hence, all of the polyphase filters are all-pass filters. Since each polyphase ρk (n ) filter has different coefficients, each may have a different phase. Therefore, these polyphase filters are the all-pass filters having possible different phases, theoretically.

Next, consider the following decimation process in Figure 12-20.


Figure 12-20. Decimation by a factor of 2 and a three-tap anti-aliasing filter.

Assuming a three-tap decimation filter, we have

H (z ) = h (0) + h (1)z –1 + h (2)z –2

w (n ) = h (0)x (n ) + h (1)x (n – 1) + h (2)x (n – 2).

The direct decimation process is shown in Table 12-2 for the purpose of comparison. Obtaining each output y (m ) requires processing filter difference equations twice, resulting in six multiplications and four additions for this particular example.

(Click to enlarge)
Table 12-2 Results of direct decimation process in Figure 12-20. (6 multiplications and 4 additions for obtaining each.

The efficient way to implement a polyphase filter is given in Figure 12-21.


Figure 12-21. Polyphase filter implementation for the decimation in Figure 12-20. (3 multiplications and 1 addition for obtaining each output y (m)).

Similarly, there are M polyphase filters. With the designed decimation filter H (z ) of the N taps, we can obtain filter bank coefficients by

ρk (n ) = h (k + nM )
for k = 0, 1, … , M – 1 and n = 0, 1, … , (N/M ) – 1. (12.13)

For our example, we see that M – 1 = 1 and (N /M ) – 1 = 1 (rounded up). Thus we have two filter banks. Since k = 0 and n = 1, k + nM = 0 + 1 – 2 = 2. The time index upper limit required for h (k + nM ) is 2 for the first filter bank ρ0 (z ). Hence ρ0 (z ) has filter coefficients h (0) and h (2).

However, when k = 1 and n = 1, k + nM = 1 + 1 – 2 = 3, the time index upper limit required for h (k + nM ) is 3 for the second filter bank, and the corresponding filter coefficients are required to be h (1) and h (3). Since our direct decimation filter h (n ) does not contain the coefficient h (3), we set h (3) = 0 to get the second filter bank with one tap only, as shown in Figure 12-21. Also as shown in that figure, achieving each y (m ) needs three multiplications and one addition. In general, the number of multiplications can be reduced by a factor of M .

The commutative model for the polyphase decimator is shown in Figure 12-22.


Figure 12-22. Commutative model for the polyphase decimation filter.

Example 12.6.
Verify y (1) in Table 12-2 using the polyphase decimation filter implementation in Figure 12-21.

Solution:

Using Figure 12-21, we write the difference equations as

y 0 (m ) = h (0)w 0 (m ) + h (2)w 0 (m – 1)

y 1 (m ) = h (1)w 1 (m ).

Assuming n = 0, n = 1, n = 2, and n = 3, we have the inputs x (0), x (1), x (2), and x (3), and

w 0 (m ): x (0) x (2) …

Delaying x (n ) by one sample and decimating it by a factor of 2 leads to

w 1 (m ): 0 x (1) x (3) … .

Hence, applying the filter banks yields the following:
m = 0, we have inputs for each filter as

w 0 (0) = x (0) and w 1 (0) = 0

then

y 0 (0) = h (0)w 0 (0) + h (2)w 0 ( – 1) = h (0)x (0)

y 1 (0) = h (1)w 1 (0) = h (1) – 0 = 0.

Combining two channels, we obtain

y (1) = y 0 (1) + y 1 (1) = h (0)x (0) + 0 = h (0)x (0)

m = 1, we get inputs for each filter as

w 0 (1) = x (2) and w 1 (1) = x (1),

then

y 0 (1) = h (0)w 0 (1) + h (2)w 0 (0) = h (0)x (2) + h (2)x (0)

y 1 (1) = h (1)w 1 (1) = h (1)x (1).

Combining two channels leads to

y (1) = y 0 (1) + y 1 (1) = h (0)x (2) + h (2)x (0) + h (1)x (1).

We note that y (1) is the same as that shown in Table 12-2. Similar analysis can be done for the commutative model shown in Figure 12-22.

Note that wavelet transform and subband coding are also in the area of multirate signal processing. We do not pursue these subjects in this book. The reader can find useful fundamental information in Akansu and Haddad (1992), Stearns (2003), Van der Vegte (2002), and Vetterli and Kovacevic (1995).

Related articles

Printed with permission form Academic Press, a division of Elsevier. Copyright 2007. “Digital Signal Processing, Fundamentals and Applications” by Li Tan. For more information about this title and other similar books, please visit www.elsevierdirect.com.

ADVANCED TECHNOLOGY, AUDIO, AUTOMOTIVE, CONSUMER ELECTRONICS & APPLIANCES, CPLD, DSP, FPGA, GAMING, HARDWARE DEVELOPMENT, LEGISLATION, MEDICAL DEVICES & SYSTEMS, MICROCONTROLLER, MICROPROCESSOR, MOBILE, PLD, POLITICS, SEMICONDUCTORS, SOFTWARE, TRANSPORTATION