5
$\begingroup$

I have a signal in white noise that has the following form: \begin{equation} r[t_i] = A e^{-t_i/\tau} \sin{(\omega t_i + \phi)} + n[t_i] \end{equation}

I would like to test whether the signal (1st term) at a known frequency $\omega$ is present against the null-hypothesis of pure noise. I know for the noise term (Poissonian with $N>>1$) the standard deviation. $\phi$ and $\tau$ are unknown, however $\tau$ should be long compared to $t_\mathrm{max}$.

Is there a statistical test for this purpose? (If there are already implemented tests, Python is preferred but any language is fine)

I simulate a signal e.g. by

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(123)

# Signal parameters
f   = 0.25 # Known
tau = 50   # Unknown
phi = 0    # Unknown
A   = 0.07 # Unknown A<1

# Signal modeling
t = np.arange(25)
s = A*np.sin(2*np.pi*f*t+phi)*np.exp(-t/tau)


# Noise parameters:
N   = 20e3 # Mean number of counts (known)
c   = 0.3  # Contrast (known)

# Noise modeling
cts = N *(c-2-2*c*s)/(c-2)
r = np.random.poisson(cts) 
plt.plot(r)
plt.title('signal with noise')
plt.show()

Exemplatory signal

$\endgroup$
3
  • $\begingroup$ Can you post some of the data, or some reasonable simulation, so we can better understand what you're working with? $\endgroup$ Feb 18 at 20:01
  • $\begingroup$ I edited the question. The simulated signal as fairly good SNR such that it is easier to understand what I am hunting for. $\endgroup$
    – P. Egli
    Feb 18 at 20:32
  • $\begingroup$ 1. Your terminology seems very strange to this statistician; in particular the use of the term "noise" for the model of $n$ seems odd there. It looks like it's not centered at 0, so it's more a bias on your damped sinusiodal "model", rather than mere noise. Note that the rest of the model (aside $n$) converges to 0, but your values are cycling around 20000. I sure wouldn't refer to that as just noise. 2. Maybe I misunderstood something but if $\tau$ is "long compared to $t_\max$", then wouldn't you only have a fraction of a cycle rather than many? $\endgroup$
    – Glen_b
    Feb 18 at 22:13

3 Answers 3

5
$\begingroup$

My intuition tells me:

  1. build an optimal filter using e.g. a correlation with your expected signal r(t_i): see first example under https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.correlate.html
  2. Run this filter over 1000 or more random noise signals, and observe the maxima of the correlation signals.
  3. Then do it for your signal. If this maximum exceeds the 95% quantile of the observed maxima under the null hypothesis, then reject the null hypothesis with significance level 5%

For why this is sound see Rogers post which relates statistical signal processing and matched filters (optimal filters).

$\endgroup$
4
$\begingroup$

It may be possible to use Bayesian statistics here.

Let's set testing aside for a moment and instead focus on estimation. If we can construct a model and assess the estimates of $A$, maybe we can say something substantive about the decay term.

A simple model in Stan would be

data {
    int n;
    array[n] int r;
    vector[n] t;
    real c;
    real N;
    real f;
}
parameters {
    real<lower=0> inv_tau;
    real<lower=-1, upper = 1> phi_scaled;
    real<lower=0, upper=1> A;
}
transformed parameters {
    vector[n] s = A .* sin(2 * pi() * f * t +  pi() * phi_scaled) .* exp(-t * inv_tau);
    vector[n] cts = N .*(c-2-2*c*s)./(c-2);
}
model{
    A ~ uniform(0, 1);
    phi_scaled ~ uniform(-1, 1);
    inv_tau ~ cauchy(0, 1);

    r ~ poisson(cts);
}
generated quantities {
   real estimated_tau = 1.0 / inv_tau;
   real estimated_phu = pi() * phi;
   
}

Here, I am modelling scaled_phi to be on the interval $[-1, 1]$ and simply multiplying it by $\pi$ in the signal part. Additionally, I am modelling the inverse of $\tau$ rather than tau because I typically find that small parameters are often modelled more easily than large ones.

When I run the model, I get good diagnostics (note, you will need to install cmdstan and cmdstanpy if you want to replicate this).


import numpy as np
import cmdstanpy
import matplotlib.pyplot as plt
import arviz as az
def simulate_data():
    rng = np.random.RandomState(123)

    # Signal parameters
    f   = 0.25 # Known
    tau = 50   # Unknown
    phi = 0    # Unknown
    A   = 0.07 # Unknown A<1

    # Signal modeling
    t = np.arange(25)
    s = A*np.sin(2*np.pi*f*t+phi)*np.exp(-t/tau)


    # Noise parameters:
    N   = 20e3 # Mean number of counts (known)
    c   = 0.3  # Contrast (known)

    # Noise modeling
    cts = N *(c-2-2*c*s)/(c-2)
    r = rng.poisson(cts) 


    return t, r, cts


def fit_model(*args, **kwargs):
    t, r, cts = simulate_data()

    model = cmdstanpy.CmdStanModel(stan_file="signal.stan")

    data = dict(t=t.tolist(), r=r.tolist(), n=r.size, f = 0.25, N = 20000, c = 0.3)

    fit = model.sample(data, **kwargs)

    return fit

fit = fit_model()

print(fit.diagnose())

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

Now, we can look at some plots of the estimated parameters

var_names = ['A', 'estimated_tau', 'estimated_phi']

fig, ax = plt.subplots(dpi = 240, ncols=3, nrows=1, figsize = (15, 5))

for vn, a in zip(var_names, ax):
    x = fit.stan_variable(vn)
    a.hist(x, edgecolor='white')
    a.set_title(vn)

enter image description here

$A$ and $\phi$ are very well estimated. Tau remains uncertain (95% credible interval: 21 - 389), perhaps that could be combatted with some better priors.

From these plots, we could conclude that $A>0$ and hence the first term is likely present.

$\endgroup$
3
$\begingroup$

In engineering context likelihood ratio test is a typical approach for signal detection. The Poissonian noise for $N\gg 1$ can be viewed as Gaussian, which allows to write the joint probability density for the observations at the selected time moments $t_1,..., t_K$. This probability could be compared to the probability that the observations are purely due to the noise. One the compares the ratio of the two probabilities to a specified threshold to decide whether the signal is present or not.

The detailed derivations can be likely found in older (pre-digital era) radio-engineering or physics textbooks. After a bit of googling I found these notes on Statistical signal processing, which seem to be pretty close to what the OP is asking.

$\endgroup$
3
  • 1
    $\begingroup$ +1 the notes talk about ”matched filter” or in other words an optimal filter which is related to my post. Nice to see that there is sound background $\endgroup$
    – Ggjj11
    Feb 19 at 11:51
  • 1
    $\begingroup$ PS: en.wikipedia.org/wiki/… also shows this connection (although less explicit than in the notes :) ) $\endgroup$
    – Ggjj11
    Feb 20 at 12:32
  • $\begingroup$ @Ggjj11 there is likely parallel terminology here, depending on whether one speaks in the context of signal processing or random/stochastic processes - even though one deals with the same problem (signal detection). I learned both in the context of detecting radio-waves (awhile ago, during my undergrad.) $\endgroup$
    – Roger V.
    Feb 20 at 12:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.