Measure Sine Wave Amplitude from ADC Signal

Question

Question

I'm trying to do some DSP that I have never done before and a nudge into the right direction would come in handy. The context is replication of this project.

The system architecture is the following:

An excitation signal is generated by a DAC. I use a 10 sample sine table to generate a 200 kHz sine wave which is subsequently filtered and the resulting SINAD is about 60 dB. The signal is then connected to the eight port passive system (see project description) via two analogue multiplexers, one carrying the excitation signal and the other the reference GND. All eight ports are connected to a 12 bit (at least 10 ENOB) ADC. The DC component is centred within the ADCs FSR by the input filter. The DAC and the ADCs are integrated into a microcomputer (Cortex-M4F).

The information I'm interested in is in the incoming signal's amplitude.

My current solution is the following:

• sample the signal with 1.25 Msps and collect 256 samples (about 40 periods of the signal with about 6 samples per period)
• estimate the signals mean and MS (mean square, like RMS but without applying the sqrt) using the 256 samples and the common equations
• calculate the MS of the AC component with $$ U_{RMS} = \sqrt{U_{DC}^2 + U_{AC,RMS}^2} \Leftrightarrow ()^2\\ U_{MS} = U_{DC}^2 + U_{AC,MS} \Leftrightarrow\\ U_{AC,MS} = U_{MS} - U_{DC}^2 $$

I'm happy with the MS of the AC component because it incorporates the information, just like the amplitude of the signal because of their direct relation.

I've read a few articles, especially this one Precision measurement of sine wave amplitude with ADC which sounded like it was exactly my problem. I couldn't wrap my head around the accepted solution though and thus asking this question.

I'm not very experienced in DSP and don't know efficient algorithms to make the most of my simple setup. What would you do to measure the signals amplitude, especially in the context of maximising the result's ENOB under the condition of 256 samples?

I captured some ADC data. Have a look at it here to get an impression of the signal's properties.

The steps in the signal amplitude hail from switching the system excitation and measurement ports to a different state. I included the switching action because it gives an impression of the signal conditions in all relevant cases, amplitude-wise. The sections with no visible signal aren't relevant, that's the case when the specific ADC channel is connected to the reference ground.

Implementation

I started with the implementation of Olli's and Cedron's idea today. Let me show you the code I wrote.

First I generated the sine and cosine tables for the DFT. First try, I didn't choose a period of 45 samples, as Olli suggested, but the theoretical length of

$$ \frac{samples}{period} = \frac{f_{sample}}{f_{signal}} = \frac{1.25~MHz}{200~kHz} = \frac{25}{4}\\ N_{samples} = \frac{25}{4} \cdot N_{periods} \mid N_{samples}, N_{periods} \in \mathbb{N}^+\\ \Rightarrow N_{periods} = 4 \cdot k \mid k \in \mathbb{N}^+ $$

which means that my tables need to be 25 entries long. I generated them with Matlab.

s1_15 = quantizer('fixed', 'round', 'saturate', [16, 15]);

x = linspace(0, 4 * 2 * pi, 25 + 1);
x(end) = [];

sin_table = num2int(s1_15, sin(x));
cos_table = num2int(s1_15, cos(x));

My previous sample buffer was 256 deep, the new sample buffer holds 250, the next smallest length that meets the requirement above.

// ch1 and ch2 alternating
u16 adc_data[500];

// sine and cosine table for single frequency DFT calculation
// s1.15 number format
const s16 dft_sin_table[25] = {0, 27667, 29649, 4107, -25248, -31164, -8149, 22431, 32188, 12063, -19261, -32703, -15786, 15786, 32703, 19261, -12063, -32188, -22431, 8149, 31164, 25248, -4107, -29649, -27667};

const s16 dft_cos_table[25] = {32767, 17558, -13952, -32510, -20887, 10126, 31739, 23887, -6140, -30467, -26510, 2058, 28715, 28715, 2058, -26510, -30467, -6140, 23887, 31739, 10126, -20887, -32510, -13952, 17558};

void evaluation(void) {
    // Accumulators for dot product
    // s21.15 number format
    s64 dot_prod_sin[2] = {0};
    s64 dot_prod_cos[2] = {0};

    // Index for sine and cosine lookup
    uint table_index = 0;

    // for each sample in sample buffer
    for (uint it_sample = 0; it_sample <= 500 - 1; it_sample += 2) {
        // First channel

        // sample in u12.0 number format
        // multiplication result in s13.15 number format
        dot_prod_sin[0] += adc_data[it_sample] * dft_sin_table[table_index];
        dot_prod_cos[0] += adc_data[it_sample] * dft_cos_table[table_index];

        // Second channel

        // sample in u12.0 number format
        // multiplication result in s13.15 number format
        dot_prod_sin[1] += adc_data[it_sample + 1] * dft_sin_table[table_index];
        dot_prod_cos[1] += adc_data[it_sample + 1] * dft_cos_table[table_index];

        // Increment and wrap table index
        table_index += 1;
        if (table_index > 24) {
            table_index = 0;
        }
    }
}

The first results looked pretty good and where lots more expressive than my previous solution. The computational complexity was totally acceptable and worked for me.

Answer 1

In a way this does not directly answer your question, because I will be suggesting that you change the section length from 256 to 225 (or 270 if you can afford it) to do coherent sampling.

Data analysis

Within each 256-sample section, your data seems to be periodic with period 45. Because the oscillator and the digital-to-analog converter (DAC) clocks originate from the same master clock via frequency dividers and/or multipliers, the discrete-time signal is bound to be periodic. That is, if we do not take into account possible time-varying qualities of the non-digital part of the system. That the frequency of interest is a known constant is a great advantage, because then only the amplitude of the oscillation at that frequency needs to be detected, not the frequency itself.

Figure 1. First 256 samples of channel 1. The horizontal axis grid is distributed every 45 samples. Lines connect the samples for easier counting of oscillation cycles of which there appear to be $7$ per period. The continuous-time frequency of interest is 194.444... kHz, $7/45$ times the sampling frequency, a bit off from the specified 200 kHz.

Fourier method

If the apparent discrete-time period $P = 45$ is true or very nearly true, it will be advantageous to do the analysis over $N = PL$ samples, where $L$ is integer. This analysis length sets up discrete Fourier transform (DFT) bins at 1) the negative frequency of interest, 2) the positive frequency of interest, and 3) zero frequency. A sinusoid consists of both a negative-frequency and a positive-frequency complex exponential: $\cos(x)$ $=$ $\tfrac{1}{2}e^{-ix}$ $+$ $\tfrac{1}{2}e^{ix}.$ The complex exponential representation is convenient, because the magnitude of either the positive or the negative frequency DFT bin directly gives half the mid-to-top amplitude of the sinusoid. The $N$ complex exponential basis functions of DFT are orthogonal, that is to say uncorrelated. So any one of the three listed input signal components present can be detected with no disturbance from the presence of the other two.

Setting up the positive frequency of interest as $\omega = 2\pi\times7/45$ and the analysis window length as $N = 45\times5 = 225,$ the value of the DFT bin of interest is calculated from data $x$ as:

$$X_\omega = \frac{1}{N}\sum_{k=0}^{N-1}e^{-i\omega k}x[k] = \frac{1}{N}\sum_{k=0}^{N-1}\big(\cos(\omega k) - i\sin(\omega k)\big)x[k]$$

If $x$ was just a sinusoid (cosine) of mid-to-peak amplitude $A$ and initial phase $\theta$ we would obtain:

$$X_\omega = \frac{1}{N}\sum_{k=0}^{N-1}e^{-i\omega k}A\cos(\omega k + \theta) = \frac{Ae^{i\theta}}{2} + \frac{Ae^{-i(2\theta - 2\omega + \pi)/2} - Ae^{-i(4N\omega + 2\theta - 2\omega + \pi)/2}}{4N\sin(\omega)}$$

If $A\cos(\omega k)$ is $2N\text{-}$periodic, then the last term above vanishes. In coherent sampling $A\cos(\omega k)$ is $N\text{-}$periodic and therefore also $2N\text{-}$periodic, in which case we have:

$$|X_\omega| = \frac{A}{2}\quad\Rightarrow\quad A = 2|X_\omega|$$

The estimated mean square power of the sinusoid at the same frequency is:

$$2|X_\omega|^2$$

This is related to the estimated mid-to-peak amplitude of the sinusoid by:

$$A = 2\times\sqrt{|X_\omega|^2} = \tfrac{2}{\sqrt{2}}\times\sqrt{2|X_\omega|^2} = \sqrt{2}\times\sqrt{2|X_\omega|^2}$$

where $\sqrt{2}$ is the crest factor of sine wave. The squared magnitude of a complex number $C$ can be calculated by $|C|^2 = \operatorname{Re}(C)^2 + \operatorname{Im}(C)^2.$

Depending on how you use the obtained results, often you can do what you want with an intermediate result. In that case it may not be needed to calculate the normalization by $\frac{1}{N}$ or $\frac{2}{N},$ to multiply by $\sqrt{2},$ or to calculate the square root to get from squared magnitude to magnitude.

You can do the estimation directly using memory-stored complex exponential or cosine and sine coefficients (like @CedronDawg), using the Goertzel algorithm, or by calculating a full DFT (not recommended).

With coherent sampling, the sinusoid mean square power estimate $2|X_\omega|^2$ is equivalent to IEEE 1057 3-parameter estimate (known frequency, unknown sinusoid phase, sinusoid amplitude and noise level), because of the orthogonality of the DFT basis vectors. $2|X_\omega|^2$ is a biased estimate of the mean square sinusoid power, as we will see in a later section. This is noted in Alegria, 2009.

Comparison

If there is no coherent sampling relationship between the analysis window length and the frequency of interest, then the analysis scheme could be modified in one of two ways (not recommended), with unfavorable consequences:

• Use the DFT bin nearest to the frequency of interest: In this case the zero-frequency bias does not disturb detection (still orthogonal to the complex exponential corresponding to the selected bin), but detection is not perfect. The detected amplitude fluctuates with the phase of the frequency of interest.
• Allow the frequency of interest $\omega$ to be of form other than $2\pi n/N,$ where $n$ is integer: The complex exponential of frequency $\omega$ is not orthogonal to the zero-frequency or to frequency $-\omega$ complex exponentials, both of which will pollute detection. The detected amplitude which fluctuates with the phase of frequency $-\omega$. (One exception is that if the scheme is modified by allowing $N = PK / 2$ where $K$ is integer, then the positive and negative frequency of interest complex exponentials will be orthogonal, but the zero frequency is not orthogonal to those and will still pollute detection).

We can confirm the validity of the above statements for your problem by testing detection on a zero-noise synthetic version of your signal by calculating the normalized detected amplitude as:

$$\frac{2}{N}\left|\sum_{k=0}^{N-1}\Big(b + \cos\left(2\pi \tfrac{7}{45} k + \alpha\right)\Big)e^{-i2\pi f k}\right|$$

where $f$ is the detector frequency which is a DFT bin frequency if it is in format $M/N,$ where $M$ is an integer, $\alpha$ is the phase of the unit-amplitude sine wave to be detected, and $b$ is the zero frequency bias. The results look like this:

Figure 2. Fourier method detected amplitude as function of phase $\alpha$ for:
black: $N = 225, f = 7/45, b = 0 \text{ or } b = 1,$
red: $N = 256, f = 7/45, b = 0,$
blue: $N = 256, f = 7/45, b=1,$ and
purple: $N = 256, f = 40/256, b = 0 \text{ or } b = 1.$
For the purple curve, $f$ chooses the nearest bin, noting that $256\times7/45 = 39.8222\ldots$

As can be seen in Fig. 2, perfect amplitude detection is achieved when $N = 225$ (black line), not when $N = 256$ (other colored lines).

SNR estimation

The residual noise $r$ is formed by subtracting from the noisy signal $x$ the fitted sinusoid $y$ and the mean of $x$:

$$X_0 = \frac{1}{N}\sum_{k=0}^{N-1}x[k],\quad X_\omega = \frac{1}{N}\sum_{k=0}^{N-1}x[k]e^{-i\omega k}$$ $$y[k] = 2\operatorname{Re}(X_\omega)\cos(\omega k) - 2\operatorname{Im}(X_\omega)\sin(\omega k)\quad\text{for all $k$}$$ $$r[k] = x[k] - X_0 - y[k]\quad\text{for all $k$}$$

Figure 3. Residual noise $r[k]$ in coherent sampling for the first 225 samples of channel 1.

The signal-to-noise ratio (SNR) can be estimated by:

$$\text{SNR} = \frac{\frac{1}{N}\sum_{k=0}^{N-1}y[k]^2}{\frac{1}{N}\sum_{k=0}^{N-1}r[k]^2} = \frac{\sum_{k=0}^{N-1}y[k]^2}{\sum_{k=0}^{N-1}r[k]^2}$$

For the first 256 samples of the channel 1 (Figs. 1 and 3), the estimated signal-to-noise ratio is 4746 (36.8 dB). With coherent sampling, the signal model components and the residual are orthogonal, and we can write the same SNR estimate simply as:

$$\mathrm{SNR} = \frac{2|X_\omega|^2}{\frac{1}{N}\sum_{k=0}^{N-1}x[k]^2 - X_0^2 - 2|X_\omega|^2}$$

Bias correction

With white additive noise, the above coherent sampling SNR estimate will be optimistically biased, because some of the noise power is subtracted from itself by hiding in the three DFT bins of the detected signal components. This implies that the detected amplitude is also biased; the more there is true noise, the larger the amplitude detected. The following (rather intuitively constructed) formula gives the spectral subtraction estimate of the mean square of the signal, if the noise is white:

$$2|X_\omega|^2 - \frac{2}{N}\left(\frac{N}{N-3}\left({\frac{1}{N}\sum_{k=0}^{N-1}x[k]^2 - X_0^2 - 2|X_\omega|^2}\right)\right)$$

The guiding principle of this is that the square of the absolute value of each of the $N$ DFT bins has on average an equal portion of the total noise power added to it. The mean square of noise is first estimated (the part of the formula in the outer parentheses) based on the $N-3$ bins that do not include frequencies $-\omega,$ $0,$ and $\omega,$ and the estimated contribution of the noise to the biased estimate of mean square of the signal, $2|X_\omega|^2,$ is subtracted from it. This subtraction may give a negative estimate, which would result in a larger estimation error than if the estimate was clamped to zero from below. This is done in the clamped spectral subtraction estimate of the mean square of the signal:

$$\max\left(2|X_\omega|^2 - \frac{2}{N-3}\left({\frac{1}{N}\sum_{k=0}^{N-1}x[k]^2 - X_0^2 - 2|X_\omega|^2}\right),\, 0\right)$$

Figure 4. Simulated coherent sampling ($N=225,$ $\omega = 2\pi\times7/45$) mean square signal estimation in presence of additive Gaussian white noise. From top to bottom: Mean of @Simon's method estimate (green), the biased estimate or equivalently the IEEE 1057 3-parameter estimate (red), the clamped spectral subtraction estimate (purple), the spectral subtraction estimate (yellow), and true (blue) signal power, through 1000 Monte Carlo estimations, as function of true signal power.

Figure 5. Mean square of estimation error as a fraction of mean square signal power in presence of additive Gaussian white noise. The same simulation and color code as in Fig. 4. (Less is better.)

Figure 6. The same as Fig. 5 but relative to the mean fractional estimation error for the biased estimate of mean square signal. (Less is better.)

Fig. 4 demonstrates the seemingly unbiased quality (all the way to the Monte Carlo noise floor) of the spectral subtraction estimator. The clamped spectral subtraction estimate is biased. Mean square of the white Gaussian noise process was kept at 1 (0 dB), so the horizontal axis is also the true SNR. As expected, the plots were found not to be affected by the amount of zero frequency offset. @Simon's method is described in a section further down. @Simon's method estimate of mean square signal is:

$$\frac{1}{N}\sum_{k=0}^{N-1}x[k]^2 - X_0^2$$

More important than the mean of multiple estimates (Fig. 4), which shows whether there is bias, is the mean square error of the estimates (Fig. 5). Compared (Fig. 6) to the biased estimator, the spectral subtraction estimator and the clamped spectral subtraction estimator become useful (somewhat more accurate) for SNR below about -10 dB, and @Simon's method becomes useful (does not give too much more error) for SNR above about 30 dB. At very low SNR, compared to the biased estimator, @Simon's method is about 38 dB worse, and the spectral subtraction estimator about 3 dB better, clamping improving this to 4 dB better, even though the clamped estimator is (more) biased. At high SNR, there appears to be no accuracy penalty from choosing a more advanced method of the four, only extra computational cost.

The implied clamped spectral subtraction SNR estimator derived from the clamped spectral subtraction estimator of mean square signal is:

$$\mathrm{SNR} = \frac{\max\left(2|X_\omega|^2 - \frac{2}{N-3}\left({\frac{1}{N}\sum_{k=0}^{N-1}x[k]^2 - X_0^2 - 2|X_\omega|^2}\right),\,0\right)}{\frac{N}{N-3}\left({\frac{1}{N}\sum_{k=0}^{N-1}x[k]^2 - X_0^2 - 2|X_\omega|^2}\right)} = \max\left(\frac{2|X_\omega|^2\left(1-\frac{3}{N}\right)}{{\frac{1}{N}\sum_{k=0}^{N-1}x[k]^2 - X_0^2 - 2|X_\omega|^2}}-\frac{2}{N},\,0\right)$$

When plotted, the SNR estimators looks nearly identical to the curves of the corresponding mean square signal estimators of Fig. 4.

Further ideas

If there is not much noise and quantization error is significant, and you want to get the best performance, then you could ensure that the shortest period in the discrete signal matches your analysis window length. This avoids repeating the same quantization error (See Fig. 1 for the repeats of the waveform) for equally repeated values of the complex exponential in the DFT bin calculation. These coherent repeats cause an increase in accumulation of quantization error. Just changing your sine table length from 10 to 11 should increase the shortest period from 45 samples to 99 samples $(7/45\times{10/11}$ $=$ $14/99).$ This changes the oscillator frequency, too, which may be a problem. A higher DAC sampling frequency and sine generator table size would give more flexibility in fine tuning.

If you know the phase of the sinusoid, you can improve the detection further. Instead of a cosine and a sine table, which you "compare" the signal against in DFT, you could have a single table already in the correct phase.

It may be that you have both phase noise and additive noise in your data.

@Simon's method

The algorithm you are currently using can also be analyzed in zero-noise conditions, in which it calculates the detected amplitude from the noise-free synthetic signal as:

$$\sqrt{2}\times\sqrt{\frac{1}{N}\sum_{k=0}^{N-1}\Big(b + \cos\left(2\pi \tfrac{7}{45} k + \alpha\right)\Big)^2 - \left(\frac{1}{N}\sum_{k=0}^{N-1}\Big(b + \cos\left(2\pi \tfrac{7}{45} k + \alpha\right)\Big)\right)^2},$$

Under ideal conditions, the method works perfectly with $N = 225$ and has better performance at $N = 256$ than any of the Fourier approaches for the same $N:$

Figure 7. @Simon's zero-frequency power subtraction method detected amplitude as function of phase $\alpha$ for:
black: $N = 225, b = 0 \text{ or } b = 1$ and
red: $N = 256, b = 0 \text{ or } b = 1.$

Your method subtracts from the estimated total power the estimated zero-frequency power. It does not have the wide-band frequency discrimination quality of the suggested Fourier approach (See: Fourier transform of the rectangular window) and will have inferior performance if there is significant wide-band noise in the signal (see Figs. 4 and 5). One wide-band noise source is analog-to-digital converter (ADC) quantization noise. On average, its power stays roughly constant over different amplitudes of the signal of interest.

Octave script for Figs. 4 and 5

Set M = 1000 (very slow) in the script to reproduce Figs. 4-6. After the main calculation, you can run the plotting parts of the script separately to get individual plots.

pkg load statistics;
graphics_toolkit("gnuplot");

N = 225;  # Input data length
omega = 2*pi*7/45;  # Signal and detection frequency, must be DFT bin freq.
M = 100;  # Monte Carlo sample size (plot was made with 1000, which is slow)
A_signal_min = sqrt(2)*10^(-60/20);  # Minimum sinusoid amplitude, must be positive
A_signal_max = sqrt(2)*10^(60/20);  # Maximum sinusoid amplitude, must be positive
S = 200;  # Number of graph points is S+1
DC = 1;  # DC offset
MS_noise = 1;  # Noise process mean square level (plot labels correct if 1)
plot_IEEE1057 = 0;  #1: plot IEEE1057 3-parameter estimate equivalent to the algorithm in Four-Parameter Sinefit version 2.0.0.0 by Marko Neitola, 0: do not plot
plot_SNR = 0;  #1: plot SNR estimates, 0: do not plot
A_signal = A_signal_min*(A_signal_max/A_signal_min).^((0:S)/S);
MS_signal = A_signal.^2/2;
k = 0:N-1;
phasor = exp(-i*omega*k);
if plot_IEEE1057
  D0=[cos(omega*k)' sin(omega*k)' ones(N,1)];
  [Q,R] = qr(D0,0);
end

mean_estimate_MS_biased = zeros(1, S+1);
mean_estimate_MS_spectral_subtraction = zeros(1, S+1);
mean_estimate_MS_clamped_spectral_subtraction = zeros(1, S+1);
mean_estimate_MS_Simon = zeros(1, S+1);
mean_estimate_MS_IEEE1057 = zeros(1, S+1);
mean_estimate_SNR = zeros(1, S+1);
mean_estimate_SNR_spectral_subtraction = zeros(1, S+1);
mean_estimate_SNR_clamped_spectral_subtraction = zeros(1, S+1);
MS_frac_err_estimate_MS_biased = zeros(1, S+1);
MS_frac_err_estimate_MS_spectral_subtraction = zeros(1, S+1);
MS_frac_err_estimate_MS_clamped_spectral_subtraction = zeros(1, S+1);
MS_frac_err_estimate_MS_Simon = zeros(1, S+1);
MS_frac_err_estimate_MS_IEEE1057 = zeros(1, S+1);
MS_frac_err_estimate_SNR = zeros(1, S+1);
MS_frac_err_estimate_SNR_spectral_subtraction = zeros(1, S+1);
MS_frac_err_estimate_SNR_clamped_spectral_subtraction = zeros(1, S+1);
for m = 1:M
  for s = 1:S+1
    x = random("normal", 0, MS_noise, [1, N]) + A_signal(s)*cos(omega*k + 2*pi*rand()) + DC;
    MS_x = mean(x.*x);
    X_0 = mean(x);
    X_omega = mean(x.*phasor);
    estimate_MS_biased = 2*abs(X_omega)^2;
    mean_estimate_MS_biased(s) += estimate_MS_biased;
    MS_frac_err_estimate_MS_biased(s) += ((estimate_MS_biased-MS_signal(s))/MS_signal(s))^2;
    estimate_MS_spectral_subtraction = 2*abs(X_omega)^2 - 2/(N-3)*(MS_x - X_0^2 - 2*abs(X_omega)^2);
    mean_estimate_MS_spectral_subtraction(s) += estimate_MS_spectral_subtraction;
    MS_frac_err_estimate_MS_spectral_subtraction(s) += ((estimate_MS_spectral_subtraction-MS_signal(s))/MS_signal(s))^2;
    estimate_MS_clamped_spectral_subtraction = estimate_MS_spectral_subtraction;
    if estimate_MS_clamped_spectral_subtraction < 0 estimate_MS_clamped_spectral_subtraction = 0; end
    mean_estimate_MS_clamped_spectral_subtraction(s) += estimate_MS_clamped_spectral_subtraction;
    MS_frac_err_estimate_MS_clamped_spectral_subtraction(s) += ((estimate_MS_clamped_spectral_subtraction-MS_signal(s))/MS_signal(s))^2;
    estimate_MS_Simon = MS_x - X_0^2;
    mean_estimate_MS_Simon(s) += estimate_MS_Simon;
    MS_frac_err_estimate_MS_Simon(s) += ((estimate_MS_Simon-MS_signal(s))/MS_signal(s))^2;
    if plot_IEEE1057 
      x0 = R\(Q.'*x');
      estimate_MS_IEEE1057 = (x0(1)*x0(1) + x0(2)*x0(2))/2;
      mean_estimate_MS_IEEE1057(s) += estimate_MS_IEEE1057;
      MS_frac_err_estimate_MS_IEEE1057(s) += ((estimate_MS_IEEE1057-MS_signal(s))/MS_signal(s))^2;
    end
    if plot_SNR
      estimate_SNR = 2*abs(X_omega)^2/(MS_x - X_0^2 - 2*abs(X_omega)^2);
      mean_estimate_SNR(s) += estimate_SNR;
      estimate_SNR_spectral_subtraction = 2*abs(X_omega)^2*(1-3/N)/(MS_x - X_0^2 - 2*abs(X_omega)^2) - 2/N;
      mean_estimate_SNR_spectral_subtraction(s) += estimate_SNR_spectral_subtraction;
      estimate_SNR_clamped_spectral_subtraction = max(2*abs(X_omega)^2*(1-3/N)/(MS_x - X_0^2 - 2*abs(X_omega)^2) - 2/N, 0);
      mean_estimate_SNR_clamped_spectral_subtraction(s) += estimate_SNR_clamped_spectral_subtraction;
    end
  endfor
endfor
mean_estimate_MS_biased /= M;
mean_estimate_MS_spectral_subtraction /= M;
mean_estimate_MS_clamped_spectral_subtraction /= M;
mean_estimate_MS_Simon /= M;
if plot_IEEE1057 mean_estimate_MS_IEEE1057 /= M; end
if plot_SNR
  mean_estimate_SNR /= M;
  mean_estimate_SNR_spectral_subtraction /= M;
  mean_estimate_SNR_clamped_spectral_subtraction /= M;
end
MS_frac_err_estimate_MS_biased /= M;
MS_frac_err_estimate_MS_spectral_subtraction /= M;
MS_frac_err_estimate_MS_clamped_spectral_subtraction /= M;
MS_frac_err_estimate_MS_Simon /= M;
if plot_IEEE1057 MS_frac_err_estimate_MS_IEEE1057 /= M; end

# First plot
hold off
plot(10*log10(MS_signal/MS_noise), 10*log10(MS_signal));
hold on
plot(10*log10(MS_signal/MS_noise), 10*log10(mean_estimate_MS_biased));
plot(10*log10(MS_signal/MS_noise), 10*log10(mean_estimate_MS_spectral_subtraction));
plot(10*log10(MS_signal/MS_noise), 10*log10(mean_estimate_MS_clamped_spectral_subtraction));
plot(10*log10(MS_signal/MS_noise), 10*log10(mean_estimate_MS_Simon));
if plot_IEEE1057 plot(10*log10(MS_signal/MS_noise), 10*log10(mean_estimate_MS_IEEE1057), "+"); end
hold off
xlabel("true signal power and SNR (dB)");
ylabel("mean estimated signal power (dB)");
xlim([10*log10(A_signal_min^2/2), 10*log10(A_signal_max^2/2)])
ylim([10*log10(A_signal_min^2/2), 10*log10(A_signal_max^2/2)])

# Second plot
hold off
plot([10*log10(A_signal_min^2/2), 10*log10(A_signal_min^2/2)], [10*log10(A_signal_min^2/2), 10*log10(A_signal_min^2/2)]);  # Dummy, to get colors right
hold on
plot(10*log10(MS_signal/MS_noise), 10*log10(MS_frac_err_estimate_MS_biased));
plot(10*log10(MS_signal/MS_noise), 10*log10(MS_frac_err_estimate_MS_spectral_subtraction));
plot(10*log10(MS_signal/MS_noise), 10*log10(MS_frac_err_estimate_MS_clamped_spectral_subtraction));
plot(10*log10(MS_signal/MS_noise), 10*log10(MS_frac_err_estimate_MS_Simon));
if plot_IEEE1057 plot(10*log10(MS_signal/MS_noise), 10*log10(MS_frac_err_estimate_MS_IEEE1057), "+"); end
hold off
xlabel("true signal power and SNR (dB)");
ylabel("mean square fractional estimation error (dB)");
xlim([10*log10(A_signal_min^2/2), 10*log10(A_signal_max^2/2)])
grid on

# Third plot
hold off
plot([10*log10(A_signal_min^2/2), 10*log10(A_signal_min^2/2)], [10*log10(A_signal_min^2/2), 10*log10(A_signal_min^2/2)]);  # Dummy, to get colors right
hold on
plot(10*log10(MS_signal/MS_noise), 10*log10(MS_frac_err_estimate_MS_biased./MS_frac_err_estimate_MS_biased));
plot(10*log10(MS_signal/MS_noise), 10*log10(MS_frac_err_estimate_MS_spectral_subtraction./MS_frac_err_estimate_MS_biased));
plot(10*log10(MS_signal/MS_noise), 10*log10(MS_frac_err_estimate_MS_clamped_spectral_subtraction./MS_frac_err_estimate_MS_biased));
plot(10*log10(MS_signal/MS_noise), 10*log10(MS_frac_err_estimate_MS_Simon./MS_frac_err_estimate_MS_biased));
if plot_IEEE1057 plot(10*log10(MS_signal/MS_noise), 10*log10(MS_frac_err_estimate_MS_IEEE1057./MS_frac_err_estimate_MS_biased), "+"); end
hold off
xlabel("true signal power and SNR (dB)");
ylabel("mean square fractional estimation error improvement (dB)");
xlim([10*log10(A_signal_min^2/2), 10*log10(A_signal_max^2/2)])
ylim([-10, 40]);
grid on

if plot_SNR
  # Fourth plot
  hold off
  plot(10*log10(MS_signal/MS_noise), 10*log10(MS_signal));
  hold on
  plot(10*log10(MS_signal/MS_noise), 10*log10(mean_estimate_SNR));
  plot(10*log10(MS_signal/MS_noise), 10*log10(mean_estimate_SNR_spectral_subtraction));
  plot(10*log10(MS_signal/MS_noise), 10*log10(mean_estimate_SNR_clamped_spectral_subtraction));
  hold off
  xlabel("true signal power and SNR (dB)");
  ylabel("mean estimated SNR (dB)");
  xlim([10*log10(A_signal_min^2/2), 10*log10(A_signal_max^2/2)])
  ylim([10*log10(A_signal_min^2/2), 10*log10(A_signal_max^2/2)])

  # Fifth plot
  hold off
  plot(10*log10(MS_signal/MS_noise), 10*log10(MS_signal));
  hold on
  plot(10*log10(MS_signal/MS_noise), 10*log10(mean_estimate_SNR));
  plot(10*log10(MS_signal/MS_noise), 10*log10(mean_estimate_SNR_spectral_subtraction));
  plot(10*log10(MS_signal/MS_noise), 10*log10(mean_estimate_SNR_clamped_spectral_subtraction));
  hold off
  xlabel("true signal power and SNR (dB)");
  ylabel("mean estimated SNR (dB)");
  xlim([10*log10(A_signal_min^2/2), 10*log10(A_signal_max^2/2)])
  ylim([10*log10(A_signal_min^2/2), 10*log10(A_signal_max^2/2)])  
end

Summary

In summary, if possible, choose a window length that fits an integer number of periods of your frequency of interest, and directly calculate the magnitude of the DFT bin at the frequency of interest. This gives no phase-dependent fluctuation and is unaffected by zero frequency bias, see Fig. 2. Otherwise this approach to detection will give error, even in zero-noise zero-drift conditions, even when using the correct frequency for detection instead of a DFT bin frequency. If your SNR may be at times poor, test if the clamped spectral subtraction estimator improves your results.

References

Francisco Alegria, "Bias of amplitude estimation using three-parameter sine fitting in the presence of additive noise", June 2009, Measurement 42(5):748-756 (doi, pdf)

Answer 2

A quick and dirty might be to simply find all the local maxima in a zone, sort them, and average the top few. Do the same for the minima in reverse. Subtract the latter from the former, divide by two.

That is not guaranteed to be precise though and you have to ensure that you don't have a whole number of sample points per cycle or they will all line up the same each cycle.

The proper way is to have two basis vectors and find your signal as a linear combination of them. Yes, you calculate a single DFT bin at the frequency of interest. This will give you amplitude, and phase if you want it. The good news is you can get good results from a single cycle, but take two or more for noise mitigation. In the sample program I chose five.

Here is the program. You will need to tweak it a bit to make it recognize your zones, file length, etc. A frame that spans an amplitude change boundary will obviously get a bad read.

You can make the process a lot more efficient by knowing where the zones are and taking fewer frame samples. Lot of ways to take this.

The idea is to give you some working code you can play with so you can learn the math.

import numpy as np


#=============================================================================
def main():

#---- Read in the Data File

        f = open( "data.csv", "r" )

        theRows = []

        for row in f:
            theRows.append( row )

        f.close()

        theRowCount = len( theRows ) 

#---- Build Signal Array

        theChannelValues  = np.zeros( theRowCount )
        theChannelIndex   = 0

        for r in range( theRowCount ) :
            theValues = theRows[r].rstrip().split( "," )
            theChannelValues[r] = theValues[theChannelIndex]

#---- Estimate Frequency

        thePeakCount = 0
        theFirstSpot = -1
        theLastSpot = -1

        for s in range( 10, 210 ) :
            if theChannelValues[s] >= theChannelValues[s-1] and \
               theChannelValues[s] >= theChannelValues[s+1] :
                print "Max at " + str( s )
                if theFirstSpot == -1 :
                   theFirstSpot = s
                else:
                   theLastSpot = s
                   thePeakCount += 1

        print theFirstSpot, theLastSpot                   

        theCyclesPerSample = float( thePeakCount ) \
                           / float(theLastSpot - theFirstSpot )

        print "Samples per cycle: " + str( 1.0 / theCyclesPerSample )

#---- Build Basis Vectors

        theCyclesPerFrame = 5

        theSamplesPerFrame = theCyclesPerFrame / theCyclesPerSample

        theRadiansPerSample = theCyclesPerSample * 2.0 * np.pi

        N = int( theSamplesPerFrame + 0.5 )

        theC = np.zeros( N )
        theS = np.zeros( N )

        for n in range( N ) :
            theAngle = theRadiansPerSample * n
            theC[n]  = np.cos( theAngle )
            theS[n]  = np.sin( theAngle )

#---- Measure Signal Frame by Frame

        theFactor = 2.0 / float( N)

        theStart = 0

        for f in range( 100 ) :
            theLimit = theStart + N

            theDotC = theChannelValues[theStart:theLimit].dot( theC )
            theDotS = theChannelValues[theStart:theLimit].dot( theS )

            theMag = theFactor * np.sqrt( theDotC*theDotC + theDotS*theDotS ) 

            print f, theMag

            theStart = theLimit


#=============================================================================
main()

Here is an output excerpt for those who don't care to or can't run it.

Samples per cycle: 6.43333333333
0 542.814867988
1 548.585610118
2 547.425033167
3 551.011654022
4 545.654363432
5 550.587903304
6 553.360391478
7 537.68546946
8 391.858803202
9 395.333637678
10 392.753814378
11 391.284856065
12 389.907641197
13 385.425898115
14 382.921216041
15 359.681851586
16 22.9742734918
17 22.4568332805
18 22.4377296186
19 21.8988868527
20 22.1896466722
21 21.955041401
22 22.1747325702
23 21.7120535962
24 244.361827764
25 243.919290421

---

This code for basis building will be slightly more efficient:

        
        theAngle = 0

        for n in range( N ) :
            theC[n]   = np.cos( theAngle )
            theS[n]   = np.sin( theAngle )
            theAngle += theRadiansPerSample

Also, it occurred to me that you might not know how to do a dot product.

        
        theDotC = theChannelValues[theStart:theLimit].dot( theC )
        theDotS = theChannelValues[theStart:theLimit].dot( theS )

These two lines can be replaced by this loop:

        
        theDotC = 0.0
        theDotS = 0.0

        for n in range( N ) :
            v = theChannelValues[theStart+n]
            theDotC += v * theC[n]
            theDotS += v * theS[n]

---

Followup to Simon's comment:

Interesting, Olli's answer is sort of an elaboration of mine. Noticed he keyed in on two points. One, the commensurability (ratios of small integers) leads to repeat patterns. Two, you want your "DFT frame" to be as close to a whole number of cycles as you can make it. This concentrates your fundamental information of interest into one bin. The DC and harmonics are orthogonal to the basis vectors so they get zeroed out in your results automatically.

If you didn't know the frequency exactly, you could set your length to a whole number plus a half cycles, then calculate the bin on either side. You can find formulas in my blog for an exact frequency calculation (not an approximation!) and exact magnitude and phase calculation independent of frame alignment. I don't think that is necessary unless somehow you need super precision. As long as you are within one sample of a whole number of cycles, your answer is going to be darn accurate. Lengthening the frame reduces any small error introduced by this. So does having more samples in each cycle.

Yes, I know you knew your zone widths, but I was trying to provide a more general solution that readily adapts to your situation. Instead of doing frame after frame, my recommendation is you center your frame in the middle of a zone, pick a number of cycles that fills most of it, but leaves a little room on the edges, and apply the approach in this answer, then onto the next zone. I don't think it can be done any better to get the best accuracy with the least effort. The phase variance between frames will be irrelevant because you only care about the one result.

Numerically, the Goertzel is an equivalent bin calculation. It builds the basis vectors on the fly with a clever trick. However, it is more efficient to precalculate and store the basis vectors if you need them repeatedly. Compare the dot product loop with the Goertzel calculations. The advantage of Goertzel is that it doesn't require the memory storage for the basis vector arrays, and it calculates the sin and cos values using the principle that when you multiply two complex values you are adding their angles.

It's a nice diagram. I'm afraid I am not really a hardware guy so I wouldn't know the particulars, like delays and stuff. So just because your know the zone widths and when the boundaries for excitation are, doesn't necessarily mean you know where the boundaries are from your return stream. I would suspect it is a constant small delay though, perhaps even negligible.

I do like Olli's diagram. It does look like a technical drawing of some sort of drive shaft though. What is common in our answers is an emphasis on selectin a frame with a whole number of cycles. He did a nice job of explaining the penalty if you don't. You should also realize since you are calculating only one bin value there is no advantage whatsoever of picking an N that is a power of two or any other particular value.

I do think that a one bin problem solution like this gives an excellent foundation for understanding the DFT as a parallel set of these. I.e. if you had two tones it would be advantageous to set your frame length so both tones have a whole number of cycles in them. That way they each get their own bin without interference from the other. Unfortunately, in this case my simple peak counting frequency estimation won't work. If you did not know the frequencies ahead of time, or they were expected to be inconsistent, you could find them by doing a DFT of an arbitrary length and then applying frequency formulas at the humps.

One final point, the DFT does not care about the units or any specifics on your situation, especially the sampling rate. It cares about how many samples there are. It works with frequencies in terms of cycles per frame. It is when you converting to real world frequencies that knowing the sampling rate is required.

Answer 3

Not an answer, just storing a representative piece (first 2048 samples from the first channel) of the question's data for easy access and so that it won't be lost. Comma separated values.

2456,1952,1526,1564,2029,2514,2584,2184,1666,1479,1794,2336,2616,2396,1883,1508,1616,2111,2552,2543,2109,1624,1507,1888,2407,2612,2324,1800,1488,1680,2195,2592,2502,2025,1564,1525,1955,2472,2612,2272,1734,1474,1731,2273,2620,2467,1960,1532,1568,2036,2524,2593,2192,1672,1482,1799,2336,2620,2400,1876,1502,1611,2111,2560,2552,2111,1617,1500,1876,2408,2623,2342,1808,1486,1667,2192,2595,2508,2034,1569,1528,1951,2465,2608,2261,1727,1478,1726,2272,2610,2464,1960,1531,1564,2032,2528,2590,2185,1666,1478,1795,2340,2623,2392,1869,1492,1616,2115,2564,2552,2107,1608,1496,1873,2410,2628,2344,1807,1476,1657,2186,2592,2520,2043,1572,1524,1952,2460,2601,2262,1742,1489,1744,2263,2600,2448,1953,1544,1579,2044,2521,2581,2175,1672,1481,1802,2345,2623,2404,1872,1493,1605,2112,2576,2555,2111,1616,1496,1870,2397,2618,2336,1809,1490,1671,2192,2584,2499,2032,1563,1526,1953,2466,2608,2255,1725,1471,1728,2272,2615,2456,1952,1524,1560,2032,2515,2592,2192,1672,1484,1797,2336,2621,2400,1876,1495,1606,2108,2556,2552,2107,1616,1497,1872,2404,2622,2329,1802,1482,1668,2198,2592,2512,2035,1568,1528,1955,2464,2608,2272,1742,1481,1732,2272,2611,2464,1959,1536,1573,2040,2524,2587,2188,1676,1488,1808,2348,2621,2400,1874,1504,1607,2110,2562,2560,2112,1616,1499,1875,2409,2626,2344,1809,1488,2268,2611,1674,1872,2219,2419,2299,1970,1706,1747,2059,2363,2387,2122,1783,1680,1907,2259,2426,2259,1914,1683,1779,2115,2395,2363,2059,1747,1699,1967,2306,2419,2211,1859,1674,1819,2171,2419,2339,2015,1707,1707,2003,2347,2426,2175,1819,1667,1858,2224,2427,2303,1971,1706,1747,2063,2363,2387,2115,1786,1686,1907,2259,2419,2259,1915,1683,1779,2115,2402,2367,2067,1747,1699,1963,2303,2423,2211,1859,1674,1819,2163,2411,2331,2011,1715,1711,2011,2331,2403,2162,1815,1667,1859,2211,2432,2299,1955,1691,1739,2067,2378,2387,2115,1779,1682,1914,2263,2419,2259,1914,1683,1779,2115,2395,2371,2067,1754,1699,1963,2307,2426,2215,1859,1673,1819,2178,2418,2331,2011,1715,1715,2015,2339,2403,2163,1819,1682,1871,2219,2419,2291,1951,1691,1747,2075,2371,2387,2111,1771,1679,1920,2274,2431,2258,1907,1675,1778,2128,2403,2371,2067,1747,1695,1955,2299,2419,2219,1867,1671,1818,2170,2411,2335,2011,1715,1715,2011,2339,2411,2163,1819,1675,1866,2215,2418,2291,1955,1691,1746,2059,2363,2387,2123,1791,1683,1907,2251,2411,2258,1911,1688,1786,2115,2394,2355,2059,1747,1703,1984,2314,2419,2211,1859,1667,1815,2170,2419,2339,2018,1715,1707,2003,2339,2419,2186,1827,1679,1866,2218,2419,2295,1962,1707,1755,2074,2371,2379,2107,1775,1683,1919,2275,2427,2258,1906,1671,1771,2122,2403,2379,2075,1747,1691,2335,2403,2162,1819,2050,2050,2050,2056,2049,2050,2050,2050,2052,2055,2050,2047,2051,2056,2051,2052,2052,2052,2052,2056,2049,2052,2051,2056,2054,2050,2050,2047,2051,2056,2051,2050,2050,2052,2056,2049,2050,2050,2050,2050,2054,2050,2047,2051,2064,2051,2050,2051,2051,2050,2056,2049,2051,2051,2056,2055,2050,2050,2047,2050,2056,2051,2052,2051,2050,2056,2049,2051,2052,2051,2052,2054,2050,2047,2052,2064,2051,2050,2050,2051,2052,2056,2049,2050,2050,2052,2055,2050,2050,2047,2050,2056,2051,2051,2050,2052,2056,2049,2050,2052,2050,2052,2055,2050,2047,2051,2064,2051,2050,2050,2050,2052,2056,2049,2052,2050,2050,2054,2050,2051,2052,2050,2064,2051,2052,2052,2052,2056,2049,2050,2051,2051,2052,2055,2050,2047,2050,2056,2051,2051,2050,2050,2050,2056,2053,2052,2051,2052,2054,2050,2052,2047,2050,2064,2051,2049,2050,2052,2056,2049,2050,2050,2050,2050,2054,2050,2047,2050,2052,2051,2052,2051,2051,2050,2056,2049,2052,2051,2056,2054,2050,2050,2047,2051,2064,2051,2050,2050,2052,2056,2049,2051,2052,2050,2052,2054,2050,2052,2051,2056,2051,2050,2052,2052,2052,2056,2049,2050,2051,2052,2054,2050,2050,2047,2051,2056,2051,2050,2050,2050,2056,2053,2052,2052,2050,2050,2054,2050,2047,2052,2064,2051,2050,2050,2052,2052,2056,2049,2050,2050,2056,2054,2052,2051,2047,2050,2056,2051,2052,2051,2050,2056,2049,2051,2064,2051,2160,1920,1788,1896,2137,2304,2248,2020,1816,1815,2024,2256,2302,2128,1887,1789,1922,2176,2320,2223,1989,1804,1836,2064,2278,2288,2090,1861,1791,1951,2195,2309,2186,1945,1790,1855,2098,2292,2272,2056,1838,1800,1990,2232,2320,2163,1916,1784,1884,2132,2312,2245,2018,1816,1816,2022,2250,2298,2124,1887,1788,1915,2168,2312,2224,1983,1801,1840,2066,2282,2290,2087,1855,1792,1955,2205,2315,2192,1948,1792,1858,2104,2293,2272,2056,1840,1799,1982,2225,2306,2156,1916,1787,1888,2132,2302,2239,2013,1816,1816,2024,2272,2295,2124,1887,1787,1919,2171,2312,2219,1983,1802,1838,2065,2275,2288,2096,1855,1790,1952,2200,2313,2192,1947,1788,1855,2098,2300,2273,2058,1840,1803,1983,2231,2312,2160,1920,1788,1891,2132,2304,2244,2017,1820,1813,2021,2251,2298,2119,1884,1784,1916,2176,2320,2219,1982,1800,1840,2062,2277,2288,2091,1861,1794,1958,2202,2316,2194,1951,1795,1864,2104,2296,2272,2060,1837,1802,1987,2232,2311,2156,1915,1785,1888,2136,2307,2246,2019,1818,1820,2025,2256,2304,2124,1888,1783,1918,2169,2312,2224,1983,1802,1836,2064,2276,2287,2089,1855,1792,1952,2199,2320,2194,1951,1790,1864,2099,2292,2272,2057,1839,1801,1983,2226,2308,2160,1919,1788,1890,2136,2303,2239,2015,1818,1819,2032,2255,2301,2125,1887,1788,1919,2168,2312,2219,1984,1808,1843,2067,2280,2288,2092,1861,1796,2232,2315,2176,1880,1732,1871,2177,2373,2288,2000,1760,1775,2040,2320,2364,2140,1843,1724,1904,2216,2385,2268,1965,1740,1792,2080,2343,2352,2104,1816,1733,1944,2251,2384,2224,1926,1736,1823,2124,2360,2323,2064,1791,1752,1990,2288,2384,2187,1886,1727,1863,2168,2380,2293,2011,1762,1776,2038,2314,2364,2144,1846,1727,1899,2209,2384,2260,1967,1745,1799,2082,2344,2352,2103,1818,1741,1945,2256,2387,2224,1926,1733,1828,2126,2361,2322,2056,1786,1751,1982,2283,2376,2184,1885,1727,1861,2168,2376,2300,2009,1758,1768,2032,2320,2367,2144,1846,1728,1904,2211,2380,2258,1966,1744,1800,2081,2340,2345,2100,1815,1736,1940,2248,2384,2224,1923,1727,1823,2128,2364,2325,2056,1786,1750,1983,2287,2380,2192,1896,1744,1867,2166,2367,2292,2016,1774,1777,2036,2313,2360,2135,1844,1727,1905,2214,2384,2259,1962,1741,1796,2088,2345,2352,2100,1812,1736,1943,2252,2384,2217,1918,1727,1828,2128,2360,2320,2047,1777,1749,1990,2288,2375,2184,1884,1726,1862,2174,2379,2299,2014,1764,1772,2033,2314,2366,2144,1848,1727,1896,2208,2384,2264,1963,1742,1794,2080,2340,2350,2097,1814,1736,1944,2247,2374,2212,1910,1713,1816,2107,2352,2314,2049,1784,1745,1983,2284,2377,2192,1894,1732,1866,2176,2384,2299,2014,1768,1776,2035,2318,2365,2144,1848,1727,1895,2200,2367,2244,1951,1727,1779,2059,2320,2336,2088,1793,1716,2288,2367,1869,439,363,1807,3489,3892,2634,947,176,1062,2848,3936,3327,1616,306,479,2066,3647,3808,2382,741,191,1289,3104,3959,3120,1380,234,652,2335,3775,3680,2120,576,258,1541,3298,3946,2882,1148,175,836,2596,3880,3519,1867,428,355,1803,3488,3900,2633,942,170,1056,2854,3936,3332,1612,309,482,2067,3647,3808,2384,751,201,1294,3088,3968,3120,1390,243,658,2340,3775,3683,2122,584,266,1545,3308,3945,2879,1152,179,838,2596,3876,3519,1864,434,359,1805,3488,3904,2638,945,177,1059,2848,3952,3335,1612,310,484,2072,3643,3800,2376,738,187,1288,3081,3962,3120,1377,231,648,2335,3782,3680,2128,579,262,1544,3300,3950,2877,1153,184,841,2608,3879,3518,1864,427,351,1803,3488,3892,2634,941,170,1055,2848,3937,3344,1620,317,489,2070,3648,3804,2379,748,197,1292,3086,3961,3120,1383,241,656,2343,3774,3680,2124,585,267,1545,3298,3945,2884,1160,189,844,2600,3888,3519,1862,429,356,1804,3487,3891,2635,947,176,1064,2849,3938,3332,1616,316,489,2070,3647,3808,2384,739,187,1288,3081,3960,3118,1373,228,647,2336,3775,3680,2120,578,260,1544,3299,3946,2879,1148,175,837,2595,3877,3519,1855,428,354,1804,3488,3900,2635,941,170,1055,2848,3948,3333,1617,318,490,2071,3646,3808,2376,747,197,1291,3082,3968,3120,1384,241,655,3888,3515,2428,1560,1104,1486,2373,2979,2751,1932,1215,1223,1968,2790,2944,2320,1467,1118,1599,2492,2992,2638,1793,1166,1311,2128,2879,2892,2169,1357,1133,1727,2619,3003,2536,1672,1119,1385,2246,2936,2825,2046,1278,1168,1849,2712,2992,2419,1560,1110,1494,2384,2984,2737,1915,1205,1227,1983,2806,2952,2298,1452,1112,1603,2495,3000,2656,1800,1157,1300,2116,2876,2912,2175,1352,1126,1721,2616,3003,2536,1674,1119,1392,2252,2937,2828,2048,1279,1175,1854,2720,2988,2432,1564,1115,1493,2378,2976,2748,1921,1211,1229,1983,2808,2951,2300,1458,1114,1599,2491,2994,2644,1796,1158,1304,2111,2876,2896,2175,1366,1136,1727,2610,2996,2544,1675,1124,1392,2256,2940,2825,2040,1269,1164,1851,2719,2986,2420,1556,1104,1491,2378,2978,2745,1919,1212,1229,1983,2816,2960,2295,1450,1107,1599,2496,3000,2643,1792,1156,1301,2124,2877,2896,2168,1354,1129,1726,2610,2998,2533,1672,1119,1388,2250,2936,2832,2046,1277,1176,1858,2720,2983,2420,1563,1114,1497,2384,2979,2741,1922,1216,1240,1989,2806,2952,2304,1456,1111,1599,2500,3000,2648,1795,1160,1302,2118,2880,2895,2173,1357,1129,1722,2614,3000,2536,1673,1119,1392,2243,2936,2827,2046,1276,1161,1846,2714,2992,2425,1559,1100,1484,2376,2992,2747,1927,1213,1227,1984,2799,2949,2301,1456,1112,1599,2494,2994,2642,1796,1168,1307,2117,2876,2896,2172,1353,1128,2736,2995,1358,1999,2685,2807,2261,1547,1255,1680,2435,2850,2551,1835,1291,1422,2118,2758,2767,2161,1471,1286,1792,2541,2855,2470,1742,1274,1499,2227,2807,2715,2057,1399,1310,1893,2623,2854,2384,1655,1268,1594,2343,2842,2649,1963,1362,1378,2015,2697,2813,2273,1565,1278,1697,2449,2856,2555,1839,1301,1427,2123,2760,2771,2162,1471,1281,1794,2545,2866,2471,1748,1281,1509,2239,2812,2715,2065,1412,1322,1901,2621,2844,2367,1644,1263,1591,2342,2840,2640,1948,1343,1363,2007,2697,2821,2277,1561,1269,1687,2440,2855,2562,1852,1312,1429,2111,2739,2751,2161,1482,1295,1793,2534,2846,2449,1727,1270,1508,2245,2815,2711,2043,1388,1307,1898,2633,2854,2365,1635,1245,1575,2327,2840,2647,1951,1341,1352,1992,2683,2806,2264,1551,1260,1677,2432,2845,2547,1831,1292,1418,2111,2749,2763,2154,1464,1272,1780,2530,2853,2465,1736,1268,1495,2226,2804,2713,2057,1404,1315,1900,2623,2848,2367,1640,1254,1583,2333,2835,2631,1946,1342,1357,1998,2683,2809,2266,1554,1260,1678,2433,2846,2547,1835,1297,1422,2111,2751,2763,2154,1465,1276,1783,2534,2856,2461,1735,1266,1495,2227,2803,2710,2047,1399,1309,1889,2616,2840,2362,1632,1245,1572,2325,2826,2623,1938,1332,1348,1991,2681,2800,2251,1534,1240,1662,2415,2837,2539,1823,1281,1404,2097,2736,2750,2141,1449,1258,1765,2516,2842,2452,1724,1257,1484,2833,2635,1944,1342