Documentation |
On this page… |
---|
Multirate Filter Bank Implementation |
The upfirdn function alters the sampling rate of a signal by an integer ratio P/Q. It computes the result of a cascade of three systems that performs the following tasks:
Upsampling (zero insertion) by integer factor p
Filtering by FIR filter h
Downsampling by integer factor q
For example, to change the sample rate of a signal from 44.1 kHz to 48 kHz, we first find the smallest integer conversion ratio p/q. Set
d = gcd(48000,44100); p = 48000/d; q = 44100/d;
In this example, p = 160 and q = 147. Sample rate conversion is then accomplished by typing
y = upfirdn(x,h,p,q)
This cascade of operations is implemented in an efficient manner using polyphase filtering techniques, and it is a central concept of multirate filtering. Note that the quality of the resampling result relies on the quality of the FIR filter h.
Filter banks may be implemented using upfirdn by allowing the filter h to be a matrix, with one FIR filter per column. A signal vector is passed independently through each FIR filter, resulting in a matrix of output signals.
Other functions that perform multirate filtering (with fixed filter) include resample, interp, and decimate.
In the case of FIR filters, it is possible to design linear phase filters that, when applied to data (using filter or conv), simply delay the output by a fixed number of samples. For IIR filters, however, the phase distortion is usually highly nonlinear. The filtfilt function uses the information in the signal at points before and after the current point, in essence "looking into the future," to eliminate phase distortion.
To see how filtfilt does this, recall that if the z-transform of a real sequence x(n) is X(z), the z-transform of the time reversed sequence x(n) is X(1/z). Consider the processing scheme.
Image of Anti Causal Zero Phase Filter
When |z| = 1, that is z = e^{j}^{ω}, the output reduces to X(e^{j}^{ω})|H(e^{j}^{ω})|^{2}. Given all the samples of the sequence x(n), a doubly filtered version of x that has zero-phase distortion is possible.
For example, a 1-second duration signal sampled at 100 Hz, composed of two sinusoidal components at 3 Hz and 40 Hz, is
fs = 100; t = 0:1/fs:1; x = sin(2*pi*t*3)+.25*sin(2*pi*t*40);
Now create a 10-point averaging FIR filter, and filter x using both filter and filtfilt for comparison:
b = ones(1,10)/10; % 10 point averaging filter y = filtfilt(b,1,x); % Noncausal filtering yy = filter(b,1,x); % Normal filtering plot(t,x,t,y,'--',t,yy,':')
Both filtered versions eliminate the 40 Hz sinusoid evident in the original, solid line. The plot also shows how filter and filtfilt differ; the dashed (filtfilt) line is in phase with the original 3 Hz sinusoid, while the dotted (filter) line is delayed by about five samples. Also, the amplitude of the dashed line is smaller due to the magnitude squared effects of filtfilt.
filtfilt reduces filter startup transients by carefully choosing initial conditions, and by prepending onto the input sequence a short, reflected piece of the input sequence. For best results, make sure the sequence you are filtering has length at least three times the filter order and tapers to zero on both edges.
Duality between the time domain and the frequency domain makes it possible to perform any operation in either domain. Usually one domain or the other is more convenient for a particular operation, but you can always accomplish a given operation in either domain.
To implement general IIR filtering in the frequency domain, multiply the discrete Fourier transform (DFT) of the input sequence with the quotient of the DFT of the filter:
n = length(x); y = ifft(fft(x).*fft(b,n)./fft(a,n));
This computes results that are identical to filter, but with different startup transients (edge effects). For long sequences, this computation is very inefficient because of the large zero-padded FFT operations on the filter coefficients, and because the FFT algorithm becomes less efficient as the number of points n increases.
For FIR filters, however, it is possible to break longer sequences into shorter, computationally efficient FFT lengths. The function
y = fftfilt(b,x)
uses the overlap add method to filter a long sequence with multiple medium-length FFTs. Its output is equivalent to filter(b,1,x).