PolyInterp.jl

import IPPDSP

function polyize{T}( h::Vector{T}, interpolation )
    hLen         = length( h )
    tapsPerPhase = int( ceil( hLen/interpolation ))
    pfbSize = tapsPerPhase * interpolation
    # check that the vector is an integer multiple of interpolation
    if hLen != pfbSize 
        hExtended             = similar( h, pfbSize )
        hExtended[1:hLen]     = h
        hExtended[hLen+1:end] = 0
        h                     = hExtended
    end
    nFilters     = interpolation
    hLen         = length( h )
    tapsPerPhase = int( hLen/nFilters )
    pfb          = reshape( h, nFilters, tapsPerPhase )'
end

function upsample{T}( x::Vector{T}, interpolation )
    buffer = zeros(T, length(x) * interpolation )
    
    for n = 1:length(x)
        buffer[ n*interpolation - interpolation + 1 ] = x[n]
    end
    
    buffer
end

function interpolate{T}( h::Vector{T}, x::Vector{T}, interpolation )
    x       = upsample( x, interpolation )
    convAns = conv( h, x )[1:end-length(h)+1]
end

function polyinterpolate{T}( pfb::Array{T}, x::Vector{T} )
    pfbSize            = size( pfb )
    numTaps            = pfbSize[1]  # numTaps/filter is the number of rows
    numFilter          = pfbSize[2]  # columns hold individual filters for the bank
    interpolationRatio = numFilter
    inputLenth         = length( x )
    outputBuffer       = similar( x, inputLenth * interpolationRatio ) # output interpRatio times larger than the inupt vector
    intermediateBuffer = similar( x, interpolationRatio ) # to store the output of each

    for n = 1:numTaps-1
        for f = 1:numFilter
            intermediate = zero(T)
            for m = 1:n                
                @inbounds intermediate += pfb[m,f] * x[n-m+1]
            end
            @inbounds outputBuffer[n*numFilter-numFilter+f] = intermediate
        end
    end
   
    for n = numTaps:inputLenth
        for f = 1:numFilter
            intermediate = zero(T)
            for m = 1:numTaps
                @inbounds intermediate += pfb[m,f] * x[n-m+1]
            end
            @inbounds outputBuffer[n*numFilter-numFilter+f] = intermediate
        end
    end

    outputBuffer
end

polyinterpolate( h, x, interpolation ) = polyinterpolate( polyize( h, interpolation ), x )



function test( numSamples, numTaps, interpolation )
    interpolation = 4
    h     = IPPDSP.lowpass( Float64, numTaps, 0.5/interpolation, true )
    
    i = 17
    
    for i = numSamples
        nx    = i
        x     = ones(nx) 
        
            
        t = time()
        polyAns = polyinterpolate( h, x, interpolation )
        tPoly = time()-t
    
        # t = time()
        # convAns = interpolate( h, x, interpolation )
        # tConv = time()-t
        
        
        t = time()
        ippAns = IPPDSP.filt( h, x, interpolation, 1 )
        tIPP = time()-t
        
        
        println( i, ":\tPoly time: ", tPoly, " IPP time: ", tIPP )

        # for i = 1:nx*interpolation
            # isapprox( polyAns[i], convAns[i] ) || error( string( "Something when wrong at index ", i ))
        # end

        # if tPoly < tConv
            # break
        # end

    end
end


#=
interpolation = 10
nx            = 2^16 #100_000
nt            = 30

x  = rand(nx)
h  = IPPDSP.lowpass( Float64, nt, 0.5/interpolation, true )

gc()
t = time()
polyAns = polyinterpolate( h, x, interpolation )
tPoly = time()-t

gc()
t = time()
ippAns = IPPDSP.filt( h, x, interpolation, 1 )
tIPP = time()-t


[ polyAns ippAns ]

for i = 1:nx*interpolation
    isapprox( polyAns[i], ippAns[i] ) || error( string( "Something when wrong at index ", i ))
end

@printf( "Poly Time: %0.3e\tIPP Time: %0.3e",  tPoly, tIPP )

=#