The folowing list contains naming conventions used in this document and Mulirate's source code. You will see references to columns and rows. Take note that Julia is column major language.
h: a vector of filter coefficients (aka taps)PFB: poly phase filter bank<symbol>Len: a quantity of something𝜙: a phase of a PFBƒ: frequency<symbol>Time: vector of time values associated with another vector A conventional polyphase FIR filter implementation takes filter taps and splits them up into phases. Why are they called phases? I still don't know. But the important thing is that this allows for efficient interpolation. In the naive approach to interpolation by a factor of L one inserts, or 'stuffs', L-1 zeros between each input sample prior to lowpass filtering:
In [78]:
interpolation = 4
xLen = 4
x = [1:xLen]
xStuffed = vec(vcat( zeros(Int,3,4), x' ))
println( "x = $(x')" )
println( "xStuffed = $(xStuffed)" )
This inefficient however. The zeros that we stuffed inbetween the input samples contribute nothing to the output of the filtering process other than placeholders to keep filter taps properly aligned with input. Remeber that FIR filter involves convolving x and h. Let's compute, using strings for readablity, the first four output samples with hypothetical vector of taps h:
y[1] = x1*h12 + x2*h11 + x3*h10 + x4*h9 + x5*h8 + x6*h7 + x7*h6 + x8*h5 + x9*h4 + x10*h3 + x11*h2 + x12*h1
y[2] = x2*h12 + x3*h11 + x4*h10 + x5*h9 + x6*h8 + x7*h7 + x8*h6 + x9*h5 + x10*h4 + x11*h3 + x12*h2 + x13*h1
y[3] = x3*h12 + x4*h11 + x5*h10 + x6*h9 + x7*h8 + x8*h7 + x9*h6 + x10*h5 + x11*h4 + x12*h3 + x13*h2 + x14*h1
y[4] = x4*h12 + x5*h11 + x6*h10 + x7*h9 + x8*h8 + x9*h7 + x10*h6 + x11*h5 + x12*h4 + x13*h3 + x14*h2 + x15*h1
In [79]:
hLen = 12
for yIdx in 1:4
hIdx = hLen
print( "y[$yIdx] = " )
xStartIdx = yIdx
xStopIdx = xStartIdx+hLen-1
for xIdx in xStartIdx:xStopIdx
# print( "$(xStuffed[xIdx])*h$hIdx" )
print( "x$xIdx*h$hIdx" )
xIdx < xStopIdx && print( " + " )
hIdx -= 1
end
println()
end
See a pettern? Let's make it more obvious by skipping miltiplications where x == 0:
In [80]:
for yIdx in 1:4
xStartIdx = yIdx
xStopIdx = xStartIdx+hLen-1
hIdx = hLen
xPrevZero = true
for xIdx in xStartIdx:xStopIdx
if xStuffed[xIdx] != 0
xPrevZero == false && print( " + " )
print( "$(xStuffed[xIdx])*h$hIdx" )
xPrevZero = false
end
hIdx -= 1
end
println()
end
Now the pattern is more obious. For the first four y's we can decompose h into four phases, 𝜙1 = [ h1, h5, h9 ], 𝜙2 = [ h2, h6, h10 ], 𝜙3 = [ h4, h7, h11 ], and 𝜙4 = [ h4, h8, h12 ].
In [81]:
using Multirate
using DSP
using PyPlot
In [82]:
N𝜙 = 32 # Number of polyphase partitions
tapsPer𝜙 = 10 # N𝜙 * tapsPer𝜙 will be the length of out prototype filter taps
resampleRate = float64(pi) # Can be any arbitrary resampling rate
polyorder = 4 # Our taps will tranformed into
Nx = 40 # Number of signal samples
t = 0:Nx-1 # Time range
xƒ1 = 0.15 # First singal frequency
xƒ2 = 0.3 # Second signal frequency
x = cos(2*pi*xƒ1*t) .+ 0.5sin(2*pi*xƒ2*t*pi) # Create the two signals and add them
tx = [0:length(x)-1] # Signal time vector
cutoffFreq = min( 0.45/N𝜙, resampleRate/N𝜙 ) # N𝜙 is also the integer interpolation, so set cutoff frequency accordingly
hLen = tapsPer𝜙*N𝜙 # Total number of filter taps
h = firdes( hLen, cutoffFreq, DSP.kaiser ) .* N𝜙 # Generate filter taps and scale by polyphase interpolation (N𝜙)
myfilter = FIRFilter( h, resampleRate, N𝜙, polyorder ) # Construct a FIRFilter{FIRFarrow} object
y = filt( myfilter, x ) # Filter x
ty = [0:length(y)-1]./resampleRate - tapsPer𝜙/2 # Create y time vector. Accout for filter delay so the plots line up
Out[82]:
In [ ]: