More on Frequency Constrained Techniques

Deriving adaptive filters within a given frequency range was demonstrated in [1]. The goal was to constrain the filter to only operate within the passband of an ADC. Outside the passband was effectively a don’t-care specification. This technique can be pushed to further to optimize various signal model parameters. While the algorithms it leads to are much less efficient than existing algorithms, it leads to the possibility of estimating parameters only within certain frequency bands. This might be valuable, for example, if there was a known corruption or interference.

Inverse Filter

From [1] and taking the complex gradient [2]:

W[n] = 2ωp, if n = 0; 2sin(ωpn)/n, if n

Rx[k,l] =

dxy[k] =

Rxh = dxy

These equations were used to derive a frequency-constrained adaptive filter. However, they can also be used to find an inverse filter. Replace x with an impulse response and y with an appropriately delayed delta. There is sample code in Appendix A. In the unconstrained version of the Wiener-Hopf equations, a regularization parameter can be used to dampen the effect of noise boosting properties of inverse filter. But it is not very effective when the filter has deep nulls. By using a combination of the regularization parameter and frequency range, there is more control over the design of the inverse filter. This is demonstrated in figure 1. The filter is an ISI filter with some corruption in the passband. With appropriately chosen parameters, the passband is restored without increasing noise in the stopband.

Figure 1: (A) Filter with corruption in passband. (B) Inverse filter with -20 dB regularization. (C) Inverse filter with -200 dB regularization. (D) Inverse filter with -20 dB regularization and constrained to half the frequency range. (E) Convolution of A and D.

Center Frequency and Symbol Timing Estimation

When the signal constellation is known, there are a number of ways to estimate the phase of the sampler so that it aligns with the symbol. One straightforward method is to use steepest descent on center frequency and delay as a function of a Farrow Filter [3]. X and F are the data matrix and Farrow Filter of length M, respectively. W is the frequency translation matrix; the Hadamard product is used to translate the rows of X. C[.] is a function mapping the argument to the nearest constellation point. It’s assumed constant when taking the derivative. V is a window (square root hanning seems to work well).

X = [xM-1, …, x1, x0

xM, …, x2, x1

....

xN+M-1, …, xN+1, xN]

W(ω) =

[ejω(M-1), … , ejω1, ejω0

ejωM, … , ejω2, ejω1

….

ejω(N+M-1), … , ejω(N+1), ejωN]

dW(ω)/dω =

[j(M-1)ejω(M-1), … , j1ejω1, j0ejω0

jMejωM, … , j2ejω2, j1ejω1

….

j(N+M-1)ejω(N+M-1), … , j(N+1)ejω(N+1), jNejωN]

The update equations follow. The update index is the iteration number, but it could also be replaced with the sample number.

ωi = ωi-1 - 𝜇ωRe

𝛼i = 𝛼i-1 - 𝜇𝛼Re

This optimization can also be done in the frequency domain. As before, the frequency variable needs to be integrated out so that the result is a function of only ω and 𝛼. The center frequency variable is now u since ω is needed for the DTFT. X is X delayed either by Farrow Filter or some other fractional delay method; it is needed only to map each xi to its nearest constellation point.

X = [x0, …, xN-2, xN-1]

W(u) = [eju0, … , eju(N-2), eju(N-1)]

Ω(ω) = [e-jω0, …, e-jω(N-2), e-jω(N-1)]

It’s the same procedure to find the update equations, except each term needs to be integrated first. Lowercase letters with subscripts are elements of the vectors represented by uppercase.

ui = ui-1 - 𝜇uRe

𝛼i = 𝛼i-1 - 𝜇𝛼Re

A Two Dimensional Example

Video compression requires computing motion vectors of macroblocks across frames. Sub-pixel motion vectors can improve compression at a given quality. One obvious method is to upsample the image. It can also be done by direct minimization in frequency. The frequency variables are ωx and ωy; the macroblocks are M1 and M2.

There are a lot of terms, but due to the separability, the result is much the same as the one dimensional update equations. Because there are so many terms, it’s written compactly in terms of the matrices Sx and Sy, which are functions of m and n. The summation without limits is used to sum over all terms.

Sx,k,l(m,n) = (ωpcos(ωp(l-n+𝛼x)-sin(ωp(l-n+𝛼x)))/(l-n+𝛼x)2*sin(ωp(k-m+𝛼y))/(k-m+𝛼y)

Sy,k,l(m,n) = (ωpcos(ωp(l-n+𝛼y)-sin(ωp(l-n+𝛼y)))/(l-n+𝛼y)2*sin(ωp(k-m+𝛼x))/(k-m+𝛼x)

𝛼x,i = 𝛼x,i-1 + 𝜇𝛼xRe

𝛼y,i = 𝛼y,i-1 + 𝜇𝛼yRe

In figure 2, a 32x32 block is taken from the corner of the Mandrill’s eye. It’s shifted ¼ pixel down and ½ pixel to the right. The horizontal and vertical offset estimates are shown.

Figrue 2: (Top) Mandrill with extracted macroblock. (Bottom) Estimates of pixel offsets. The step size is .0001 and ωp=𝜋/2.

Appendix A: Inverse Filter Code

function get_constrained_inverse(h,N,ωp,noise_dB)

M = length(h)

σ²ᵥ = 10^(noise_dB/10)

x = [zeros(M-1); h; zeros(N-M+1)]

y = zeros(Complex128,N)

y[M] = 1.

n = N-1:-1:-N+1

b = zeros(Complex128,M)

A = zeros(Complex128,M,M)

for i = 1:M

for k = 1:N

c = 2.0./n[N-k+1:2N-k].*sin.(ωp*n[N-k+1:2N-k])

c[k] = 2.0ωp

b[M-i+1] += conj(x[k+i-1])*sum(y[1:N].*c)

end

for i = 1:M

for k = 1:N

for j = 1:N

if k==j

c = 2.0ωp

else

c = 2.0/(j-k)*sin.(ωp*(j-k))

end

A[M-i+1,M:-1:1] += conj(x[j+i-1])x[k:k+M-1]*c

end

A += eye(M)*σ²ᵥ

A\b

end

References

[1] https://drive.google.com/open?id=1ktpy_ZDrcPVVEA373AqYvAZSRFqycCbd

[2] https://mediatum.ub.tum.de/doc/631019/631019.pdf

[3] http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=15483