Sampling rate conversion is a common operation in digital signal processing systems allowing first to adapt the original sampling rate to target requirements, e. g. in audio systems from 48 kHz to 44.1 kHz or in software defined radios (SDR) from analog to digital converter (ADC) sampling rate to exact multiples of signal bandwidth like for DAB+ sampling rate of 2048 MHz and secondly to significantly decrease computational demands on hardware. However, correct and efficient resampling may not be as straightforward as one would think. Let’s have a look at a little bit of theory and key points during implementation step by step.
Decimation
With the decimation in signal processing we mean changing the sample rate down by an integer factor. This is mostly the first conversion block in the SDR receiver systems to efficiently decrease hardware demands and, finally, significantly lower the cost for a customer. The main idea is very simple, just to keep only every \( M^{th}\) sample from the incomming signal and others throw away. This idea can be written as
$$x_d[n]=x[nM]=x_c(nMT)$$
where \(x_d\) is \(n^{th}\) sample of original signal \(x\) after decimation by a factor \(M\) and \(x_c\) is its continuous-time representation with sampling period \(T\). This implies that the new sampling period is \(T_d=MT\) which obviously decreases the Nyquist rate and to prevent aliasing the signal needs to be properly filtered1. The image aliasing by decimation can be expressed in spectral domain using discrete-time Fourier transform (DTFT)
$$X_d\left(e^{j\omega}\right)=\frac{1}{M}\sum_{k=0}^{M-1}X\left(e^{\frac{j\omega}{M}-\frac{j2\pi k}{M}}\right)$$
where the uppercase \(X\) stands for the frequency domain representation of \(x\) and \(\omega=2\pi f\), where \(f\) is a frequency. Graphically the Nyquist zone compression with decimation is expressed as follows.
Finally, the decimation is the cascade of (anti-aliasing) filter and decimator or compressor.
So far theprocess was straightforward without any trick that we even needn’t to show the process of decimation in time domain graphically, but let’s do this now.
There is an arbitrary wave sampled with sampling period \(T\) depicted as black and after decimation with \(T_d\) as red. Thinking about implementation one will notice the filter located prior compressor is the most demanding block within a decimator and, actually, it is calculating samples that are thrown away later in the compressor. Implemented, the compressor doesn’t provide any complex task than ignoring some samples. The trick to do here is not to calculate the ignored samples already in the filter. So let’s have a look at the filter implementation closer. Consider a finite impulse response (FIR) filter with an original impulse response \(h\) that can be rewritten as successively delayed decomposed versions
$$h[n]=\sum_{k=0}^{M-1}h_k[n-k]$$
where
$$h_k[n]=\begin{cases}h[n+k] & n=\text{integer multiple of }M \\ 0 & \text{otherwise}\end{cases}$$
The filter has then decomposed structure
where \(H_k(z^M)\) is the coresponding transfer function of \(h_k(n)\). From the implementation point of view the filtering or a convolution is repeated multiplications and additions of correct samples of the signal and filter impulse response, thus it is wise to rewrite each step in a table. Let’s create \(h_k(n)\) first.
The original impulse response is decimated with factor \(M\) and upsampled with factor \(M\), i.e. interleaved with zeros to its original length. There are \(M\) different impulse responses created from successively shifted samples. Now, let’s show the convolution with the input signal in a table view.
Since the filter is used prior decimator, we keep only one output sample from the table diagram above, let’s say we keep only \(y_5\) and then \(y_6\) and \(y_7\) would be thrown away, so there is no need to calculate them. The samples \(\vec{x}\) are coming in the rate \(T_s\), but \(\vec{y}\) samples are leaving with a rate \(\frac{T_s}{M}\), thus multiplication with \(\vec{h_k}\) can be made sequentially. To summarize, the multiply-and-accumulate (MAC) operations performed at \(T_s\) is reduced by the factor \(M\), but the coefficient for a single MAC unit is not constant and it must process \(M\) coefficients periodically.
This is a significant saving that can be combined with the specific properties of the impulse response. One typically used is of even or odd impulse responses, i.e. filters with linear phase, and the implementation is straightforward. A multirate processing can be another technique, which leads to a cascade of filters with relaxed sampling rate and shorter impulse response. The length of impulse response is the property that significantly complicates the single-rate implementations of decimators. The filter length grows with the factor \(M\) while keeping the same filter specifications. The example in figure below shows the two impulse responses simply generated as an inverse DTFT of rectangular windows for different cut-off frequencies. The one of the wider window vanishes \(\frac{T_1}{T_2}\) faster. Well, if an interested reader would try to replicate these curves, he will find there is a scaling factor \(M\) applied inversely to the impulse responses causing gain of the transfer function. As we will show later the central sample (\(x_0\)) equals \(\frac{1}{M}\), therefore the impulse responses actually vanishes to the same level at the edges, but starting at different level. On the other hand, scaling the coefficients is wise to keep the processing accuracy.
Eventhough the multirate processing uses multiple filters, the overall hardware demands can be lower. Cascading multiple decimators of factor 2 is quite common due to advantageous property of the so caled half-band filters. First, let’s look at the theory to give a better understanding how it works.
Looking again at the typical impulse response of the FIR filters used for decimator one can notice it crosses zero many times. Generally the digital filter coefficients are the samples of the impulse response sometimes close to the zero crossing point. Let’s think about the sampling of the impulse response right at the time of crossing intuitively every second sample, with advice of the impulse response \(T_c=\frac{1}{2}\) above. What will generally be the resulting transfer function, or rule?
First, let’s the DTFT into more specific form for the FIR filter design, it means we assume real time domain impulse response.
$$X_k=\sum_{n=0}^{N-1}x_ne^{-j2\pi\frac{k}{N}n}=x_0+(-1)^kx_{\frac{N}{2}}+\sum_{n=1}^{\frac{N}{2}-1}x_ne^{-j2\pi\frac{k}{N}n}+x_{N-n}e^{-j2\pi\frac{k}{N}(N-n)}$$
where we isolated the real and complex expressions. Now, assuming \(x_n=x_{N-n}\), which is the condition of spectral symmetry we get
$$X_k=x_0+(-1)^kx_{\frac{N}{2}}+2\sum_{n=1}^{\frac{N}{2}-1}x_n\cos\left(2\pi\frac{k}{N}n\right)$$
which is in fact the cosine transform. Now, assume the even samples of the impulse response are zero excluding \(x_0\).
$$X_k=x_0+2\sum_{n=0}^{\frac{N}{4}-1}x_{2n+1}\cos\left(2\pi\frac{k}{N}(2n+1)\right)$$
Further examining, we can express the equation for special values of \(k\) and get the rules for the transfer function.
$$X_k=\begin{cases} x_0 & k=\frac{N}{4}\\ x_0+2\sum_{n=0}^{\frac{N}{4}-1}x_{2n+1}\cos\left(2\pi\frac{m}{N}(2n+1)\right) & k=m\\ x_0-2\sum_{n=0}^{\frac{N}{4}-1}x_{2n+1}\cos\left(2\pi\frac{m}{N}(2n+1)\right) & k=\frac{N}{2}-m \end{cases}$$
To rephrase this we see that the transfer function has to be symmetrical around sample \(k=\frac{N}{4}\) , which is the global offset value, and in our case for the normalized tranfer function the value is \(\frac{1}{2}\). Generally it is normalized absolute integral of the transfer function.
To implement the half-band filter lets apply the polyphase decomposition the same way as earlier to the general FIR filter. First, look at the table view of the impulse response.
The second vector \(\vec{h_1}\) contains only the central sample, which leads to the following signal diagram of the half-band decimating filter.
The sub-filter \(H_0(z)\) runs at half the input sample rate which again leads to resouce savings, when MAC units are shared between two coefficients, or channels. Additional savings can be reached implementing odd type of the HB-filter2.
So far it was the theoretical basis allowing you to write the script or a simple program in C to start testing and initiate to pave the way to functional efficient implementation in a real hardware. We spent much time working on optimal implementation for high speed processing on FPGA and, finally, created IP cores enabling all highlited featured. Contact us for more information.
Interpolation
Unlike decimation, the interpolation refers to increasing the sample rate by an integer factor. As a counterpart to the decimator an interpolator is commonly employed as a last block in the SDR on the transmitting side before DAC to save hardware demands and again lower the final cost of the system. The task of interpolation is to create new samples in between the original and somehow estimate its value. In radio (or audio) signal processing we estimate the samples with linear convolution, i.e. filtration, and therefore band-limiting thespectrum of the signal. Let’s relate the sampled signal \(x[n]\), continuous signal \(x_c[t]\) and interpolated sampled signal \(x_i[n]\) as
$$x_i[n]=x\left[\frac{n}{L}\right]=x_c\left[\frac{nT}{L}\right]$$
where \(T\) is a sampling period of the original signal and \(L\) is the interpolation factor. The signal diagram of the interpolator system is as follows
First, the incomming samples are expanded with zeros in upsampler or expander. The expanded signal can be expressed as
$$x_e[n]=\sum_{k=-\infty}^{\infty}x[k]\delta[n-kL]$$
which in frequency domain shows the expanded spectrum as frequency-scaled original signal spectrum.
$$X_e\left(e^{j\omega}\right)=\sum_{k=-\infty}^{\infty}x[k]e^{-j\omega Lk}=X\left(e^{j\omega L}\right)$$
In other words, the expansion increases the sampling frequency and the Nyquist zone with it, and as a result there is \(L\) original zones shrinked in a new single nyquist zone as depicted in the figure.
The blue dotted signals in an expanded spectrum should be removed by filtering after expansion to finish interpolation3. The red spectra in the Nyquist zone 1′ in the expanded spectrum are baseband images created as a result of time-discretization. For the sake of sampled signal band-limitation we cannot talk about Nyquist zones as commonly used with continuous-time signals, therefore we count them unsing apostrophe.
The interpolation system is not very complicated, but as well as classic decimator use filtration on the side of high sampling rate, i.e. after expander, and thus the filter itself is the most demanding component during implementation. Looking closer at the signal, the expander simply adds \(L-1\) zeros between samples as depicted in figure.
As the samples are 0 we don’t need to include them into convolution. Looking carefuly at the polyphase decomposed filter as derived in the decimator section we can see that each time only one filter branch results non-zero sample because the input sample pulse train is expanded with zeros.
From that we can state that there is needed only a single filter with set of decimated impulse responses, which are loaded on every outcomming sample. How to createe the set of impulse responses is crear from the table view below.
And the implementation of the filter is then depicted in table view.
The filter demands can be the significantly decreased, up to pre-expander filter level, which is great benefit at almost no design cost. Unfortunately, there is no gain employing symmetry of the FIR filters, because the ommited samples lies in the signal vector and not in the impulse response vector as it was in the decimator case. But multirate processing using half-band filters can have additonal effect.
As well as in the case of decimator we worked hard to prepare optimized module of the fast and efficient interpolator for FPGA. If you don’t want to spent months on developing your own, do not hesitate to contact us either to quote the IP itself or a training.
Fractional resampling
So far we discussed the decreasing and increasing the sampling rate by an integer value. Unfortunately this may not satisfy many specific needs in signal processing and rational factor is needed. The straightforward solution to this problem is to employ a combination of the interpolator and decimator modules right in this order. Sadly, this is very common approach in a systems around us. Well, why not to do it thsi simple way? Starting to think about how to put the things together one will reveal the main issue with this sub-optimal strategy – the need to increase the sampling rate (much) higher than is needed for the output signal. This often brings a hard limit in the resampling factors, you will hardly implement the resampling factor 7/9 with the ADC running at 180 MHz. Let’s have a look at the rational resampling a little bit closer an I promise you that you can do it at 180 MHz sample rate maximum.
In the signal diagram we joined the interpolator and decimator
and the filters merged. Its transfer function must satisfy the reconstruction anti-aliasing conditions of the interpolator and decimator respectively, which in common case is the low-pass filter with the lower cut-off frequency of the two. Remember that optimization of the decimator is based on the idea not to filter samples that are immediately thrown away. Let’s extend this idea to not even create the samples at the interpolator that are thrown away in the decimator. See the table view of the resampling filter below.
There is made an elimination combined from decimation by the factor 3 and from interpolation by the factor 2. Since the resampling filter combines the two elimination rules it is even more complicated to be implemented in logic, one may need to draw even more tables to find the best strategy for target platform.
To be complete, we designed fast resampling filter as an FPGA IP core with which we can provide much stranger resampling factors with unprecedented harware demands. Contact us, if you want to know more or simply quote the IP core for your design.
- “Properly filtered” doesn’t necessarily mean low-pass filtering, the decimation may be wisely used for digital downconversion just by aliasing and selecting desired Nyquist zone.
- The half-band filter can be designed as even, but there won’t be the desired zeros in the impulse response.
- The filter needn’t necessarily to remove the dotted spectra. A designer can select any of the signal within the 1st nyquist zone and take advantage of an oversampling upconversion.
Leave a Reply