Sampling rate conversion

Not Applicable

Not Applicable

Not Applicable

This invention relates to converting an input data stream having a first sampling rate to an output data stream having a second sampling rate.

Sampling rate conversion is performed when a discrete-time system needs to combine data streams with differing sampling rates. This is a common occurrence in audio and video applications. For example, audio sources may occur at different sampling rates, such as 44.1 kHz (kilohertz) for compact disks (CDs) and 48 kHz for digital audio tape (DAT).

One way of performing sampling rate conversion is to use a combination of up-sampling, filtering, and down-sampling. This can work if the ratio of the two sampling rates equals a rational number. In many cases, however, the two sampling rates are not related by a simple rational number. This occurs frequently when trying to synchronize data streams that are arriving from separate digital networks. For example, in audio applications, it is common to have a discrete-time system operating at 44.1 kHz while its input data stream is arriving at roughly 44.09 kHz with time variations. The receiving portion of the discrete-time system must adjust for these time variations during sampling rate conversion.

One way of performing sampling rate conversions on time-varying data streams uses finite impulse response (FIR) filters and linear weighted interpolation based on neighboring samples.

It is an important object to provide improved sampling rate conversion.

In general, in one aspect, the invention is directed to converting an input signal having a first sampling rate to an output signal having a second sampling rate. The invention features obtaining an intermediate sampling value from the input signal and filtering the intermediate sampling value to obtain the output signal. The intermediate sampling value corresponds to a sample taken at the second sampling rate on a continuous-time representation of the input signal.

This aspect of the invention may include one or more of the following features. Obtaining the intermediate sampling value includes obtaining a coefficient that corresponds to the sample from the continuous-time representation of the input signal and determining the intermediate sampling value based on the coefficient and an impulse value. The impulse value corresponds to the input signal. The coefficient is a value of the continuous-time representation of the input signal taken at the sample. The coefficient is determined from a previous coefficient by multiplying the previous coefficient by a constant. The constant corresponds to a difference between the sample and a previous sample.

Filtering the intermediate sampling value includes adding the intermediate sampling value to a previous value that corresponds to a previous sample taken at the second sampling rate. Filtering is performed by a first order digital filter.

A second intermediate sampling value is obtained from the input signal and the second intermediate sampling value is filtered to obtain the output signal. The second intermediate sampling value corresponds to a second sample taken at the second sampling rate on the continuous-time representation of the input signal. In this case, filtering is performed using a second order digital filter.

One of the first and second sampling rates may be a compact disk sampling rate and the other of the first and second sampling rates may be a digital audio tape sampling rate. The continuous-time representation of the input signal may be a decaying exponential function or a sinusoidal function.

Other features, objects and advantages of the invention will become apparent from the following description and drawings.

FIG. 1A is a system diagram of a prior art sampling rate conversion process; and FIG. 1B shows graphs of the signals included in the system diagram of FIG. 1A;

FIG. 2A is a system diagram of the prior art sampling rate conversion process; FIG. 2B shows graphs of the signals included in the system diagram of FIG. 2A; FIG. 2C shows an intermediate sampling value obtained from an input signal according to the invention; and FIG. 2D is a system diagram of a sampling rate conversion process that uses the intermediate sampling value according to the invention;

FIG. 3 is a flowchart showing the sampling rate conversion process that uses the intermediate sampling value according to the invention;

FIG. 4 shows the impulse response h(t) broken down into separate components h w (t) and h i (t);

FIG. 5 is a system diagram of a sampling rate conversion process that includes the separate components h w (t) and h i (t);

FIG. 6 a system diagram of a sampling rate conversion process that is derived from the system diagram of FIG. 5 and that includes the intermediate sampling value according to the invention;

FIG. 7 shows signals present in the system of FIG. 6 when a first order digital filter is used;

FIG. 8 is a close-up view of signal w c (t) from FIG. 7;

FIGS. 9A through 9F are graphs showing how w c (t) varies over time;

FIG. 10 shows signals present in the system of FIG. 6 when a second order digital filter is used;

FIG. 11 is a graph showing a 9 th order filter;

FIG. 12 includes graphs showing the impulse response, h(t), of the 9 th order filter and the individual partial fraction expansion components of the impulse response; and

FIG. 13 is a block diagram of a system implementing the process of FIG. 3 according to the invention.

With reference now to the drawings, FIGS. 1A and 1B show prior art sampling rate conversion process 10 . In those figures, a digital signal x(n) 11 having a first sampling rate T 0 is converted to a digital signal y(n) 12 having a second sampling rate T 1 . The conversion process 10 includes generating a time-based signal x c (t) 14 using discrete-to-continuous (D/C) converter 15 , where x c (t) is: x c  ( t ) = ∑ k = - ∞ ∞     x  [ k ]  δ  ( t - kT 0 )

δ is the impulse function, and k is a sample.

The signal x c (t) comprises impulses 14 a to 14 c which occur at the same points as the data 11 a in digital signal x(n). Filter H(s) 16 converts x c (t) to a continuous-time representation y c (t) 17 of x(n). In this example, filter 16 is a first order, low-pass filter with a cutoff frequency that is the smaller of n/T 0 and n/T 1 (the Nyquist rate). Continuous-to-discrete (C/D) converter 19 then samples y c (t) at rate T 1 to obtain digital signal y(n).

FIGS. 2A and 2B show how a single data point 20 sampled at T 0 {x(n)} is converted to digital signal y(n) 21 sampled at T 1 using the system of FIG. 1 A. The continuous-time representation y c (t) 22 for a single data point is a decaying exponential for a first order filter H(s). One or both of T 0 and T 1 may be time varying relative to the other. That is, T 0 and T 1 may be not be related by a rational number.

The sampling rate conversion process 25 according to the invention is shown in FIG. 3 . Process 25 obtains ( 301 ) an intermediate sampling value w(n) 26 (FIG. 2C) from input digital signal x(n). Value w(n) corresponds to a sample taken at T 1 (the output rate). Value w(n) is then digitally filtered ( 302 ) to obtain y(n).

A system diagram 27 for implementing process 25 is shown in FIG. 2 D. In system diagram 27 , block 29 includes transfer function H(s) for obtaining w(n) from x(n) and H d (z) 30 , where H d (z) 30 is the filter used to obtain y(n). The following describes how to obtain w(n) from x(n) and then to obtain y(n) from w(n).

Assume that H(s) in block 29 is a first order filter with a system function of H  ( s ) = 1 s + a .

H(s) thus has a pole at s=−a and a zero at infinity. The impulse response, h(t), of a system with function H(s) is

h ( t )= e −at u ( t ),

where u(t) is the unit step function defined by u  ( t ) = 1     for     t ≥ 0 = 0     otherwise .

If “a” is real and positive then, by definition, the impulse response h(t) is real and the system is stable. Also assume that T 1 =1. This does not restrict sampling rate conversion process 25 , since any change in the sampling rate can be achieved by choosing an appropriate value for T 0 .

The impulse response h(t) of system 27 (FIG. 4A) can be represented as a finite-duration or time-limited window h w (t) (FIG. 4B) convolved with a train of decaying impulses h i (t) (FIG. 4 C). More specifically, h  ( t ) =  - at  u  ( t ) = α t  u  ( t ) = h w  ( t ) * h i  ( t ) ,

where α=e −a , h w  ( t ) = α t  ( u  ( t ) - u  ( t - 1 ) ) = α t  r  ( t ) , and h i  ( t ) = ∑ k = 0 ∞     α k  δ  ( t - k )

The resulting system diagram 31 is shown in FIG. 5, in which H(s) from FIGS. 1A and 2A is replaced by h w (t) 32 and h i (t) 34 .

If the sampling performed by C/D converter 35 and the impulses in h i (t) are both spaced by one second in time, C/D converter 35 and h i (t) 34 can be interchanged. The resulting system diagram 36 is shown in FIG. 6, in which h i (t) is replaced by its equivalent discrete time filter h d (n) 37 , where

h d ( n )=α n u ( n ).

The output of C/D converter 35 is the intermediate sampling value w(n). The value of w(n) may be determined as follows.

Assuming that x(n)=δ(n−n 0 ) and x c (t)=δ(t−n 0 T 0 ), then w c (t)=h w (t−n 0 T 0 ), which starts at time n 0 T 0 and has a duration of one second. FIG. 7 shows w c (t) 39 and x c (t) 40 for a given x(n) 41 . The value for w(n) 42 is sampled from w c (t) at the second/output sampling rate. The function w c (t) is a continuous-time representation of the input signal x(n). For a first order filter, as in this example, w c (t) is a decaying exponential function that corresponds to x(n).

C/D converter 35 (FIG. 6) has a single nonzero output at a time n=┌n 0 T 0 ┐, where the notation “┌n 0 T 0 ┐” indicates the smallest integer value that is greater than or equal to n 0 T 0 . This is so because C/D converter 35 samples at integer time values only. The amplitude of the sample of w(n) will thus be

h w (┌ n 0 T 0 ┐−n 0 T 0 ),

as shown in FIG. 8 (which is a close-up view of w c (t) from FIG. 7 ). If the input x(n) is the impulse function δ(n−n 0 ), then the value of w(n) is h w multiplied by x(n), or

w ( n )= h w (┌ n 0 T 0 ┐−n 0 T 0 )δ( n−┌n 0 T 0 ┐),

where δ(n−┌n 0 T 0 ┐) is the impulse function at time ┌n 0 T 0 ┐. Defining τ(n) to be “┌n 0 T 0 ┐−n 0 T 0 ” (see FIG. 8 for a graphical representation of τ(n)) makes

w ( n )= h w (τ( n 0 ))δ( n−┌n 0 T 0 ┐).

To determine the response of the system to an input x(n), where x(n) contains multiple nonzero samples m, the equation is as follows ω  [ n ] = ∑ m = - ∞ ∞     x  [ m ]  h ω  ( r  [ m ] )  δ  [ n - ⌈ mT 0 ⌉ ]

Evaluating this sum may involve stepping through x(n) one sample at a time, multiplying the resulting value by the coefficient h w (τ(n)), and then adding the resulting product to the output. The coefficient h w (τ(n)) is obtained by sampling the continuous-time representation h w at τ(n) (see FIG. 8 ), as described below. The coefficient h w (τ(n)) corresponds to the sample x(n), since h w is the continuous-time representation of x(n).

To obtain y(n) from w(n), w(n) is filtered using the digital filter h d (n) 37 (FIG. 6 ). As noted above,

h d ( n )=α n u ( n ).

Taking the Z-transform of h d (n) yields

H d ( z )= y ( z )/ w ( z )=1/(1 −αz −1 ).

Rearranging the terms and taking the inverse Z transform of H d (z) leads to the difference equation

y ( n )= w ( n )+α y ( n −1),

which defines the filter for obtaining y(n) from w(n).

The following describes determining the values of h w (τ(n)) that are used in determining w(n).

The value of τ(n) changes as h w (t) changes for each x(n). To determine how τ(n) changes over time, Δ(n) (or simply “Δ”) is defined to be

Δ( n )=τ( n )−τ( n −1).

Noting that τ(n)=┌n 0 T 0 ┐−n 0 T 0 , results in the following Δ  [ n ] =    ⌈ nT 0 ⌉ - nT 0 - ⌈ ( n - 1 )  T 0 ⌉ + ( n - 1 )  T 0 =    ⌈ nT 0 ⌉ - ⌈ ( nT 0 - T 0 ⌉ - T 0 =    { ⌈ T 0 ⌉ - T 0 , τ  [ n - 1 ] < T 0 - ⌊ T 0 ⌋ ⌊ T 0 ⌋ - T 0 , τ  [ n - 1 ] ≥ T 0 - ⌊ T 0 ⌋

Thus, Δ(n) can take on two values: τ(n) either increases by a known quantity, corresponding to a positive value of Δ, or τ(n) decreases by a known quantity, corresponding to a negative value of Δ. So, τ(n) remains within the range of

0≦τ( n )<1.

FIGS. 9A to 9 F show shifted and scaled copies of h w (t). For sample n=1, nT 0 =0.9 (FIG. 9 B); for sample n=2, nT 0 =1.8 (FIG. 9 C); for sample n=3, nT 0 =2.7 (FIG. 9 D); and so on. The value of Δ in this example is either 0.1 or −0.9, depending on the sampling time. That is, for nT 0 =0.9, ┌nT 0 ┐ is 1 and ┌nT 0 ┐−nT 0 (i.e., τ(n)) is 0.1; for nT 0 =1.8, ┌nT 0 ┐ is 2 and ┌nT 0 ┐−nT 0 is 0.2; for nT 0 =2.7, ┌nT 0 ┐ is 3 and ┌nT 0 ┐−nT 0 is 0.3; and so on. The difference between a current value of ┌nT 0 ┐−nT 0 (τ(n)) and a previous value of ┌nT 0 ┐−nT 0 (τ(n−1))is 0.1 (Δ). Thus, in this example, the value of Δ increases by 0.1 as h w (t) changes over time. Alternatively, when nT 0 =8.1, ┌nT 0 ┐ is 9, and ┌nT 0 ┐−nT 0 is 0.9. Following this progression leads to a Δ value of −0.9.

Knowing the value of Δ, it is possible to determine τ(n) from τ(n−1) and, thus, to determine h w (τ(n)) from h w (τ(n−1)) and Δ. More specifically, from above, h w  ( τ  ( n ) ) =    α τ  ( n ) =    α ( τ  ( n - 1 ) + Δ ) =    h w  ( τ  ( n - 1 ) )  α Δ .

Thus, at each sample, the new coefficient h w (τ(n)) can be calculated from the previous coefficient h w (τ(n−1)) in real-time by multiplying h w (τ(n−1)) by the constant α Δ . The value of α Δ can be determined beforehand, stored in memory, and then retrieved for the multiplication. Using α Δ to update h w (τ(n)) eliminates the need to store a table of values for h w (τ(n)).

The foregoing describes sampling rate conversion using a first order digital filter h d (n). The invention, however, is not limited to use with first order filters. For example, a second order filter may be used. In the case of a second order filter, assume that the system function H(s) 29 in FIG. 2D has one complex conjugate pole pair, one zero on the real axis, and one zero at infinity. For such a system,

h ( t )=γ t sin( wt +φ) u ( t ).

Representing the impulse response of the second order system as the convolution of a time limited portion h w (t) and an infinite train of impulses yields h  ( t ) =    γ t  sin  ( ωt + o )  u  ( t ) =    h ω  ( t ) × h i  ( t ) .

where

h, r ( t )=γ t {sin( wt −φ)τ( t )−sin( w ( t −2)+φ)τ( t −1)}.

and h i  ( t ) = ∑ k = 0 ∞     γ k  sin  ( ω  ( k ÷ 1 ) ) sin     ω  δ  ( t - k ) .

As above, interchanging the filtering and sampling operations yields the system shown in FIG. 6, where h d  [ n ] = γ n  sin  ( ω  ( n + 1 ) ) sin     ω  u  [ n ] .

When x(n)=δ(n−n 0 ), w c (t) has a time-width of two seconds for the second order case. This is shown in FIG. 10, which shows the second order case 45 in relation to the first order case 46 from FIG. 7 . As shown, in the second order case 45 , w c (t) is a portion (time-limited window) of a sinusoidal function and the two samples 49 and 50 are taken in respective phases 51 and 52 of that sinusoidal function to obtain w(n).

The first of the samples 49 occurs at t=┌n 0 T 0 ┐ and the second of these samples 50 occurs at t=┌n 0 T 0 ┐+1. The amplitude of first sample 49 is h w (τ(n 0 )) and the amplitude of second sample 50 is h w (τ(n 0 )+1). So, for an input of x(n)=δ(n−n 0 ),

w ( n )= h w (τ( n 0 ))δ( n−┌n 0 T 0 ┐)+ h w (τ( n 0 )+1)δ( n−┌n 0 T 0 ┐−1),

This is analogous to the equation for w(n) described for first order filters above. For an arbitrary input x(m), the intermediate sampling value w(n) is as follows ω  [ n ] = ∑ m = - ∞ ∞  x  [ m ]  { h ω , τ  [ m ] }  δ  ⌈ n - [ mT 0 ⌉ ] + h ω  ( τ  [ m ] - 1 )  δ  [ n - ⌈ mT 0 ⌉ - 1 ] }

The equations for determining the first coefficient h w (τ(n)) and the second coefficient h w (τ(n)+1) are as follows:

h w (τ[ n ])= A·h w (τ[ n −1])+ B·h w (τ[ n −1]+1)

h w (τ[ n ]+1)= C·h w (τ[ n −1])+ D·h w (τ[ n −1]+1),

where A =    γ Δ  ( cos     Δω + sin     Δωcotω ) B =    γ ( Δ - 1 )  sin     Δω sin     ω C =    - γ ( Δ - 1 )  sin     Δω sin     ω D =    γ Δ  ( cos     Δω - sin     Δωcotω )

In the foregoing equations, Δ is as described above and γ and ω are known second order filter properties.

As above, the constants A, B, C and D can be determined beforehand, stored in memory, and retrieved when it comes time to perform the necessary multiplication. For each sample, new coefficients h w (τ(n)) and h w (τ(n)+1) are determined in real-time from the old coefficients h w (τ(n−1)) and h w (τ(n−1)+1) using only four multiplication operations, i.e., by A, B, C and D. Since there are two possible values for Δ (as above), two different version of A, B, C and D should be determined, stored, and applied as appropriate.

To obtain the digital filter h d (n) for the second order case, we take the Z transform of h d  [ n ] = γ n  sin  ( ω  ( n + 1 ) ) sin     ω  u  [ n ]

to obtain H d  ( z ) = 1 1 - ( 2  γcosω )  z - 1 + γ 2  z - 2

Rearranging the terms and taking the inverse Z transform of H d (Z) leads to the difference equation

y=[n]=w[n ]+(2γ cos w ) y[n −1]−γ 2 y[n −2].

This involves only two multiplication operations per sample to obtain y(n) from w(n).

The invention can be used with higher order filters as well. For example, assume that H(s) has N (N>2) poles and N−1 zeros in the finite s-plane. A partial fraction expansion can be performed on H(s) to break down H(s) into the sum of several first and/or second order pieces, each of which can be solved as described above. Each of these pieces can be implemented in parallel and the results added together to find a solution to the Nth order case. To implement the Nth order case, the impulse response of the system is defined as h  ( t ) = ∑ i = 1 N p     β i  α i t  u  ( t ) + ∑ i = 1 N p     θ i  γ i t  sin  ( ω i  t + o i )  u  ( t ) ,

where N r is the number of real poles and N p is the number of complex conjugate pole pairs. A similar process may be used for the case where H(s) has repeated poles.

To summarize, the sampling rate conversion process for a given H(s) having N r real poles and N p pairs of complex conjugate poles (for a total N poles) is as set forth in FIG. 3 . That is, each of the intermediate sampling values w(n) are obtained ( 301 ) by obtaining ( 301 a ) the coefficients, h w , that are needed from the continuous-time representation of the input signal. For N (N≧1) total poles, N such coefficients will be needed. Values for w(n) are determined ( 301 b ) by stepping through the input samples one at a time and, at each step (i) multiplying the value of the sample by the coefficients h w and adding the results to w(n), and (ii) updating the coefficients h w in preparation for the next iteration. After values for w(n) have been obtained, w(n) is filtered ( 302 ) using an appropriate filter and the results of the filtering summed to obtain output y(n).

The following describes process 25 as performed on two digital audio streams that are to be combined and output through the same DAC (digital-to-analog converter). Each input is a two-channel audio stream. However, one of the inputs is arriving at 44.1 kHz and the other is arriving at 44.09794 kHz. Process 25 is used to up-sample the second input by a fraction T 0 =0.9995328.

To do this, H(s) is chosen to be an elliptic filter with a passband edge frequency of 2π(20 kHz/44.1 kHz) radians per second and a stopband edge frequency of 2π(24.1 kHz/44.1 kHz) radians per second. The passband is allowed ±0.2 dB (decibels) of ripple and the stopband is allowed at least 80 dB of attenuation. To match these specifications, a 9 th order elliptic filter having the response shown in FIG. 11 is used. The impulse response h(t) and its individual partial fraction components h 1 (t) to h 5 (t) are shown in FIG. 12 .

For generating the values of w(n), process 25 performs one multiplication operation per sample for the real pole. Process 25 performs two multiplication operations for each of the four-pole pairs. These are performed twice because there are two channels. Thus, a total of eighteen multiplication operations are performed for every sample to generate values of w(n). To update the h w coefficients, one multiplication operation is performed for the real pole and four multiplication operations are performed for each of the pole pairs, for a total of seventeen multiplication operations. To filter the values of w(n) to obtain the output y(n), nine multiplication operations are performed for each of the two channels, adding a total of eighteen more multiplication operations. Thus, process 25 performs a total of 53 multiplication operations. In real time, these multiplication operations are performed for every sample, i.e., 44,100 times per second, for a total of 2.34 million multiplication operations per second.

In a standard FIR sampling conversion technique, a time-varying impulse response having a length of twenty-seven is required (thirty-one if it is a linear phase). So,twenty-seven multiplication operations per channel are needed. This leads to about 2.38 million multiplication operations per second, which is comparable to the number of multiplication operations required by process 25 . The difference, however, is that the FIR technique requires a large table of impulse responses/coefficients to be stored. Such a table is on the order of seven thousand points. For process 25 , only about fifty-two values (e.g., Δ, A, B, C and D above) are stored. Thus, relatively little memory is needed to perform the sampling rate conversion process of the invention.

FIG. 13 shows a device 55 (e.g., a computer) for performing sampling rate conversion process 25 . Device 55 includes a processor 56 , a memory 57 , and a,storage medium 59 , e.g., a hard disk (see view 60 ). Storage medium 59 stores computer-executable instructions 61 for performing process 25 . Processor 56 executes computer-executable instructions 61 out of memory 57 to perform process 25 .

Process 25 , however, is not limited to use with the hardware/software configuration of FIG. 13; it may find applicability in any computing or processing environment. Process 25 may be implemented in hardware, software, or a combination of the two (e.g., using an ASIC (application-specific integrated circuit) or programmable logic). Process 25 may be implemented in one or more computer programs executing on programmable computers that each includes a processor, a storage medium readable by the processor (including volatile and non-volatile memory and/or storage elements), at least one input device, and one or more output devices. Program code may be applied to data entered using an input device to perform process 25 and to generate output information. The output information may be applied to one or more output devices.

Each such program may be implemented in a high level procedural or object-oriented programming language to communicate with a computer system. However, the programs can be implemented in assembly or machine language. The language may be a compiled or an interpreted language.

Each computer program may be stored on a storage medium or device (e.g., CD-ROM, hard disk, or magnetic diskette) that is readable by a general or special purpose programmable computer for configuring and operating the computer when the storage medium or device is read by the computer to perform process 25 . Process 25 may also be implemented as a computer-readable storage medium, configured with a computer program, where, upon execution, instructions in the computer program cause the computer to operate in accordance with process 25 .

Other embodiments not described herein are also within the scope of the following claims. For example, an input signal whose sampling rate is to be converted may be partitioned into blocks. After each block, new h w coefficients are determined. This block-based approach permits changing of the sampling rate T 0 for each block. The invention can be used to mix digital signals from any asynchronous sources in real-time, not just digital audio.

Representative Matlab code for performing sampling rate conversion in accordance with process 25 using an N th (N>1) order filter is shown in the attached Appendix.