2D Sinusoidal Parameter Estimation with Offset Term

TL;DR

This paper addresses 2D sinusoidal parameter estimation including an offset term, deriving the maximum likelihood estimator and the Cramer-Rao lower bound to improve understanding of estimation accuracy.

Contribution

It introduces a novel model with an offset term for 2D sinusoidal estimation and derives the MLE and CRLB for this model.

Findings

Derived the MLE solution for the 2D sinusoid with offset

Established the CRLB for estimator variance

Provides theoretical bounds for estimation accuracy

Abstract

We consider the parameter estimation of a 2D sinusoid. Although sinusoidal parameter estimation has been extensively studied, our model differs from those examined in the available literature by the inclusion of an offset term. We derive both the maximum likelihood estimation (MLE) solution and the Cramer-Rao lower bound (CRLB) on the variance of the model's estimators.

Equations

$$f(x,y)=A\sin (2\pi (f_{0}x+f_{1}y)+\varphi )+B,$$

$$\frac{1}{N}\sum _{n=0}^{N-1}{e}^{i(\omega n+\varphi )}=\frac{e^{i(\pi /2-\omega /2+\varphi )}+e^{-i(\pi /2-\omega (N-1/2)-\varphi )}}{2N\sin (\omega /2)}$$

$$\lim_{N\to\infty}\frac{1}{N^{k+1}}\sum_{n=0}^{N-1}n^{k}e^{i(2\pi fn+\phi)}=\left\{\begin{array}[]{@{}c@{\quad}l@{}}\frac{1}{k+1}e^{i\phi}&\text{$f=0,1$}\\ 0&\text{$f\neq 0,1$}\\ \end{array}\right.$$

$$\text{var}(\widehat{{\theta }_i})\ge [{\mathbf{\eta }}^{-1}(\mathbf{\theta }){]}_{ii},$$

$$[\mathbf{\eta }(\mathbf{\theta }){]}_{ij}=-E\left[\frac{{\partial }^2\ln p(\mathbf{x};\mathbf{\theta })}{\partial {\theta }_i\partial {\theta }_j}\right]$$

$$s(x,y)=f(x,y)+w(x,y),x=0,\dots ,N-1;y=0,\dots ,N-1$$

$$\ln p(\mathbf{s};\mathbf{\theta })=$$

$$-2s(x,y)B+A^2{\sin }^2(2\pi (f_{0}x+f_{1}y)+\varphi )$$

$$+2AB\sin (2\pi (f_{0}x+f_{1}y)+\varphi )+B^2)$$

$$\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial A^2}=$$

$$=$$

$$-\frac{1}{2{\sigma }^2}\sum _x\sin (4\pi f_{0}x+2\varphi )\sum _y\sin (4\pi f_{1}y)$$

$$\approx$$

$$\boxed{E\left[\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial A^2}\right]\approx -\frac{N^2}{2{\sigma }^2}}$$

$$\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial A\partial B}$$

$$\boxed{E\left[\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial A\partial B}\right]\approx 0}$$

$$\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial A\partial \varphi }=$$

$$+2\cos (2\pi (f_{0}x+f_{1}y)+\varphi )(B-s(x,y)))$$

$$E\left[\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial A\partial \varphi }\right]=$$

$$+2\cos (2\pi (f_{0}x+f_{1}y)+\varphi )(B-E[s(x,y)]))$$

$$=$$

$$\approx$$

$$\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial A\partial f_0}=$$

$$+4\pi x\cos (2\pi (f_{0}x+f_{1}y)+\varphi )(B-s(x,y)))$$

$$E\left[\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial A\partial f_0}\right]=$$

$$+4\pi x\cos (2\pi (f_{0}x+f_{1}y)+\varphi )(B-E[s(x,y)]))$$

$$=$$

$$\approx$$

$$\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial {\varphi }^2}=$$

$$+2A\sin (2\pi (f_{0}x+f_{1}y)+\varphi )(s(x,y)-B))$$

$$E\left[\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial {\varphi }^2}\right]=$$

$$+2A\sin (2\pi (f_{0}x+f_{1}y)+\varphi )(E[s(x,y)]-B))$$

$$=$$

$$\approx$$

$$\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial \varphi \partial f_0}=$$

$$+4\pi Ax\sin (2\pi (f_{0}x+f_{1}y)+\varphi )(s(x,y)-B))$$

$$E\left[\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial \varphi \partial f_0}\right]=$$

$$+4\pi Ax\sin (2\pi (f_{0}x+f_{1}y)+\varphi )(E[s(x,y)]-B))$$

$$=$$

$$\approx$$

$$\frac{{\partial }^2\ln p(\mathbf{s};\mathbf{\theta })}{\partial f_0^2}=$$

Full text

I Introduction

In this paper, we examine the problem of parameter estimation of a 2D sinusoid. Although sinusoidal parameter estimation has been extensively studied [1, 2, 3, 4, 5], our model differs slightly from those examined in the available literature by the inclusion of an offset term. We derive both the MLE solution and the Cramer-Rao lower bound (CRLB) on the variance our model’s estimators, and then implement our approach on several fingerprint images of varying quality.

We specifically consider the discrete 2D sinusoidal signal

[TABLE]

where $x=0,...,N-1$, $y=0,...,N-1$; $\boldsymbol{\theta}=(A\;B\;\phi\;f_{0}\;f_{1})$ is the vector of parameters to be estimated: $A$ is the amplitude of the sinusoid, $B$ is its offset, $\phi$ is its phase shift, and ${\bf f}=(f_{1}\;\;f_{0})^{T}$ is its frequency. Such a model could describe the signal intensity at pixel $(x,y)$ of an $N\times N$ image. The main difference between Eq. (1) and those studied in [1, 2, 3, 4, 5] is the inclusion of the offset term $B$.

Eq. (1) arises in fingerprint biometrics. Fingerprint texture is characterized by the periodic flow of ridges and furrows, so it contains both frequency and orientation information; the frequency content is due to the inter-ridge spacing present in the fingerprint, while the orientation is due to the flow pattern exhibited by the ridges. If an acquired (gray-level) 2D fingerprint image is partitioned into sub-blocks, where each sub-block contains a ridge segment, the gray level intensity variations can be modeled via Eq. (1), whose parameters characterize the enclosed ridge’s frequency and orientation within the sub-block.

II Theory

The following theorems will prove useful in our derivations of both the CRLB and MLE of $\boldsymbol{\theta}$.

Lemma II.1

For $\omega\in[0,2\pi]$,

[TABLE]

Corollary II.1.1

For* $f\in[0,1]$ and integer $k\geq 0$,*

[TABLE]

II-A Crammer-Rao Lower Bound (CRLB) of Estimator $\boldsymbol{\theta}$

Consider the $p\times 1$ vector parameter $\boldsymbol{\theta}=(\theta_{1}\;\dots\;\theta_{p})$. We will assume that the estimator $\boldsymbol{\hat{\theta}}$ is unbiased. The CRLB gives a lower bound on the variance of any unbiased estimator, and the CRLB of estimator $\hat{\theta_{i}}$ is

[TABLE]

where $\boldsymbol{\eta(\theta)}$ is the $p$ x $p$ Fisher information matrix; it is defined as

[TABLE]

for $i=1,2,\dots,p$ and $j=1,2,\dots,p$.

We consider the signal

[TABLE]

where $f(x,y)$ is given by Eq. (1) and $w(x,y)$ is the noise. Since we assume the noise is white Gaussian, i.e. $w(x,y)=\frac{1}{2\pi\sigma^{2}}\exp(-\frac{x^{2}+y^{2}}{2\sigma^{2}})$, we have $s(x,y)\sim\mathcal{N}(f(x,y),\sigma^{2})$.

Denote ${\bf z}=\text{vec}\{f(x,y);\;x=0,\dots,N-1,\;y=0,\dots,N-1\}$ and ${\bf w}=\text{vec}\{w(x,y);\;x=0,\dots,N-1,\;y=0,\dots,N-1\}$; both are of dimension $N^{2}$ x $1$. Then Eq. (4) can be rewritten in vector form as ${\bf s}={\bf z}+{\bf w}$, where the signal measurements ${\bf s}\sim\mathcal{N}_{N^{2}}({\bf z}(\boldsymbol{\theta}),\sigma^{2}{\bf I}_{N^{2}\times N^{2}})$. Then the log-likelihood function of $\boldsymbol{\theta}$ (ignoring the fixed term) is

[TABLE]

We now derive the elements forming the Fisher information matrix, given by Eq. (3).

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial A^{2}}\right]$:

[TABLE]

where we have used the approximation that $\frac{1}{N}\sum_{x=0}^{N-1}\sin(4\pi f_{0}x+2\phi)\approx 0$ for large $N$ and $f_{0}\neq 0,1/2,1$.

[TABLE] 2. 2.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial A\partial B}\right]$:

[TABLE]

where we have employed the approximation that $\frac{1}{N}\sum_{x=0}^{N-1}\sin(2\pi f_{0}x+\phi)\approx 0$ for large $N$ and $f_{0}\neq 0,1$.

[TABLE] 3. 3.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial A\partial\phi}\right]$:

[TABLE]

[TABLE] 4. 4.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial A\partial f_{0}}\right]$:

[TABLE]

[TABLE]

where we have used the approximation that $\frac{1}{N^{2}}\sum_{x=0}^{N-1}x\sin(4\pi f_{0}x+2\phi)\approx 0$ for large $N$ and $f_{0}\neq 0,1/2$. 5. 5.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial A\partial f_{1}}\right]\approx\boxed{0}$ (calculation similar to (4)) 6. 6.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial B^{2}}\right]=\boxed{-\frac{N^{2}}{\sigma^{2}}}$ 7. 7.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial B\partial\phi}\right]\approx\boxed{0}$ (calculation similar to (2)) 8. 8.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial B\partial f_{0}}\right]\approx\boxed{0}$ (calculation similar to (4)) 9. 9.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial B\partial f_{1}}\right]\approx\boxed{0}$ (calculation similar to (4)) 10. 10.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial\phi^{2}}\right]$:

[TABLE]

[TABLE]

where we have employed the identity $\cos^{2}(x)=\frac{1+\cos(2x)}{2}$. 11. 11.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial\phi\partial f_{0}}\right]$:

[TABLE]

[TABLE] 12. 12.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial\phi\partial f_{1}}\right]\approx\boxed{-\frac{\pi A^{2}N^{2}(N-1)}{2\sigma^{2}}}$ (calculation similar to (11)) 13. 13.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial f_{0}^{2}}\right]$:

[TABLE]

[TABLE]

where we have used the approximation that $\frac{1}{N^{3}}\sum_{x=0}^{N-1}x^{2}\cos(4\pi f_{0}x+2\phi)\approx 0$ for large $N$ and $f_{0}\neq 0,1/2$ 14. 14.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial f_{0}\partial f_{1}}\right]$:

[TABLE]

[TABLE] 15. 15.

$E\left[\frac{\partial^{2}\ln p({\bf s};\boldsymbol{\theta})}{\partial f_{1}^{2}}\right]\approx\boxed{-\frac{\pi^{2}A^{2}N^{2}(N-1)(2N-1)}{3\sigma^{2}}}$ (calculation similar to (13))

Noting that the determinant is $|\boldsymbol{\eta}(\boldsymbol{\theta})|=\frac{\pi^{4}A^{6}N^{10}(N^{2}-1)^{2}}{144\sigma^{10}}$, matrix inversion yields

[TABLE]

Hence, the CRLB of our estimator $\boldsymbol{\hat{\theta}}$ in Eq. (1), under the assumption of white Gaussian noise $\mathcal{N}(0,\sigma^{2})$, is

[TABLE]

The CRLB of the amplitude and offset terms depend on known values, i.e. the dimension of image sub-block and the variance of the noise, while that of the frequencies and phase depend on an unknown parameter, i.e. the amplitude.

II-B Maximum Likelihood Estimation (MLE) of Sinusoidal Parameters

Recall that the log-likelihood function of $\boldsymbol{\theta}$ is

[TABLE]

In order to maximize the likelihood, we need to minimize the squared error:

[TABLE]

The estimator $\boldsymbol{\hat{\theta}}$ that minimizes the squared error $J$ is the maximum likelihood estimator.

We will return to the squared error, but let’s rewrite Eq. (1) as

[TABLE]

where $x=0,...,N-1;\;y=0,...N-1$. Let ${\bf u}$ and ${\bf v}$ be the $N^{2}\times 1$ vectors, respectively, denoting the array of $\sin$ and $\cos$ terms in Eq. (1). Denote $\alpha_{1}=A\cos\phi$, $\alpha_{2}=A\sin\phi$, and $\boldsymbol{\alpha}=(\alpha_{1}\;\alpha_{2}\;B)$. Further let $H=[{\bf u}\;{\bf v}\;{\bf 1}]$, which is $N^{2}\times 3$. Now we can rewrite the squared error as

[TABLE]

Optimizing $J$ with respect to $\boldsymbol{\alpha}$ yields

[TABLE]

so that

[TABLE]

Minimizing $J$ is now equivalent to maximizing ${\bf s}^{T}{\bf H}({\bf H}^{T}{\bf H})^{-1}{\bf H}^{T}{\bf s}$, or equivalently,

[TABLE]

Noting that

[TABLE]

and simplifying yields

[TABLE]

Recall that the Fourier transform (FT) of a function $f(x,y)$ is $F(f_{x},f_{y})=\sum_{x}\sum_{y}f(x,y)e^{-2\pi i(f_{x}x+f_{y}y)}$. Denoting the FT of $s(x,y)$ as $S(f_{0},f_{1})$, we lastly obtain

[TABLE]

where $|S(f_{0},f_{1})|^{2}$ denotes the periodogram of $s(x,y)$. Since the second term in Eq. (6) is fixed, the expression ${\bf s}^{T}{\bf H}({\bf H}^{T}{\bf H})^{-1}{\bf H}^{T}{\bf s}$ is maximized when the periodogram of the signal is maximized.

The frequencies at which the periodogram is maximized have to be found numerically. Denote the optimal frequencies as $\hat{f_{0}}$ and $\hat{f_{1}}$; now Eq. (5) becomes

[TABLE]

Hence, the maximum likelihood estimators of $\boldsymbol{\hat{\theta}}$ in Eq. (1) are

[TABLE]

Here, $F(f_{0},f_{1})=\sum_{x}\sum_{y}f(x,y)e^{-2\pi i(f_{0}x+f_{1}y)}$, i.e. the 2D discrete Fourier transform (FT) of Eq. (1), and $|F(f_{0},f_{1})|^{2}$ denotes the periodogram of $f(x,y)$. The frequencies at which the periodogram is maximized, $(\hat{f_{0}},\hat{f_{1}})$, have to be found numerically. Note that the maximum likelihood estimator of the offset term $B$ is simply the mean of the signal measurements, while the maximum likelihood estimator of the amplitude is the magnitude of the FT of the signal evaluated at the optimal frequencies.

III Size of $N$

To get an idea of how big the dimension $N$ should be in order for the approximation

[TABLE]

to hold, we look at the plots of two functions: 1) $y(f)=\frac{1}{N}\sum_{x=0}^{N-1}\sin(4\pi fx+\phi)$ and 2) $y(f)=\frac{1}{N}\sum_{x=0}^{N-1}\sin(2\pi fx+\phi)$ for $\phi=0,\pi/4$ and $N=20$ measurements. These plots are shown in Figs. 1 and 2.

For the $\omega=4\pi f$ case, shown in Fig. 1, if $f$ is not near 0, 0.5, or 1, the summation is approximately zero. For the $\omega=2\pi f$ case, shown in Fig. 2, if $f$ is not near 0 or 0.5, the summation is approximately zero; however, it has a slower approximation to zero than the $\omega=4\pi$ case. The plots illustrate that $N=20$ measurements is adequate for Eq. (7) to be valid.