Trouble Generating a Sine Wave in C/Matlab from Basics

So I've managed to fix it.  My problem was I wasn't properly taking into account the Sampling Rate of the produced raw wav file, and that by using the loop variable, there is a discontinuity for each iteration of the loop.

The final working code is as follows:

#include <sys/types.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>

#include <unistd.h>

#define M_PI 3.1415926535897932384626433832795L
#define NUM_CHANNELS 2
#define NUM_SAMPLES 512
#define SAMPLE_RATE 44100

double frequency = 500.0;

int main(int argc, char *argv[])
{
    if (argc >= 2){
        frequency = atof(argv[1]);
    }

    int16_t samples[NUM_CHANNELS * NUM_SAMPLES];
    unsigned cbBuffer=sizeof(samples);
    int counter = 0;

    while(1){
        for (int i = 0; i < NUM_CHANNELS * NUM_SAMPLES; i+=NUM_CHANNELS){
            samples[i] = (int16_t) 3000 * sin(frequency * 2 * M_PI * (double)counter /(double)(SAMPLE_RATE));
            samples[i+1] = samples[i];
            counter++;
        }

        int done=write(STDOUT_FILENO, samples, cbBuffer);
        if(done<0){
            fprintf(stderr, "%s : Write to stdout failed, error=%s.", argv[0], strerror(errno));
            exit(1);
        }else if(done!=cbBuffer){
            fprintf(stderr, "%s : Could not read requested number of bytes from stream.\n", argv[0]);
        }
    }

    return 0;
}


The OP is implementing the Hogenauer Filter (thank you Eugene! http://read.pudn.com/downloads163/ebook/744947/123.pdf), also called a CIC Filter, as an efficient equivalent of the moving average filter, and is getting a noise error result that is 10x more than expected.
Short Answer
The reason for the additional error in the OP's case is due to not having an extended precision accumulator.
Long Answer
We will show what the predicted noise is, for both the properly designed Moving Average and CIC filters, and then simulation results of the various structures as confirmation.
Both structures are shown below with the optional scaling for normalization, properly located at the output. The upper drawing as the Moving Average Filter is a moving average over 11 samples, and the lower drawing is mathematically equivalent as the Hogenauer or Cascade-Integrator-Comb (CIC) Filter. (For details on why these are equivalent, see CIC Cascaded Integrator-Comb spectrum)

What Is the Expected Noise?
We will first detail what noise due to numerical precision we should expect in a properly designed moving average filter. Fixed and floating point systems will be limited by the finite quantization levels given by the precision of the number. The difference between floating point and fixed point is that with fixed point the designer (or good compiler) needs to be extra careful of overflow and underflow conditions at every output (nodes) in the design, and scale the nodes accordingly such as with bit-shifting to prevent such things from happening. With floating point, this scaling happens for us automatically by the floating point processor, with overhead stored in each number. (If time to market is important, the floating point is the way to go- but if cost and power are the primary metrics, then fixed point should be strongly considered). The diagram below details the single precision floating point representation to illustrate this. The exponent of the number is equivalent to a left or right shift, scaling the number to the ranges as shown on the left side of the diagram. So even though floating point can handle an extremely large numerical range- for any given instance, the closest number we can get to that number will always be within the precision set by the mantissa. As the exponent increases, the range of the numbers available for that given exponent increases, but we will still only have the precision of the mantissa and sign bit for the quantity of numbers we can choose from.

Single precision floating point has 25 bits of precision as given by the 23 bit mantissa, plus the sign bit, plus the Robert B-J "hidden-1" bit. Double precision floating point equivalently has 54 bits of precision.
For the implementation of CIC filters, a key item that fixed point brings is the modulo arithmetic during overflow such that the subtraction in the comb is still numerically accurate regardless of the overflow. This is not how floating point arithmetic occurs where instead an overflow (or underflow) error will ultimately occur in the expected random walk behavior of the accumulator output. (So further implementation details would be required to reset the accumulator as the accumulator output exceeds certain levels toward overflow and underflow conditions). Further the quantization error in a floating point CIC implementation will be proportional to the level at the accumulator output (which is visible in the more bursty patterns of noise in the plots below for the floating point implementations).
Related is this post on the dynamic range of floating point systems: More Simultaneous Dynamic Range with Fixed Point or Floating Point? and this excellent presentation @RBJ has made at the 2008 AES Conference https://www.aes.org/events/125/tutorials/session.cfm?code=T19 which I am not sure is available anywhere online (Robert can comment). At that other post RBJ educated me about the additional hidden bit in the dynamic range result that I had confirmed with the results in my answer there.
Quantization Noise in an Accumulator
Regardless of fixed or floating point, the noise due to the accumulation that is present in both structures (Moving Average Filter and CIC Filter) is specific to any accumulator so worth-while providing full details of that operation.
For the case of the Moving Average Filter where the accumulation is done over a fixed number of iterations, the resulting noise due to precision is stationary, ergodic, band-limited and will approach a Gaussian distribution.
In contrast, for the output of the accumulator in the CIC Filter (not the final output but the internal node) is a non-stationary non-ergodic random walk random process with otherwise similar qualities as what we will detail below for accumulator noise.
Noise due to quantization is reasonably approximated as a white noise process with a uniform distribution. The variance of a uniform distribution is $r^2/12$, where $r$ is the range; thus resulting in the $q^2/12$ variance for quantization noise with $q$ being a quantization level. What occurs as this noise is accumulated is demonstrated in the diagram below, where for any addition, the distribution at the output of the adder would be the convolution of the distributions for the noise samples being summed. For example, after one accumulation the uniform distribution at the input would convolve with the uniform distribution of the previous sample resulting in a triangular distribution also with a well known variance of $q^2/6$. We see through successive convolutions after each iteration of the accumulator that the variance grows according to:
$$\sigma_N^2 = \frac{Nq^2}{12}$$
Which is the predicted variance both at the output just prior to the scaling of the Moving Average Filter where $N$ is fixed (11 in the OP's example), and at the output of the accumulator ("Integrator") in the CIC filter, where N is a counter that increases with every sample of operation. Consistent with the Central Limit Theorem, the distribution after a fixed number of counts $N$ quickly approaches a Gaussian, and due to the obvious dependence between samples introduced in the operation will no longer be white (and given the structures themselves are low-pass filters).   The scaling by dividing by $N$, appropriately placed at the Moving Average Filter output, returns the variance to be $\sigma = q^2/12$, thus having the same variance as the input but now with a band-limited nearly Gaussian distribution.  Here we see the critically of allowing filters to grow the signal (extended precision accumulators), and if we must scale, reserve scaling for the output of the filter. Never scale by scaling the input, or scaling the coefficients! Scaling in these alternate approaches will result in an increased noise at the output.

Thus we see that the predicted noise variance due to precision at the output of the Moving Average Filter is $q^2/12$, and is a Gaussian, band-limited, ergodic and stationary noise process.
Noise at Output of CIC Filter
The noise at the output of the accumulator in the CIC implementation has a variance that increases with every sample, so is a non-stationary, non-ergodic random walk process. It is itself a low-pass filter structure, creating dependence between the samples so that they are no longer independent. We would almost at this point declare it to be unusable but then in the following differencing structure we see where the magic happens: similar to using the 2-sample Variance to measure random systems with divergent properties, the output of the delay and subtraction as done in the "Comb" is a stationary, ergodic, nearly Gaussian random process!
Specifically given the difference of the two random walk signals, namely the signal and the same random walk signal as it as $N$ samples prior, we see that the result of this difference would be the same as we achieved for the Moving Average Filter output; specifically, before scaling:
$$\sigma_N^2 = \frac{Nq^2}{12}$$
And with the final scaling operation results in the same $q^2/12$ result for the CIC Filter as was obtained for the Moving Average Filter, with all the same properties regarding stationarity, ergodicity and band-limiting.
Also to be noted here is that the accumulator output noise, as a random walk noise process, grows in variance without bound at rate $N$; this means that inevitably the accumulator output will over/under flow due to error alone. For a fixed point system this is of no consequence as long as the operation rolls over on such an overflow or underflow condition; the subsequent subtraction, as long as only one such over/under flow occurred between the signals being subtracted, would be the same result (modulo arithmetic). However in floating point, an over/under flow error will occur. We see that the very low likelihood for this to occur given the error growth rate of $N\sigma^2$ unless our signal itself is operating continuously with a minimum or maximum exponent scale. For example, with single-precision floating point, and considering a probability of occurrence bound by as large as $5\sigma$ to say "unlikely", it would take $12 \times 2^{25}/5$ which is approximately 80.5M samples for the error to traverse every exponent to then reach over/underflow. This would be a good justification to only do the CIC filter in fixed point implementations, unless it is known that both the signal magnitude and total processing duration would prohibit this condition from occurring.
Simulation Results
The first simulation is to confirm the noise characteristics and variance of the accumulator output. This was done with a uniform white noise with $q = 1$, accumulated and differenced over 11 samples following the CIC structure (no output scaling was done). The upper plot below shows the noise at the output of the accumulator as well as the delayed version of this same signal from within the comb structure prior to being differenced. We see the unbounded wandering result of this random walk signal, but we also see that due to the correlation / dependency introduced in the accumulator that the difference between these two signals is stationary and bounded as shown in the middle plot. The histogram over a longer sequences confirms the Gaussian shape, and the variance of this result, with $q=1$ in the simulation was measured to 0.907 as predicted by $Nq^2/12$ with $N = 11$. (Which is the predicted variance of the output of the CIC prior to the final divide by $11$ that is shown in the first diagram).

An FFT of the differenced signal that was in the histogram above confirms the expected band-limited result:

Finally the various implementation were compared using single precision floating point so that we could use a double precision reference model as representative of "truth" for the desired moving average computation, and allow for the ability to extend the precision appropriately in the fixed point result to confirm best practice for implementation.
For this simulation, the following models were compared with names used and descriptions below:
base: Baseline double precision moving average filter used as a reference: I compared using filter and conv with identical results, and ultimately used:
base = filter(ones(11,1),11,x);

I also confirmed that the scaling of 11 shown is effectively done at the end per the diagram.
base SP: Moving average filter same as baseline with single precision floating point, which will confirm the noise growth by a factor of $N$ due to not having an extended precision accumulator:
base_SP = y_filt_sp = filter(ones(windLen,1, "single"),single(windLen),single(x));

OP: Single Precision implementation for Hogenauer done as a for loop like the OP had done, but is significantly faster than OP's actual approach. I confirmed that the result is cycle and bit accurate to his with using a double precision variant of this. I confirmed what is shown below is functionally identical to scaling after the loop. The issue is the accumulator is not extended precision.
y_mod_sp = nan(testLen,1);
xBuff = zeros(windLen+1, 1, "single");
accum = single(0);
for a = 1:testLen
  # acccumulate
  accum += single(x(a));
  #shift into buffer
  xBuff = shift(xBuff,1);
  xBuff(1)= accum;

  # comb and scale (works same if scale moved to after loop)
  y_mod_sp(a) = (xBuff(1) - xBuff(windLen + 1)) / single(windLen);
endfor

CIC: Single Precision Floating Point CIC Implementation without extended precision accumulator:
# hogenauer with filter command
y_hog_sp = filter(single([1 0 0 0 0 0 0 0 0 0 0 -1]), single([windLen -windLen]), single(x));

CIC_ext: Single Precision Floating Point CIC with extended precision Acccumulator:
# hogenauer with filter command extended precision (demonstrating
# the benefit of scaling only at output
y_hog_sp2 = single(filter([1 0 0 0 0 0 0 0 0 0 0 -1], [windLen -windLen], x));

With the following results as presented in the plot below, showing the differences from baseline in each case (given as "base - ....").
In summary, we expect the error signal from baseline at the output of the single precision CIC filter to have a standard deviation of $\sigma = q/\sqrt{12}$ where $q = 1/2^{25}$, resulting in $\sigma = 8.6e-9$.
From the simulation, the actual results for standard deviations were (for the stationary cases):
base - OP: $\sigma = 2.1e-7$
base - CIC: (not stationary)
base - base SP: $\sigma = 2.5e-8$
base - CIC ext: $\sigma = 7.8e-9$

I do not yet understand why the precision limitation in the CIC approach using the filter command results in a random walk error and this requires further investigation. However we do see by using an extended precision accumulator as shown in the "base-CIC ext" case, the best possible performance is achieved for numerical error. Extending the precision in the OP's method would certainly result in similar performance (at a much larger run time in MATLAB but may illuminate approaches in other platforms which I suspect was the motivation for coding it in a loop).
The 'base-base SP' result demonstrates how the standard deviation will grow by $N$ if an extended precision accumulator is not used in the standard Moving Average Filter, where the result of $\sigma = 2.5e-8$ which is in close agreement to this prediction given by $\sigma = \sqrt{11/12}/2^{25} = 2.85e-8$.
The OP's result is an order of magnitude larger than expected and is quite bursty, although does appear to be stationary. The explanation for the "burstiness" of the errors for the OP model is clearer after observation of the plot of the actual signal (not the difference signal) at the accumulator output plotted below. The floating point error is proportional to this signal depending on which exponent we are in, and for each the associated error or minimum quantization level will be, for single precision floating point, $1/2^{25}$ smaller. We see from the plot of the simulation result above that the error magnitude in the output for the OP case is generally proportional to the absolute magnitude of the accumulator output, which is an unbounded random walk! It is for this reason imperative that the precision at the accumulator be extended such that the maximum deviation of the random walk result between the resulting signal and its delayed copy in the comb not exceed the final precision desired. This is the reason the OP is seeing 10x more noise in that implementation!

CODE COMPARISON IN OP'S QUESTION:
The OP's comparative code for the option using filter() should not be inside a loop! (Observe that the exact same y2 result that is itself $10^4$ samples long is simply being computed $10^4$ times.)
This would be the correct comparison below showing the Hogenauer filter (CIC) structure simulated with the filter command (y2) and compared to the OP's code for the same function (y). The filter line y2 executes the entire $10^4$ data set in 0.854 seconds on my machine, while the other code took as along as me to write this and is still crunching-- so I canceled that and reduced testLen to 3000 samples to get a quicker comparison (97.08 seconds vs 0.039 seconds):
clc
clear
windLen = 11;
testLen = 10^4;
normCoeff = 1 / windLen;
xBuff = zeros(windLen, 1);
x = randn(testLen, 1);

tic
for a = 1:testLen
    varState = 0;
    y = nan(size(x));
    xBuff(windLen + 1:windLen + length(x)) = x;
    for ind=1:length(x)
        varState = varState + xBuff(windLen + ind) - xBuff(ind);
        y(ind) = varState * normCoeff;
    end
end
toc

tic
    y2 = filter([1 0 0 0 0 0 0 0 0 0 0 -1], [11 -11], x);
toc

And the resulting error difference y-y2:

A quicker implementation in MATLAB of the Hogenauer in a loop form (in case that was really needed to be consistent with a C implementation for example) but without yet addressing the "mysterious" error contribution, would be as follows:
tic
y = nan(testLen, 1);
xBuff = zeros(windLen+1, 1);
accum = 0;
for a = 1:testLen
    # acccumulate
    accum += x(a);

    #shift into buffer
    xBuff = shift(xBuff,1);
    xBuff(1)= accum;

    # comb and scale
    y(a) = (xBuff(1) - xBuff(windLen + 1)) / windLen;
 endfor
 toc

tic
    y2 = filter([1 0 0 0 0 0 0 0 0 0 0 -1], [11 -11], x);
toc

For this case I  was able to quickly process the full $10^4$ samples resulting in comparative elapsed time of 0.038 seconds for the filter() approach (y2) vs 2.385 seconds for the loop approach (y). The difference between the two results y-y2 is plotted below:


Is the FFT magnitudes of multiple audio sources additive?

The FFT algorithm efficiently computes the complex valued DFT frequency vector $X[k]$ for a given time domain sequence $x[n]$.

The magnitude of the FFT output $X[k]$ for a given bin $m$ is the magnitude of the complex number $X[m]$, so like this :

$$X[m] = X_r[m]  + j ~ X_i[m] \implies |X[m]| = \sqrt{ X_r^2[m] + X_i^2[m] }$$

Now assuming two signals $x_1[n]$ and $x_2[n]$ with their corresponding complex valued DFT vectors $X_1[k]$ and $X_2[k]$, we can see from the magnitudes of complex numbers that

$$| X_1[k] + X_2[k] | \neq | X_1[k] | + | X_2[k] | $$

hence the magnitude of FFT/DFT is not additive...


Envelope detection in frequency domain with MATLAB

you can use envspectrum.
a good example is shown in this example.


Depends on on what kind of feature(s) you want your CNN trained.  With 2D convolution, convolution in the F dimension of a T-F spectrogram will produce something like a real cepstrum, which has proved useful in speech recognition and musical pitch detection.


FIR Filter Design - obtaining FIR coefficients from frequency response

here's some very old source in C to do an FFT.  it's nowhere near as good as FFTW, but it's free and self-contained.

if you're doing this inverse DFT and windowing (i would recommend the Kaiser window) to design an FIR, probably it doesn't need to be an FFT.  you can write a very simple DFT (no Cooley-Tukey radix-2 optimizations) and run it.  i don't think you'll care too much if it takes a couple of seconds longer to run.

/*
A set of utility programs to compute the Fast Fourier Transform (FFT):

                       N-1
                X[k] = SUM { x[n]exp(-j2(pi)nk/N) }
                       n=0

and inverse Fast Fourier Transform (iFFT):

                           N-1
                x[n] = 1/N SUM { X[k]exp(+j2(pi)nk/N) }
                           k=0

To speed things up, multiplication by 1 and j are avoided.  The input, x[],
is an array of complex numbers (pairs of doubles) of length N = 2^p.  The
calling program supplies p = log2(N) not the array length, N.  The output
is placed in BIT REVERSED order in x[].  Call bitreverse(x, p) to swap back
to normal order, if needed.  However, iFFT(X, p, trigtbl) requires its input,
X[], to be in bit reversed order making bit reversing unnecessary in some
cases, such as fast convolution.  trigtbl is an array of doubles of length
N/4 containing the sin function from 0 to pi/2 used to compute the FFT.
Call sintable(trigtbl, p) ONLY ONCE before either FFT(x, p, trigtbl) or
iFFT(X, p, trigtbl) to set up a sin table for FFT computation.

Written in Megamax C (for the Mac) by Robert Bristow-Johnson (1985).
*/

/*
#include <complex.h>
*/

typedef struct {
double real;
double imag;
} complex;

#define Re(z) (z).real
#define Im(z) (z).imag

#define PI (3.14159265359)

double sin();

FFT (x, p, trigtbl)
complex x[];
int p;
register double *trigtbl;
{
register long length, step, stepsize, size;
register complex *top, *bottom, *end;               /* top & bottom of FFT butterfly */
complex temp;

size = 1L<<p;
end = x + size;

length = size>>1;
size >>= 2;
stepsize = 1L;
while ( length > 1L) {
    top = x;
    while (top < end) {
        bottom = top + length;

        Re(temp) = Re(*top) - Re(*bottom);          /* butterfly: twiddle = 1 */
        Im(temp) = Im(*top) - Im(*bottom);
        Re(*top) += Re(*bottom);
        Im(*top) += Im(*bottom);
        *bottom = temp;
        top++;
        bottom++;

        for (step = stepsize; step < size; step += stepsize) {
            Re(temp) = Re(*top) - Re(*bottom);      /* butterfly: twiddle = exp(-jÌ) */
            Im(temp) = Im(*top) - Im(*bottom);
            Re(*top) += Re(*bottom);
            Im(*top) += Im(*bottom);
            Re(*bottom) = Re(temp)*trigtbl[size - step] + Im(temp)*trigtbl[step];
            Im(*bottom) = Im(temp)*trigtbl[size - step] - Re(temp)*trigtbl[step];
            top++;
            bottom++;
        }

        Re(temp) = Im(*top) - Im(*bottom);          /* butterfly: twiddle = -j */
        Im(temp) = Re(*bottom) - Re(*top);
        Re(*top) += Re(*bottom);
        Im(*top) += Im(*bottom);
        *bottom = temp;
        top++;
        bottom++;

        for (step = stepsize; step < size; step += stepsize) {
            Re(temp) = Im(*top) - Im(*bottom);      /* butterfly: twiddle = -j*exp(-jÌ) */
            Im(temp) = Re(*bottom) - Re(*top);
            Re(*top) += Re(*bottom);
            Im(*top) += Im(*bottom);
            Re(*bottom) = Re(temp)*trigtbl[size - step] + Im(temp)*trigtbl[step];
            Im(*bottom) = Im(temp)*trigtbl[size - step] - Re(temp)*trigtbl[step];
            top++;
            bottom++;
        }
        top = bottom;
    }
    stepsize <<= 1;
    length >>= 1;
}

top = x;
bottom = x + 1;
while (top <  end) {
    Re(temp) = Re(*top) - Re(*bottom);              /* butterfly: twiddle = 1 */
    Im(temp) = Im(*top) - Im(*bottom);
    Re(*top) += Re(*bottom);
    Im(*top) += Im(*bottom);
    *bottom = temp;
    top += 2;
    bottom += 2;
}
}

iFFT (X, p, trigtbl)
complex X[];
int p;
register double *trigtbl;
{
register long length, step, stepsize, size;
double scale;
register complex *top, *bottom, *end;               /* top & bottom of FFT butterfly */
complex temp;

size = 1L<<p;
end = X + size;

scale = 1.0/(double)size;
top = X;
bottom = X + 1;
while (top <  end) {
    Re(temp) = (Re(*top) - Re(*bottom))*scale;      /* butterfly: twiddle = 1/N */
    Im(temp) = (Im(*top) - Im(*bottom))*scale;
    Re(*top) = (Re(*top) + Re(*bottom))*scale;
    Im(*top) = (Im(*top) + Im(*bottom))*scale;
    *bottom = temp;
    top += 2;
    bottom += 2;
}

length = 1L;
size >>= 2;
stepsize = size;
while ( stepsize >= 1L) {
    length <<= 1;
    top = X;
    while (top < end) {
        bottom = top + length;

        temp = *bottom;                             /* butterfly: twiddle = 1 */
        Re(*bottom) = Re(*top) - Re(temp);
        Im(*bottom) = Im(*top) - Im(temp);
        Re(*top) += Re(temp);
        Im(*top) += Im(temp);
        top++;
        bottom++;

        for (step = stepsize; step < size; step += stepsize) {
                                                    /* butterfly: twiddle = exp(+jÌ) */
            Re(temp) = Re(*bottom)*trigtbl[size - step] - Im(*bottom)*trigtbl[step];
            Im(temp) = Im(*bottom)*trigtbl[size - step] + Re(*bottom)*trigtbl[step];
            Re(*bottom) = Re(*top) - Re(temp);
            Im(*bottom) = Im(*top) - Im(temp);
            Re(*top) += Re(temp);
            Im(*top) += Im(temp);
            top++;
            bottom++;
        }

        Re(temp) = -Im(*bottom);                    /* butterfly: twiddle = +j */
        Im(temp) = Re(*bottom);
        Re(*bottom) = Re(*top) - Re(temp);
        Im(*bottom) = Im(*top) - Im(temp);
        Re(*top) += Re(temp);
        Im(*top) += Im(temp);
        top++;
        bottom++;

        for (step = stepsize; step < size; step += stepsize) {
                                                    /* butterfly: twiddle = +j*exp(+jÌ) */
            Re(temp) = -Im(*bottom)*trigtbl[size - step] - Re(*bottom)*trigtbl[step];
            Im(temp) = Re(*bottom)*trigtbl[size - step] - Im(*bottom)*trigtbl[step];
            Re(*bottom) = Re(*top) - Re(temp);
            Im(*bottom) = Im(*top) - Im(temp);
            Re(*top) += Re(temp);
            Im(*top) += Im(temp);
            top++;
            bottom++;
        }
        top = bottom;
    }
    stepsize >>= 1;
}
}

sintable(trigtbl, p)
register double *trigtbl;
int p;
{
register long size, i;
double theta;

size = 1L<<(p-2);
theta = (PI/2.0)/(double)size;

for (i = 0; i < size; i++)
    trigtbl[i] = sin(theta*(double)i);
}

bitreverse (x, p)
register complex *x;
int p;
{
complex temp;
register long k, i, r, size, count;

size = (1L<<p) - 1L;
for (k = 1L; k < size; k++) {
    i = k;
    r = 0;
    for (count = size; count > 0; count >>= 1) {
        r <<= 1;
        r += (i & 0x00000001L);
        i >>= 1;
    }
    if (r > k) {
        temp = x[r];
        x[r] = x[k];
        x[k] = temp;
    }
}
}


I have solved this problem by adding certain padding to overlap area.

There are pixels belonging to image1, image2, background and overlap. A pixels in overlap are successively relabeled to either image1 or image2 depending on neighborhood.

This padding will make some space for the blurring so that sharp edges of the overlap have no change of appearing.

Another treat is photometric calibration, i.e. gain compensation and vignetting removal. This minimizes or even eliminates the differences in low frequency band.

Finally, the blurring have to avoid background pixels (usually black) from bleeding into the image. This can be done by alpha blending or using binary mask as described here.

It is also possible to reduce amount of blurring near edges and corners of the overlap. This, unfortunately, reduces performance of the blending algorithm and the results are somewhat uncontrolled (as the blurring is non-uniform).

Finally, a gradient domain blending improves the quality near edges as it does not rely on Gaussian blurs of fixed sizes.


Rapid changes in signal correspond to high frequencies - proof

One might argue that the Dirac delta $\delta(t)$ , is a signal that exhibits maximum change because all frequencies are uniformly increased, not just the high frequencies


scipy's cwt is primitive and error prone; below is via ssqueezepy.cwt:

Code:
Note that if you seek to code the wavelet yourself, you'll need to take it to the frequency domain first (preferably analytically via Fourier transform and then sampling it, rather than via FFT), then pass it like wavelet = Wavelet(my_func); cwt(x, wavelet).
import numpy as np
from ssqueezepy import cwt
from ssqueezepy.visuals import plot, imshow

#%%# Helper fn + params #####################################################
def exp_am(t, offset):
    return np.exp(-pi*((t - offset) / .1)**10)

pi = np.pi
v1, v2, v3 = 64, 128, 32

#%%# Make `x` & plot #########################################################
t = np.linspace(0, 1, 2048, 1)
x = (np.sin(2*pi * v1 * t) * exp_am(t, .2) +
     (np.sin(2*pi * v1 * t) + 2*np.cos(2*pi * v2 * t)) * exp_am(t, .5)  +
     (2*np.sin(2*pi * v2 * t) - np.cos(2*pi * v3 * t)) * exp_am(t, .8))
plot(x, title="x(t) | t=[0, ..., 1], %s samples" % len(x), show=1)

#%%# Take CWT & plot #########################################################
Wx, scales = cwt(x, 'morlet')
imshow(Wx, yticks=scales, abs=1,
       title="abs(CWT) | Morlet wavelet",
       ylabel="scales", xlabel="samples")


Improving an auto-correlation based guitar pitch detector

I realized I had forgotten to make an update to this.

I ended up following @robert bristow-johnson's recommendations in the comments.

I used time-domain McLeod Pitch Method in a live application (records audio in a loop) which works very well.

You can see the source code here:
https://github.com/sevagh/Pitcha


Variance of an Implicit Function of Kalman State Vector

I think I have the solution.
I'd be happy to hear others' thought.

Defining $ F \left(r, v, a, {T}_{tth} \right) = r + v {T}_{tth} + \frac{a {{T}_{tth}}^{2}}{2} $ which is the implicit function which connects all variables.
Since we're dealing with non linear function the variance is given by:

$$ var \left( {T}_{tth} \right) = J P {J}^{T} $$

Where $ P $ is the state vector covariance at a given time and $ J = \left[ \frac{\partial {T}_{tth} }{\partial r} \frac{\partial {T}_{tth} }{\partial v} \frac{\partial {T}_{tth} }{\partial a} \right] $

Using Total Derivative Law on the Implicit Function:
$$ \frac{\partial {T}_{tth} }{\partial r} = -\frac{{F}_{r}}{{F}_{{T}_{tth}}}, \frac{\partial {T}_{tth} }{\partial v} = -\frac{{F}_{v}}{{F}_{{T}_{tth}}}, \frac{\partial {T}_{tth} }{\partial a} = -\frac{{F}_{a}}{{F}_{{T}_{tth}}} $$

Where $ {F}_{x} $ is the partial derivative of $ F $ with respect to $ x $.

Each individual component would be:

$$
\frac{\partial {T}_{tth} }{\partial r}=-\frac{1}{v+aT_{tth}}, \frac{\partial {T}_{tth} }{\partial v}=-\frac{T_{tth}}{v+aT_{tth}}, \frac{\partial {T}_{tth} }{\partial a}=-\frac{T^2_{tth}}{2\left(v+aT_{tth}\right)}
$$

Transform the state covariance matrix with $J P {J}^{T}$ and you'd have the right solution.

Given a diagonal state covariance matrix, we should have:

$$
P_{tth}=\frac{\sigma^2_r}{\left(v+aT_{tth}\right)^2}+\frac{T^2_{tth}\sigma^2_v}{\left(v+aT_{tth}\right)^2}+\frac{T^4_{tth}\sigma^2_a}{4\left(v+aT_{tth}\right)^2}
$$


Detail in images require higher frequency basis functions.  The frequency in this case is measuring fluctuations in intensity as a distance is traversed.  With a lot of detail, a lot of fluctuations, thus higher frequency.
Tamp down the higher frequency and you lose detail, i.e. the image blurs.
The dampening is measured (described) best on a log scale.
Consider this answer Honey Mustard.


Is the sum of discrete sinusoids an eigenfunction of an LTI system?

Given an LTI system with impulse response $h(n)$, a complex sinusoid of the form $e^{j\omega_0n}$ is always an eigenfunction:

$$\mathcal{T}\{ e^{j\omega_0 n} \} = \sum_{k=-\infty}^{\infty} h(k) e^{j \omega_0 (n-k)} = e^{j \omega_0 n} \left( \sum_{k=-\infty}^{\infty} h(k) e^{-j \omega_0 k}\right) = H\left(e^{j\omega_0}\right) \cdot e^{j \omega_0 n}$$

Note, however, that the sum of eigenfunctions corresponding to different eigenvalues is not an eigenfunction. You can see this question for a demonstration in a generic vector space. In this case it can be easily seen that, because the system is LTI (in particular, linear), the sum of the respective eigenfunctions will appear at the output:

$$\mathcal{T}\{ e^{j\omega_0 n} + e^{j2\omega_0 n} \} =
H\left(e^{j\omega_0}\right) \cdot e^{j \omega_0 n} + H\left(e^{j2\omega_0}\right) \cdot e^{j 2\omega_0 n} $$

As you can see, there is no way to express the output as $K\cdot\left(e^{j\omega_0 n} + e^{j2\omega_0 n}\right)$, except in the case where the two eigenvalues (as stated above, $H\left(e^{j\omega_0}\right)$ and $H\left(e^{j2\omega_0}\right)$ in this case) are the same.


I have dealt with similar problems before. Probably the first you need to know is how does the noise spectrum look like.

First approach (mostly empirical): I see your signal of interest (SoI) is a train of rectangular pulses, thus the main frequency components are around DC. You can low-pass filter it to remove the main noise components (probably a moving-average filter will do). Then set a detection threshold (empirically), to a value that will be always surpassed by the filtered SoI, and rarely or never by the filtered noise.

Second approach (mostly theoretical): Check the noise components in the frequency domain (apply FFT). If the noise is nearly flat from DC to $f_s/2$, you can apply the optimal solution to this detection problem, which is implementing an energy detector, as follows:

Assuming we're in presence of (zero mean) Gaussian noise, estimate its statistics (noise variance will do).

How to do it?
Since your SoI is a train of rectangular pulses, its main frequency components are around DC. You can high-pass filter it to remove the SoI. If the filter's cut-off frequency is low enough, the filtered noise will keep a similar noise variance after filtering. Otherwise, multiply the after-filter variance by $(f_s/2)/(f_s/2 - W)$, where $W$ is the filter's bandwidth, and $f_s$ is the sampling rate.
Apply the optimal threshold for detecting a signal in AWGN:
$$
\lambda = 2~\sigma^2_\omega ~ \sqrt{N}~ \left(Q^{-1}(P_{fa}) + \sqrt{N}\right)~~~~[1]
$$

$\sigma^2_\omega$ is the noise variance, $N$ is the signal length (number of samples of the sliding buffer, see step 3), $P_{fa}$ is the design false-alarm probability (i.e., how likely you would admit a false alarm to be triggered).
Use a sliding buffer to analyze the whole signal by pieces, and what you need to compare to $\lambda$ is the energy contained in such a buffer: $$E = \sum^{N-1}_{n = 0} |x[n]|^2$$

More details on the second approach can be found, for example, here:

[1]  S. Atapattu, C. Tellambura, and H. Jiang, Energy Detection for Spectrum Sensing in Cognitive Radio, Springer 2014. Chapter 2:
Conventional Energy Detector.


Loss of precision due to low-pass filtering in physical units

In general, no.  In fact, if the low-pass filter characteristics were well-chosen in the sense of getting the best post-filtering signal to noise ratio, then the amount that you might degrade any such measurements would be minimal.

In part this conclusion comes from my "well-chosen" qualifier.  Roughly speaking (and I can only speak roughly because your problem is so vaguely stated), if you've got a true feature in the object that you're imaging that resolves to an image feature as sharp as your 1-pixel line then unless you have absolutely tremendous noise your $\sigma = 5$ Gaussian filter is poorly chosen.

I would assume that for looking at tissues, you're looking for features that are significantly larger than 0.05mm, and in x-rays you're getting significant white noise -- hence, a filter that smooths out the noise without smoothing out the features that you can actually reasonably expect to be there is only going to make things -- including measurement precision -- better, or at least no worse.

(But if I had a shop accident and had some AWG 36 gauge copper wire embedded in me, then I would presume that your $\sigma = 0.25\mathrm{mm}$ filter may blur the edges too much -- but that's a different problem, and thank goodness that copper wire shows up very strongly on an x-ray).

In general, what you can do is to set out the problem with more specifics: i.e., what's a typical feature you're looking at, how do its edges naturally merge into the surroundings, what's its contrast in the image, and how much noise is there in the image.  With all of these in hand, then you can start asking how much a given filter will impact tasks such as identifying the location of an edge (or, since the subject is mammography, presumably identifying the size of a lesion).


a single complex number in the FFT result table is, in fact, a representation of some (small and determined by the sampling rate and length of signal) freq band

Not really. The DFT is first and foremost a mathematical procedure that simply allows representing a time domain sequence as a sum of complex exponentials. The physical interpretation of that operation depends a lot on the circumstances but in general it's NOT true that FFT coefficient at some frequency $\omega$ represents accurately what's happening physically at frequency $\omega$. They are related but the exact nature of the relationship depends a lot on the details.
The notable exception here are time sequences that are periodic with the FFT length, but periodicity implies infinite length which doesn't exist in the real world.

want to estimate a phase (single number) for a frequency band of 9-11 Hz.

That would require you to first define what a single number phase for a frequency band is. Typically the phase is different at each frequency, so typically a single number doesn't
work well.


Check the scaling factors used in your DFT/IDFT implementation. There is an arbitrary scaling factor that differs by implementation. Specifically, if you take an input vector $x$ and round-trip it via a DFT followed by an IDFT:

$$
y = IDFT(DFT(x))
$$

Then you may see a couple different results in common implementations:

$y = x$
$y = Nx$, where $N$ is the transform length

Other cases that I wouldn't be too surprised to see:

$y = \sqrt{N}x$
$y = \frac{1}{N}x$
$y = \frac{1}{\sqrt{N}}x$

Which one of these that you see depends on which scaling factors are used in the forward and reverse transforms. If you want the norm of the DFT of $x$ and the norm of $x$ to be the same, then you want the unitary DFT definition:

$$
X[k] = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1}x[n]e^{\frac{-j2\pi k n}{N}}
$$


How to form Kalman filtering matrices for a problem with variable acceleration?

Take the first model and train your model noise matrix $\bf{Q}$ to account for the variable acceleration. A popular model is the Singer model from 1990 (see [1]).

[1] Singer, R. A., “Estimating Optimal Tracking Filter Performance for Manned Maneuvering Targets,” IEEE Transactions on Aerospace and Electronic Systems, vol. AES-6, no. 4, .pp. 473-483, July 1990.


How does BRIEF produce consistent keypoints with random sampling?

It looks to me as if they came up with a random pixel patch vs pixel patch set of points and they used this to train their system.  The system then recognizes different transformations by looking at the database of transformations and seeing how the transformations would affect the match.  It then picks the transformation with the best score (lowest Hamming distance).  The random pixel locations are only random once during the system design.


Studying spectral leakage due to DFT length and frequency resolution

It looks like what you are doing is only changing the DFT size (using N1 variable), this is equivalent to padding your signals with zeros up to a given length. Your signal (for N1=200) looks like:

It can be immediately seen that the periodic extension is not really a continuous waveform. This discontinuity is introducing the leakage effect that you are observing.

Below a quick and dirty code that demonstrates that. I did both padding of input signal in time domain, as well as calling the fft function with twice the signal size. Plots are in dB scale.

F0=10; %Frequency of the sinusoid
Fs=100; %Sampling Frequency
observationTime = 1; %observation time in seconds
t=0:1/Fs:observationTime-1/Fs; %time base

%% Ideal case
x=sin(2*pi*F0*t);%sampled sine wave

N1=100; %DFT length same as signal length
X1 = abs(1/N1*fftshift(fft(x,N1)));%N-point complex DFT of x
f1=(-N1/2:1:N1/2-1)*Fs/N1; %frequencies on x-axis
X1 = 10*log10(X1/max(X1));
subplot(2,1,1)
plot(f1,X1, 'o-');
title('Ideal case')
grid on

%% Padding signal with zeros
x2 = [x zeros(1, length(x))];

N2=200; %DFT length same as signal length
X2 = abs(1/N2*fftshift(fft(x, N2)));%N-point complex DFT of x
f2=(-N2/2:1:N2/2-1)*Fs/N2; %frequencies on x-axis
X2 = 10*log10(X2/max(X2));
subplot(2,1,2)
plot(f2,X2, 'o-');
hold on
grid on

%% Padding with zeros implicitly
N2=200; %DFT length same as signal length
X3 = abs(1/N2*fftshift(fft(x, N2)));%N-point complex DFT of x
f3=(-N2/2:1:N2/2-1)*Fs/N2; %frequencies on x-axis
X3 = 10*log10(X3/max(X3));
plot(f3,X3, 'r-');
title('Padding with zeros')
legend({'Padded vector', 'Twice the FFT size'})


Derivation of the lowpass to bandpass transformation

You're right that the multiplication of a low pass and a high pass filter results in a band pass filter, as long as the cut-off frequency of the low pass is higher than the cut-off frequency of the high pass. The problem with that approach is that low pass and high pass filters with magnitude responses that are optimal according to some chosen criterion (Butterworth, Chebyshev, Cauer) will not result in an optimum band pass filter.
On the other hand, mapping a single optimal filter will result in another optimal filter. Using $\omega_l\omega_u=\omega_0^2$, where $\omega_l$ and $\omega_u$ are lower and upper band edges, respectively, and $\omega_0$ is the center frequency of the band pass filter, and leaving out constants for the sake of simplicity, the transformation can be written as
$$s\longmapsto \frac{s^2+\omega_0^2}{s}\tag{1}$$
[Note that $\omega_l$ and $\omega_u$ are denoted as $\omega_1$ and $\omega_2$ in the OP, but they're used in a different way in the figure below.]
The mapping $(1)$ maps DC ($\omega=0$) to the desired center frequency $\omega_0$. Furthermore, $s=\pm\infty$ is mapped to $s=0$ and $s=\infty$. So the whole frequency axis of the low pass filter is mapped to the positive frequency axis of the band pass filter. (The same is true for the negative half axis of the band pass filter):

(from: Digital Filter Design by Parks and Burrus)


I think you are referring to the symbol shown at the left of the structure: the arm that moves between different phases?  See slide 4 of this PDF.

The point the Polyphase representation is trying to get across is that we give different samples of the input to different Phases of the Polyphase filter.  Polyphase Filtering is a Digital Signal Processing approach, so your input would only be defined at discrete locations.  Between the samples, we do not know where the signal was (unless it is bandlimited...).

So you can think of it as a memoryless switch.


How do wavelet levels depend on one another?

I may be a little confused by your question (intuitive dependency?) but I think these two sources will help you wrap your head around the relationship between levels in a DWT:

https://mil.ufl.edu/nechyba/www/eel6562/course_materials/t5.wavelets/intro_dwt.pdf
http://users.rowan.edu/~polikar/WTpart4.html

I think one way to think of it is because a DWT can be inverted, you can take a signal and perform a DWT on it and then you can reverse that DWT to get the original signal. So there's a clear dependance between levels in a DWT. I don't expect this answer to be accepted but i hope it can point you in the right direction to what you're looking for.

Small quote from link #1

"In wavelet analysis, the Discrete Wavelet Transform (DWT) decomposes a signal into a set of mutually orthogonal wavelet basis functions. These functions diﬀer from sinusoidal basis functions in that they are spatially localized – that is, nonzero over only part of the total signal length. Furthermore, wavelet functions are dilated, translated and scaled versions of a a common function φ, known as the mother wavelet. As is the case in Fourier analysis, the DWT is invertible, so that the original signal can be completely recovered from its DWT representation.
"


Why is image blurring called convolution and not cross correlation?

Broadly speaking, image blurring replaces each pixel value by a weighted sum of adjacent pixel values. Now, this weighted sum can be represented as convolution or as a correlation
for the simple reason that a convolution of $I$ and $g$  (where $I$ is the image and $g$ the
kernel) is the same as the correlation of $I$ and $\hat{g}$ where $\hat{g}$ is just
$g$ "flipped over". As fhucho has already pointed out, if the kernel is symmetric, the
flipping over has no effect.

It is easier to understand this idea in 1D rather than 2D.  For each $n$, suppose that
the result (weighted sum of adjacent values) that we need to get is
$$\begin{align}
\hat{I}[n] &= \quad I[n]a[0]\\
&\quad +\ I[n+1]a[1] + I[n+2]a[2] + \cdots + I[n+N]a[N]\\
&\quad +\ I[n-1]a[-1] + I[n-2]a[-2] + \cdots + I[n-N]a[-N]\\
&= \sum_{m=-N}^N I[n+m]a[m]
\end{align}$$
We readily recognize as the correlation of $\mathbf I$ and the sequence
$$\mathbf a =\left(a[-N], a[-N+1], \cdots, a[0], \cdots a[N-1], a[N]\right).$$
Now, define the sequence $\mathbf g$ as $\mathbf  a$ flipped over, that is,
$g[m] = a[-m]$ for all $m$. Note that $\mathbf a$ is just $\mathbf g$ flipped
over (and is thus the $\hat{g}$ mentioned earlier). Then we have that
$$\begin{align}
\hat{I}[n] &= \sum_{m=-N}^N I[n+m]a[m]\\
&= \sum_{m=-N}^N I[n+m]g[-m]\\
&= \sum_{-i=-N}^N I[n-i]g[i] &\text{replace}~ m~\text{by} ~ -i\\
&= \sum_{i=-N}^N I[n-i]g[i] & \text{write sum in reverse order}\\
&= \mathbf I \star \mathbf g\bigr|_n
\end{align}$$
A similar calculation can be carried out for 2D signals (images) but
the details are messier and sometimes it is hard to see the forest
because of the trees.

In summary, whether we choose to regard image blurring as
correlation or convolution
is entirely a matter of convenience and nomenclature.
We can convolve the image and the "impulse response"
or correlate the image with the "esnopser eslupmi",
whichever we prefer, and in either case, we are replacing
each pixel value with the same weighted sum of adjacent pixel values.


Correction of amplitude after zero padding for upsampling purposes

If you use an 1/N normalized DFT, then you shouldn't have any problems with your amplitude when you take the inverse DFT (with no normalization factor/factor of 1).  Consider the case of a pure tone with a whole number of cycles in the frame.  Only one bin pair will be non-zero.  No matter how many extra zeroes you insert, the inverse DFT will reconstruct the same signal in the same duration, just appropriately upsampled.  The amplitude of the signal will be twice the magnitude of of each bin value independent of your sample count.

If the normalization factors are already hard coded in your routines as "1" forward, "1/N" inverse, you can simply use them in reverse roles.  Take the "inverse" of your signal, zero pad it (splitting the Nyquist), then take the "forward".
Alternatively, just multiply your final results by a factor of $\frac{N_{new}}{N_{old}}$, which confirms your observations in the update.


How to derive the bin width / maximum frequency of DFT/STFT

The maximum frequency for a real waveform is $F_s/2$, while the maximum frequency for a complex waveform is $F_s$.  The DFT is a transformation of N samples in time to N samples in frequency, with bins numbers $0$ to $N-1$, with bin $0$ representing "DC", and bin $N-1$ representing 1 bin less than the sampling rate.

Complex signals are resented as phasors rotating in time in either a counterclockwise ($e^{j\omega t}$) or clockwise ($e^{-j\omega t}$) direction representing positive and negative frequencies specifically. Know that each bin in the DFT represents a single complex frequency ($e^{j\omega t}$), and given Euler's Identity you need two complex frequencies to represent a sinusoidal tone:

$$2\cos(\omega t) = e^{j\omega t} + e^{-j\omega t}$$

Hence the reason why you only need $F_s/2$ for real waveforms since the second tone is redundant (always the complex conjugate for real signals).

Also importantly, the DFT is cyclical as it is a sampled system, no different than how aliasing occurs when we use an A/D converter.

The simplest way to see all of this is to consider a bicycle wheel representing a complex tone, and a strobe light as our sampling clock.

The wheel can spin counter-clockwise or clock-wise, and we will refer to counter-clockwise as a negative frequency, and clockwise as a positive frequency.

If the wheel is spinning 2 cycles per second clockwise, and our strobelight is going at 10 cycles per second, we will witness the 2 cycles per second spin rate. However if our wheel is spinning 12 cycles per second clockwise, it will appear just as if it is spinning 2 cycles per second.  If our wheel was spinning 10 cylces per second, it would appear as it is not spinning at all (DC!).

This is just like the DFT, and why bin N-1 is the highest frequency. The next bin N (if we had more bins) would be the sampling frequency itself and would be no different than bin 0. This is the cyclical nature of the DFT, and how for complex signals we can map them either from $-F_s/2$ to $+F_s/2$ or as is done in the DFT from $0$ to $F_s$. (since it is cyclical either can be a unique span).

So we have N bins with each frequency as $0/F_s, 1/F_s, 2/Fs ... (N-1)/F_s$, we see that the bin width is $F_s/N$.

Ultimately and more universally a time sample of duration $T$ has a frequency resolution of $1/T$ (Time Frequency Duality). So if our sampling rate is $F_s$ samples / sec and we have $M$ samples, then the duration $T$ will be $N/F_s$ and then our frequency resolution is the inverse of this which is $F_s/M$. If we zero pad that out to $N$ total samples, the frequency spacing of each bin will decrease to $F_s/N$ but we have not done anything to improve the frequency resolution, or our ability to discern two closely spaced frequencies- it will still be $F_s/M$! So here we note the interesting point that our frequency resolution is really not to do with the sampling rate but to do with the total duration in time of our actual (non-zero padded) data set.  Here $T$ is the length of time for the actual data within the window. With zero padding we increase the number of DFT bins which therefore increases the number of frequency samples but it does not increase the frequency resolution!  See this post for more details on that with a specific example:

What happens when N increases in N-point DFT

This assumes a rectangular window, which has the narrowest frequency resolution as given by the relationship above (and specifically this is the equivalent noise bandwidth). Any other window used will decrease the frequency resolution further, as detailed in many other posts on here with regards to windowing specifically.

Here is an example demonstrating real and complex frequencies and their DFT:

N = 10;        # Total number of samples
n = 0:9;       # sample index
k = 3;         # frequency index (Fs = N)
y = cos(2*pi*k*n/N);
stem(n,fft(y))

The plot above is for the DFT of the function $y = \cos(3\omega_n)$ where $\omega_n$ is $2\pi/N$. Note that this is the same as $\frac{1}{2}e^{j3\omega_n} +\frac{1}{2}e^{-j3\omega_n}$. Because of the cyclical nature of the DFT described above, $e^{-j3\omega_n} = e^{j(10-3)\omega_n} = e^{j7\omega_n}$, so we see the two DFT bins representing this at bin 3 and bin 7.

Think of the bicylce wheel and the strobe light, if it was rotating 7 cycles per second in the positive direction with the strobe light at 10 cycles per second, it would visually appear to be rotating 3 cycles per second in the negative direction! The upper bins of the DFT represent the negative frequencies. For real signals the negative frequencies are the complex conjugate (equal magnitude, opposite phase) of the positive frequencies, hence we only need to show the frequencies from $0$ to $F_s/2$ in the case of real signals, as we already know what the negative frequencies are.

Below shows the same result for a complex tone at at bin 3 only.

y2 = exp(j*2*pi*k*n/N);
stem(n,fft(y2))

As for the statement by the OP:

  halving the window duration would remove every second frequency bin,
  but keeping the same window duration and removing every second sample
  (halving the sample rate) would correspond to removing the right half
  of the bins.

Halving the DFT duration (not the window duration but the window duration will be modified if the DFT duration is less than double the window) would remove every second frequency bin and not impact the remaining bins only under the condition the condition that the contents remaining in the window are either cyclically invariant to the original data set (cascading the remainng window with itself results in the same original data) or there is nothing in the portion of the data that was removed (removing zeros only).  Please this post which details this further, knowing what is shown in frequency at that post can equally be described as the time domain.

Effect of changing sample rate, window duration and zero padding on DTFT and DFT


Can deep neural networks achieve real-time video analysis?

Nvidia seems to have published some white papers comparing DNN inference performance between high-powered CPUs and (of course) Nvidia GPUs.  (one example)

Ballpark seems to be that some systems can meet or exceed typical video frame rate thru-put for some class of DNN classification tasks.  Whether those image sizes and/or DNN architectures and hyperparameters meets your needs is one question, and whether desktop/server-class GPUs fit your desired power envelope is another question.


How to plot noise shaped spectrum of First order Incremental Sigma Delta ADC's output?

Update: as Sundar commented his question is specific to the “incremental Sigma Delta ADC” of which I only have a vague understanding of. To the extent the incremental implementation is a traditional Sigma Delta ADC that is reset periodically I would believe this model still applies as described similar to using a moving average frequency response to model an integrate and dump. Comments welcome.

Below is how I use an equivalent loop equation to model the noise shaping in first and higher order Sigma Delta architectures. This is qualitative since the quantizer is highly non-linear (So would a apply from a small signal perspective, and can be further more accurately approximated with a piece-wise linear approximation to the quantizer gain.

Below shows continuous time equivalents, with $N_q$ as the quantization noise input from the 1 bit converter. The discrete time equivalent is easily derived by mapping the integrator. An integrator using either Backward Euler or Impulse Invariance mappings both result in:

$$\frac{1}{s} \leftrightarrow \frac{T}{1-z^{-1}}$$


Low pass and High pass filter Coefficient

To elaborate a bit on Fat32's answer: the most straightforward thing to do is to compute (or estimate) the following two sums:

$$H(e^{j 0})=\sum_nh[n]\tag{1}$$

and

$$H(e^{j \pi})=\sum_n(-1)^nh[n]\tag{2}$$

where $(1)$ is the value of the frequency response at DC (i.e., $\omega=0$), and $(2)$ is the value of the frequency response at Nyquist (i.e., at $\omega=\pi$).

A low pass filter should have a relatively large value for $(1)$ and a very small value (ideally zero) for $(2)$. For a high pass filter the opposite is the case. If both values are small (and if $h[n]$ is not zero) then it's probably a band pass filter, and if both values are relatively large, it's probably a band stop filter. This of course only applies if you can assume that the filter approximates some standard frequency selective filter characteristic.


It will match for circular convolution modulo $N$, where $N$ is 4 here. For finite length sequences product of DFT of 2 sequences is equivalent to DFT of circular-convolution of the 2 sequences.

>> cconv([0,1,2,1], [0,1,-1,1], 4)

ans =

 0     1     2     1

>>ifft(fft([0,1,2,1]).*fft([0,1,-1,1]))

ans =

 0     1     2     1


How is time used in this autocorrelation expression?

The autocorrelation function of a random process $x(t)$ is defined as $$R(t,\tau) = {\mathbb E}\{x(t) x(t+\tau)\}.$$
For a stationary random process, this function does not depend on $t$, i.e., we have $R(t,\tau) = R(\tau) \; \forall t$. The correlation only depends on the time difference between two samples, not the "absolute" time $t$ when the correlation is measured.
For a stationary and ergodic process, ensemble averaging can be replaced by time averaging and we have $$R(\tau) = \lim_{T\rightarrow \infty} \int_T x(t) x(t+\tau) {\rm d}t,$$ which in the discrete case collapses to the sum you were showing.
The expression mentioned by Brown is a function of $t$ and $\tau$, no doubt about that. However, Brown argues (I haven't followed through the arguments, but they should be in the second paper he was citing) that the expression can be used to compute a "narrowed autocorrelation function" (what ever that is) and that if viewed as a function of $\tau$ it has peaks when $\tau$ is equal to multiples of the period $T$ of the periodic signal $f(t)$. This sounds like the particular choice of $t$ is irrelevant, which makes sense if we consider $f(t)$ to be stationary and ergodic. Then again, the notation if not really rigorous. It is worth reading some more though because a few equations later he introduces the bracket notation $\langle \cdot \rangle$ for an average over a period $T$. That would make a lot of sense for his definition of $S(\tau)$ too. In fact, if we do open Brown et al 1989 in equation (1b) he does include the $\langle \cdot \rangle$ notation. Maybe it was forgotten in the other paper?


The function is based on Matlab's deconv, so reading that page should help understand it.

Here's a docstring I wrote for SciPy's deconvolve, but haven't submitted yet because I'm not sure it's 100% correct: https://github.com/scipy/scipy/pull/430#issuecomment-13675004

The input to deconvolve is signal and divisor, and your output is quotient and remainder, where signal was originally produced by signal = convolve(divisor, quotient) + remainder.

original = [0, 1, 0, 0, 1, 1, 0, 0]
impulse_response = [2, 1]
recorded = scipy.signal.convolve(impulse_response, original)
print recorded
# [0 2 1 0 2 3 1 0 0]
recovered, remainder = scipy.signal.deconvolve(recorded, impulse_response)
print recovered, remainder
# [ 0.  1.  0.  0.  1.  1.  0.  0.] [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]

I don't understand what you're trying to do with PMFs, though, so I can't give a more specific example.


How many samples per cycle(time period) is actually enough to reproduce the sine wave sound properly?

It depends on how good your interpolation between samples is.  If you have really good interpolation, anything more than 2 samples per cycle will suffice.  This is not just a theoretical fact, it is my experience in practice using polyphase interpolation that combines 64 adjacent table-lookup samples.

If you're using linear interpolation (which combines 2 adjacent table-lookup samples in the interpolation), I would use no less than 1024 samples per cycle for a high-quality sine wave (about 120 dB S/N).


Auto Detection of Rotation Angle on Arbitrary Image with Orthogonal Features

If I understand your method 1 correctly, with it, if you used a circularly symmetrical region and did the rotation about the center of the region, you would eliminate the region's dependency on the rotation angle and get a more fair comparison by the merit function between different rotation angles. I will suggest a method that is essentially equivalent to that, but uses the full image and does not require repeated image rotation, and will include low-pass filtering to remove pixel grid anisotropy and for denoising.

Gradient of isotropically low-pass filtered image

First, let's calculate a local gradient vector at each pixel for the green color channel in the full size sample image.

I derived horizontal and vertical differentiation kernels by differentiating the continuous-space impulse response of an ideal low-pass filter  with a flat circular frequency response that removes the effect of the choice of image axes by ensuring that there is no different level of detail diagonally compared to horizontally or vertically, by sampling the resulting function, and by applying a rotated cosine window:

$$\begin{gather}h_x[x, y] = \begin{cases}0&\text{if }x = y = 0,\\-\displaystyle\frac{\omega_c^2\,x\,J_2\left(\omega_c\sqrt{x^2 + y^2}\right)}{2 \pi\,(x^2 + y^2)}&\text{otherwise,}\end{cases}\\
h_y[x, y] = \begin{cases}0&\text{if }x = y = 0,\\-\displaystyle\frac{\omega_c^2\,y\,J_2\left(\omega_c\sqrt{x^2 + y^2}\right)}{2 \pi\,(x^2 + y^2)}&\text{otherwise,}\end{cases}\end{gather}\tag{1}$$

where $J_2$ is a 2nd order Bessel function of the first kind, and $\omega_c$ is the cut-off frequency in radians. Python source (does not have the minus signs of Eq. 1):

import matplotlib.pyplot as plt
import scipy
import scipy.special
import numpy as np

def rotatedCosineWindow(N):  # N = horizontal size of the targeted kernel, also its vertical size, must be odd.
  return np.fromfunction(lambda y, x: np.maximum(np.cos(np.pi/2*np.sqrt(((x - (N - 1)/2)/((N - 1)/2 + 1))**2 + ((y - (N - 1)/2)/((N - 1)/2 + 1))**2)), 0), [N, N])

def circularLowpassKernelX(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.fromfunction(lambda y, x: omega_c**2*(x - (N - 1)/2)*scipy.special.jv(2, omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  kernel[(N - 1)//2, (N - 1)//2] = 0
  return kernel

def circularLowpassKernelY(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.fromfunction(lambda y, x: omega_c**2*(y - (N - 1)/2)*scipy.special.jv(2, omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  kernel[(N - 1)//2, (N - 1)//2] = 0
  return kernel

N = 41  # Horizontal size of the kernel, also its vertical size. Must be odd.
window = rotatedCosineWindow(N)

# Optional window function plot
#plt.imshow(window, vmin=-np.max(window), vmax=np.max(window), cmap='bwr')
#plt.colorbar()
#plt.show()

omega_c = np.pi/4  # Cutoff frequency in radians <= pi
kernelX = circularLowpassKernelX(omega_c, N)*window
kernelY = circularLowpassKernelY(omega_c, N)*window

# Optional kernel plot
#plt.imshow(kernelX, vmin=-np.max(kernelX), vmax=np.max(kernelX), cmap='bwr')
#plt.colorbar()
#plt.show()

Figure 1. 2-d rotated cosine window.

Figure 2. Windowed horizontal isotropic-low-pass differentiation kernels, for different cut-off frequency $\omega_c$ settings. Top: omega_c = np.pi, middle: omega_c = np.pi/4, bottom: omega_c = np.pi/16. The minus sign of Eq. 1 was left out. Vertical kernels look the same but have been rotated 90 degrees. A weighted sum of the horizontal and the vertical kernels, with weights $\cos(\phi)$ and $\sin(\phi)$, respectively, gives an analysis kernel of the same type for gradient angle $\phi$.

Differentiation of the impulse response does not affect the bandwidth, as can be seen by its 2-d fast Fourier transform (FFT), in Python:

# Optional FFT plot
absF = np.abs(np.fft.fftshift(np.fft.fft2(circularLowpassKernelX(np.pi, N)*window)))
plt.imshow(absF, vmin=0, vmax=np.max(absF), cmap='Greys', extent=[-np.pi, np.pi, -np.pi, np.pi])
plt.colorbar()
plt.show()

Figure 3. Magnitude of the 2-d FFT of $h_x$. In frequency domain, differentiation appears as multiplication of the flat circular pass band by $\omega_x$, and by a 90 degree phase shift which is not visible in the magnitude.

To do the convolution for the green channel and to collect a 2-d gradient vector histogram, for visual inspection, in Python:

import scipy.ndimage

img = plt.imread('sample.tif').astype(float)
X = scipy.ndimage.convolve(img[:,:,1], kernelX)[(N - 1)//2:-(N - 1)//2, (N - 1)//2:-(N - 1)//2]  # Green channel only
Y = scipy.ndimage.convolve(img[:,:,1], kernelY)[(N - 1)//2:-(N - 1)//2, (N - 1)//2:-(N - 1)//2]  # ...

# Optional 2-d histogram
#hist2d, xEdges, yEdges = np.histogram2d(X.flatten(), Y.flatten(), bins=199)
#plt.imshow(hist2d**(1/2.2), vmin=0, cmap='Greys')
#plt.show()
#plt.imsave('hist2d.png', plt.cm.Greys(plt.Normalize(vmin=0, vmax=hist2d.max()**(1/2.2))(hist2d**(1/2.2))))  # To save the histogram image
#plt.imsave('histkey.png', plt.cm.Greys(np.repeat([(np.arange(200)/199)**(1/2.2)], 16, 0)))

This also crops the data, discarding (N - 1)//2 pixels from each edge that were contaminated by the rectangular image boundary, before the histogram analysis.

$\pi$
$\frac{\pi}{2}$
$\frac{\pi}{4}$
$\frac{\pi}{8}$
$\frac{\pi}{16}$
$\frac{\pi}{32}$
$\frac{\pi}{64}$
—$0$Figure 4. 2-d histograms of gradient vectors, for different low-pass filter cutoff frequency $\omega_c$ settings. In order: first with N=41: omega_c = np.pi, omega_c = np.pi/2, omega_c = np.pi/4 (same as in the Python listing), omega_c = np.pi/8, omega_c = np.pi/16, then: N=81: omega_c = np.pi/32, N=161: omega_c = np.pi/64 . Denoising by low-pass filtering sharpens the circuit trace edge gradient orientations in the histogram.

Vector length weighted circular mean direction

There is the Yamartino method of finding the "average" wind direction from multiple wind vector samples in one pass through the samples. It is based on the mean of circular quantities, which is calculated as the shift of a cosine that is a sum of cosines each shifted by a circular quantity of period $2\pi$. We can use a vector length weighted version of the same method, but first we need to bunch together all the directions that are equal modulo $\pi/2$. We can do this by multiplying the angle of each gradient vector $[X_k,Y_k]$ by 4, using a complex number representation:

$$Z_k = \frac{(X_k + Y_k i)^4}{\sqrt{X_k^2 + Y_k^2}^3} = \frac{X_k^4 - 6X_k^2Y_k^2 + Y_k^4 + (4X_k^3Y_k - 4X_kY_k^3)i}{\sqrt{X_k^2 + Y_k^2}^3},\tag{2}$$

satisfying $|Z_k| = \sqrt{X_k^2 + Y_k^2}$ and by later interpreting that the phases of $Z_k$ from $-\pi$ to $\pi$ represent angles from $-\pi/4$ to $\pi/4$, by dividing the calculated circular mean phase by 4:

$$\phi = \frac{1}{4}\operatorname{atan2}\left(\sum_k\operatorname{Im}(Z_k), \sum_k\operatorname{Re}(Z_k)\right)\tag{3}$$

where $\phi$ is the estimated image orientation.

The quality of the estimate can be assessed by doing another pass through the data and by calculating the mean weighted square circular distance, $\text{MSCD}$, between phases of the complex numbers $Z_k$ and the estimated circular mean phase $4\phi$, with $|Z_k|$ as the weight:

$$\begin{gather}\text{MSCD} = \frac{\sum_k|Z_k|\bigg(1 - \cos\Big(4\phi - \operatorname{atan2}\big(\operatorname{Im}(Z_k), \operatorname{Re}(Z_k)\big)\Big)\bigg)}{\sum_k|Z_k|}\\
= \frac{\sum_k\frac{|Z_k|}{2}\left(\left(\cos(4\phi) - \frac{\operatorname{Re}(Z_k)}{|Z_k|}\right)^2 + \left(\sin(4\phi) - \frac{\operatorname{Im}(Z_k)}{|Z_k|}\right)^2\right)}{\sum_k|Z_k|}\\
= \frac{\sum_k\big(|Z_k| - \operatorname{Re}(Z_k)\cos(4\phi) - \operatorname{Im}(Z_k)\sin(4\phi)\big)}{\sum_k|Z_k|},\end{gather}\tag{4}$$

which was minimized by $\phi$ calculated per Eq. 3. In Python:

absZ = np.sqrt(X**2 + Y**2)
reZ = (X**4 - 6*X**2*Y**2 + Y**4)/absZ**3
imZ = (4*X**3*Y - 4*X*Y**3)/absZ**3
phi = np.arctan2(np.sum(imZ), np.sum(reZ))/4

sumWeighted = np.sum(absZ - reZ*np.cos(4*phi) - imZ*np.sin(4*phi))
sumAbsZ = np.sum(absZ)
mscd = sumWeighted/sumAbsZ

print("rotate", -phi*180/np.pi, "deg, RMSCD =", np.arccos(1 - mscd)/4*180/np.pi, "deg equivalent (weight = length)")

Based on my mpmath experiments (not shown), I think we won't run out of numerical precison even for very large images. For different filter settings (annotated) the outputs are, as reported between -45 and 45 degrees:

rotate 32.29809399495655 deg, RMSCD = 17.057059965741338 deg equivalent (omega_c = np.pi)
rotate 32.07672617150525 deg, RMSCD = 16.699056648843566 deg equivalent (omega_c = np.pi/2)
rotate 32.13115293914797 deg, RMSCD = 15.217534399922902 deg equivalent (omega_c = np.pi/4, same as in the Python listing)
rotate 32.18444156018288 deg, RMSCD = 14.239347706786056 deg equivalent (omega_c = np.pi/8)
rotate 32.23705383489169 deg, RMSCD = 13.63694582160468 deg equivalent (omega_c = np.pi/16)

Strong low-pass filtering appear useful, reducing the root mean square circular distance (RMSCD) equivalent angle calculated as $\operatorname{acos}(1 - \text{MSCD})$. Without the 2-d rotated cosine window, some of the results would be off by a degree or so (not shown), which means that it is important to do proper windowing of the analysis filters. The RMSCD equivalent angle is not directly an estimate of the error in the angle estimate, which should be much less.

Alternative square-length weight function

Let's try square of the vector length as an alternative weight function, by:

$$Z_k = \frac{(X_k + Y_k i)^4}{\sqrt{X_k^2 + Y_k^2}^2} = \frac{X_k^4 - 6X_k^2Y_k^2 + Y_k^4 + (4X_k^3Y_k - 4X_kY_k^3)i}{X_k^2 + Y_k^2},\tag{5}$$

In Python:

absZ_alt = X**2 + Y**2
reZ_alt = (X**4 - 6*X**2*Y**2 + Y**4)/absZ_alt
imZ_alt = (4*X**3*Y - 4*X*Y**3)/absZ_alt
phi_alt = np.arctan2(np.sum(imZ_alt), np.sum(reZ_alt))/4

sumWeighted_alt = np.sum(absZ_alt - reZ_alt*np.cos(4*phi_alt) - imZ_alt*np.sin(4*phi_alt))
sumAbsZ_alt = np.sum(absZ_alt)
mscd_alt = sumWeighted_alt/sumAbsZ_alt

print("rotate", -phi_alt*180/np.pi, "deg, RMSCD =", np.arccos(1 - mscd_alt)/4*180/np.pi, "deg equivalent (weight = length^2)")

The square length weight reduces the RMSCD equivalent angle by about a degree:

rotate 32.264713568426764 deg, RMSCD = 16.06582418749094 deg equivalent (weight = length^2, omega_c = np.pi, N = 41)
rotate 32.03693157762725 deg, RMSCD = 15.839593856962486 deg equivalent (weight = length^2, omega_c = np.pi/2, N = 41)
rotate 32.11471435914187 deg, RMSCD = 14.315371970649874 deg equivalent (weight = length^2, omega_c = np.pi/4, N = 41)
rotate 32.16968341455537 deg, RMSCD = 13.624896827482049 deg equivalent (weight = length^2, omega_c = np.pi/8, N = 41)
rotate 32.22062839958777 deg, RMSCD = 12.495324176281466 deg equivalent (weight = length^2, omega_c = np.pi/16, N = 41)
rotate 32.22385477783647 deg, RMSCD = 13.629915935941973 deg equivalent (weight = length^2, omega_c = np.pi/32, N = 81)
rotate 32.284350817263906 deg, RMSCD = 12.308297934977746 deg equivalent (weight = length^2, omega_c = np.pi/64, N = 161)

This seems a slightly better weight function. I added also cutoffs $\omega_c = \pi/32$ and $\omega_c = \pi/64$. They use larger N resulting in a different cropping of the image and not strictly comparable MSCD values.

1-d histogram

The benefit of the square-length weight function is more apparent with a 1-d weighted histogram of $Z_k$ phases. Python script:

# Optional histogram
hist_plain, bin_edges = np.histogram(np.arctan2(imZ, reZ), weights=np.ones(absZ.shape)/absZ.size, bins=900)
hist, bin_edges = np.histogram(np.arctan2(imZ, reZ), weights=absZ/np.sum(absZ), bins=900)
hist_alt, bin_edges = np.histogram(np.arctan2(imZ_alt, reZ_alt), weights=absZ_alt/np.sum(absZ_alt), bins=900)
plt.plot((bin_edges[:-1]+(bin_edges[1]-bin_edges[0]))*45/np.pi, hist_plain, "black")
plt.plot((bin_edges[:-1]+(bin_edges[1]-bin_edges[0]))*45/np.pi, hist, "red")
plt.plot((bin_edges[:-1]+(bin_edges[1]-bin_edges[0]))*45/np.pi, hist_alt, "blue")
plt.xlabel("angle (degrees)")
plt.show()

Figure 5. Linearly interpolated weighted histogram of gradient vector angles, wrapped to $-\pi/4\ldots\pi/4$ and weighted by (in order from bottom to top at the peak): no weighting (black), gradient vector length (red), square of gradient vector length (blue). The bin width is 0.1 degrees. Filter cutoff was omega_c = np.pi/4, same as in the Python listing. The bottom figure is zoomed at the peaks.

Steerable filter math

We have seen that the approach works, but it would be good to have a better mathematical understanding. The $x$ and $y$ differentiation filter impulse responses given by Eq. 1 can be understood as the basis functions for forming the impulse response of a steerable differentiation filter that is sampled from a rotation of the right side of the equation for $h_x[x, y]$ (Eq. 1). This is more easily seen by converting Eq. 1 to polar coordinates:

$$\begin{align}h_x(r, \theta) = h_x[r\cos(\theta), r\sin(\theta)] &= \begin{cases}0&\text{if }r = 0,\\-\displaystyle\frac{\omega_c^2\,r\cos(\theta)\,J_2\left(\omega_c r\right)}{2 \pi\,r^2}&\text{otherwise}\end{cases}\\
&= \cos(\theta)f(r),\\
h_y(r, \theta) = h_y[r\cos(\theta), r\sin(\theta)] &= \begin{cases}0&\text{if }r = 0,\\-\displaystyle\frac{\omega_c^2\,r\sin(\theta)\,J_2\left(\omega_c r\right)}{2 \pi\,r^2}&\text{otherwise}\end{cases}\\
&= \sin(\theta)f(r),\\
f(r) &= \begin{cases}0&\text{if }r = 0,\\-\displaystyle\frac{\omega_c^2\,r\,J_2\left(\omega_c r\right)}{2 \pi\,r^2}&\text{otherwise,}\end{cases}\end{align}\tag{6}$$

where both the horizontal and the vertical differentiation filter impulse responses have the same radial factor function $f(r)$. Any rotated version $h(r, \theta, \phi)$ of $h_x(r, \theta)$ by steering angle $\phi$ is obtained by:

$$h(r, \theta, \phi) = h_x(r, \theta - \phi) = \cos(\theta - \phi)f(r)\tag{7}$$

The idea was that the steered kernel $h(r, \theta, \phi)$ can be constructed as a weighted sum of $h_x(r, \theta)$ and $h_x(r, \theta)$, with $\cos(\phi)$ and $\sin(\phi)$ as the weights, and that is indeed the case:

$$\cos(\phi) h_x(r, \theta) + \sin(\phi) h_y(r, \theta) = \cos(\phi) \cos(\theta) f(r) + \sin(\phi) \sin(\theta) f(r) = \cos(\theta - \phi) f(r) = h(r, \theta, \phi).\tag{8}$$

We will arrive at an equivalent conclusion if we think of the isotropically low-pass filtered signal as the input signal and construct a partial derivative operator with respect to the first of rotated coordinates $x_\phi$, $y_\phi$ rotated by angle $\phi$ from coordinates $x$, $y$. (Derivation can be considered a linear-time-invariant system.) We have:

$$\begin{gather}x = \cos(\phi)x_\phi - \sin(\phi)y_\phi,\\
y = \sin(\phi)x_\phi + \cos(\phi)y_\phi\end{gather}\tag{9}$$

Using the chain rule for partial derivatives, the partial derivative operator with respect to $x_\phi$ can be expressed as a cosine and sine weighted sum of partial derivatives with respect to $x$ and $y$:

$$\begin{gather}\frac{\partial}{\partial x_\phi} = \frac{\partial x}{\partial x_\phi}\frac{\partial}{\partial x} + \frac{\partial y}{\partial x_\phi}\frac{\partial}{\partial y} = \frac{\partial \big(\cos(\phi)x_\phi - \sin(\phi)y_\phi\big)}{\partial x_\phi}\frac{\partial}{\partial x} + \frac{\partial \big(\sin(\phi)x_\phi + \cos(\phi)y_\phi\big)}{\partial x_\phi}\frac{\partial}{\partial y} = \cos(\phi)\frac{\partial}{\partial x} + \sin(\phi)\frac{\partial}{\partial y}\end{gather}\tag{10}$$

A question that remains to be explored is how a suitably weighted circular mean of gradient vector angles is related to the angle $\phi$ of in some way the "most activated" steered differentiation filter.

Possible improvements

To possibly improve results further, the gradient can be calculated also for the red and blue color channels, to be included as additional data in the "average" calculation.

I have in mind possible extensions of this method:

1) Use a larger set of analysis filter kernels and detect edges rather than detecting gradients. This needs to be carefully crafted so that edges in all directions are treated equally, that is, an edge detector for any angle should be obtainable by a weighted sum of orthogonal kernels. A set of suitable kernels can (I think) be obtained by applying the differential operators of Eq. 11, Fig. 6 (see also my Mathematics Stack Exchange post) on the continuous-space impulse response of a circularly symmetric low-pass filter.

$$\begin{gather}\lim_{h\to 0}\frac{\sum_{N=0}^{4N + 1} (-1)^n f\bigg(x + h\cos\left(\frac{2\pi n}{4N + 2}\right), y + h\sin\left(\frac{2\pi n}{4N + 2}\right)\bigg)}{h^{2N + 1}},\\
\lim_{h\to 0}\frac{\sum_{N=0}^{4N + 1} (-1)^n f\bigg(x + h\sin\left(\frac{2\pi n}{4N + 2}\right), y + h\cos\left(\frac{2\pi n}{4N + 2}\right)\bigg)}{h^{2N + 1}}\end{gather}\tag{11}$$

Figure 6. Dirac delta relative locations in differential operators for construction of higher-order edge detectors.

2) The calculation of a (weighted) mean of circular quantities can be understood as summing of cosines of the same frequency shifted by samples of the quantity (and scaled by the weight), and finding the peak of the resulting function. If similarly shifted and scaled harmonics of the shifted cosine, with carefully chosen relative amplitudes, are added to the mix, forming a sharper smoothing kernel, then multiple peaks may appear in the total sum and the peak with the largest value can be reported. With a suitable mixture of harmonics, that would give a kind of local average that largely ignores outliers away from the main peak of the distribution.

Alternative approaches

It would also be possible to convolve the image by angle $\phi$ and angle $\phi + \pi/2$ rotated "long edge" kernels, and to calculate the mean square of the pixels of the two convolved images. The angle $\phi$ that maximizes the mean square would be reported. This approach might give a good final refinement for the image orientation finding, because it is risky to search the complete angle $\phi$ space at large steps.

Another approach is non-local methods, like cross-correlating distant similar regions, applicable if you know that there are long horizontal or vertical traces, or features that repeat many times horizontally or vertically.


Simplify equation of single pole IIR transfer function

$$H(z)\cdot H(z^{-1})= \frac{K}{1-\alpha \cdot e^{-j \omega}} \cdot \frac{K}{1-\alpha \cdot e^{+j \omega}}$$

$$= \frac{K^2}{1-\alpha \cdot (e^{-j \omega}+e^{+j \omega})+ \alpha^2} = \frac{K^2}{1 + \alpha^2- 2 \cdot \alpha \cdot Re\{e^{j \omega}\}} = \frac{K^2}{1 + \alpha^2- 2 \cdot \alpha \cdot \cos( \omega)} $$


Solving for impulse response h[n] given input-output pairs

Without any generalizations, this particular problem could be solved with the following: Let the causal FIR impulse response of the LTI system be $h[n]$ of length $L$, which can be inferred from the first set that $L = 4$, hence let the impulse response be $$ h[n] = \{a,b,c,d\} $$

Then from the first set, it can be seen that
$$y[n] = h[n]+h[n-2] = \{a,b,c+a,d+b,c,d\} = \{0,1,0,2,0,1\}$$

from which you can deduce that $$a=0 ~,~ b=1 ~,~ c=0 ~,~ d=1 $$ which confirms that your result was right. For a longer impulse response you would need to solve a more involved set of (linear) equations...


mathematical transfer for different camera attributes

This Process is Known as Color Correction.It is Explained in more detail here.

To Give a short Explanation:

Lets suppose

You have a Input Color Matrix represented as

    | R_O1 G_O1 B_O1 |
A = | R_O2 G_O2 B_O2 |
    | R_On G_On B_On |

and a Reference Target Color Matrix represented as

    | R_T1 G_T1 B_T1 |
B = | R_T2 G_T2 B_T2 |
    | R_Tn G_Tn B_Tn |

Color Correction Matrix X is computed as follows:

            **A X = B**

    | R_O1 G_O1 B_O1 |     | A11 A12 A13 |   =  | R_T1 G_T1 B_T1 |
    | R_O2 G_O2 B_O2 |   * | A21 A22 A23 |      | R_T2 G_T2 B_T2 |
    | R_On G_On B_On |     | A31 A32 A33 |      | R_Tn G_Tn B_Tn |

X = (B.A^T)(A.A^T)^-1

Some one please help me editing the formulas!


It is quite very difficult to give a correct answer, because you do not give the full information. But I will try.
One of standard way is using "matching filter":

http://en.wikipedia.org/wiki/Matched_filter

http://www.dspguide.com/ch17/3.htm

If length of your signal pattern is quite small (<50 usec) you can use this method.


Finding minimum in the middle of a histogram

To me, this looks like you want something like

start from right. Add up bins' values, as long as

the added bin is larger than the last or
the value added by that bin divided by the current sum is still above a fixed threshold

Adjust the threshold to fit your needs.

Often, things like smoothing the histogram (effectively: applying a low pass filter to it) help. Or you do something like finding a function fit, something like a polynomial of fixed degree (let's say 5) that minimizes some error function (for example, sum square absolute error in all bins), and then analytically find local extrema.


how can the channel fading model (Rayleigh, or Rice, or Log normal distribution) infer the way channel estimation is done?

My personal point of view: what matters much more than the distribution is the second-order statistics of the channel. In particular things like delay spread / coherence bandwidth and Doppler spread / coherence time. These parameters tell you a thing or two that you must take into account when estimating your channels:

The coherence bandwidth (roughly equal to the inverse of the delay spread) tells you for what bandwidth you can consider the channel approximately flat. If the bandwidth you want to transmit is below that, you can consider the channel a scalar quantity, which is simple to estimate. If your bandwidth is higher, you either need to go to frequency-domain techniques like OFDM and then estimate the channel on several frequencies (spaced no more than the coherence bandwidth) or estimate channel taps in time domain, e.g., when using a Rake receiver.
The coherence time (roughly equal to the inverse Doppler spread) tells you how quickly the channel changes in time, i.e., for how long you can consider it approximately constant. After the coherence time has passed, your channel estimate is no longer valid and you need to train again (or employ tracking if you can).

That said, the channel conditions (like LOS/NLOS) do influence these quantities as well. If you have a strong LOS, chances are your coherence time is higher since the LOS will not change rapidly in time.


Is it theoretically possible to perfectly quantize a continuous signal?

Depends on assumptions you are willing to make and what type of signals are you trying to sample, but in theory I think that sampling rate equal to the Planck time would be a gold standard for anything...

This translates to sampling frequency of $1.855 \times 10 ^ {43} \mathtt{Hz}$ ($18.55$ tredecillion hertz). Personally I believe that machines will never be so fast. Obviously we can limit our assumptions to more realistic ones.

Obviously there is more constraints, such as: what sensor is capable of measuring displacement/velocity with such accuracy? Obviously this assumption is applicable even for subatomic scale, and most likely we don't need such accuracy in measurements of sound signals. Additionally when frequency is increasing, sound attenuation in media (air in this case) is increasing to enormous amounts and ultra-sonic sound waves are not propagating on large distances. At some point you will start measuring the Brownian motion. You also mentioned about possibility of measuring each particle separately and it is another, very difficult problem to tackle.


Is it necessary to taper two ends before fft a time series?

Tapering the two ends of a short audio signal allows it to be played at an arbitrary time after and before silence without a speaker "pop" or loud click.

Using a raised cosine taper reduces potential high frequency spectral content inherent in any sudden sharp startup transient, which might not be present in the rest of the data due to prior low pass filtering.


Getting a sub-sample resolution

A very cheap (in terms of code size) solution is just to upsample your signal. In matlab, this can be done with interp(y ,ratio). A slightly more complicated solution consists in naively detecting peaks ; and for each peak, fitting a parabola through y[peak - 1], y[peak], y[peak + 1] ; then using the point at which this parabola is maximal as the true peak position.

Regarding peak detection

A bunch of techniques which help:

As suggested by Hilmar, convolving the signal by a Gaussian or Hann window, the width of which is roughly equal to half the minimum interval you want to see between detected peaks. Since temporal accuracy seems essential to your application, make sure that you take into account the time delay introduced by the filtering, though!
Subtract to your signal a median filtered version of itself (with a fairly large observation window) ; and divide the result by a standard-deviation filtered version of itself. This gets rid of trends and allows the thresholds to be expressed in units of standard deviations.
For peak-picking, I formulate that using a "top-hat" filter. Define the top-hat filtered version of your signal as yt[n] = max(y[n - W], y[n - W + 1], ..., y[n + W - 1], y[n + W]) ; and use as peaks the points where y[n] == yt[n] and y[n] > threshold.

All this can be very efficiently implemented in Matlab with a few passes of nlfilter.


Deriving SER & BER for 4QAM, 16QAM and 32QAM

I will try to give a relatively simple answer to a complex question. I will not give any exact expressions for the error rates for two reasons. The first I stated in my initial sentence, the second is that from your formulas I can see that you only looked for (or at least found) approximate expressions.

Judging whether an expression for the error rate is correct - either as an exact expression or as a valid approximation - is complicated by several facts, the most important of which are:

different authors use different functions to express the cdf of the Gaussian distribution. The following relations hold:
$$Q(x)=\frac12\textrm{erfc}\left(\frac{x}{\sqrt{2}}\right)=\frac12\left[1-\textrm{erf}\left(\frac{x}{\sqrt{2}}\right)\right]$$
where erf$(x)$ is called the error function, erfc$(x)$ is the complementary error function, and $Q(x)$ is simply called the Q-function.
The noise power spectral density $N_0$ is usually defined as the one-sided density, but some authors use the two-sided density. See this answer to a related question for a brief discussion of this issue.
The signal energy can be expressed by the energy per symbol, the energy per bit, by the minimum distance between symbols, etc.

For these reasons the same error rate can be expressed by several different formulas and it takes some patience to verify that they are actually all the same.

In the following I will use:

the Q-function $Q(x)$
energy per bit $E_b$
the one-sided noise power spectral density $N_0$
$M=2^k$ is the number of symbols, where each symbol represents $k=\log2 M$ bits
the assumption that all symbols are equally likely

The probability of symbol error $P_s$ for an AWGN channel can be approximated by

$$P_s\approx a(M)Q\left(\frac{d(M,E_b)}{\sqrt{2N_0}}\right)\tag{1}$$

This approximation
is valid at high and moderate SNRs, and is often used in practice. The function $a(M)$ is the average number of nearest neighbors in the symbol constellation. The function $d(M,E_b)$ is the minimum distance between two symbols. It depends of course on the number of symbols $M$ as well as on the energy per bit $E_b$. For even $k$ (i.e. for quadratic constellations) $a(M)$ is given by

$$a(M)=4\left(1-\frac{1}{\sqrt{M}}\right)\tag{2}$$

Let's look at the case $M=16$. The 4 inner points have 4 neighbors, the 8 edge points have 3 neighbors, and the 4 corner points have 2 neighbors, i.e. the average number of neighbors is

$$\frac14\cdot 4+\frac12\cdot 3+\frac14\cdot 2=3$$

which of course equals $a(16)=3$.
In the case of a cross-constellation the formula (2) is not valid anymore. For a cross-constellation with $M=32$ you have 16 inner points with 4 neighbors, 8 edge points with 3 neighbors, and 8 "corner" points with 2 nearest neighbors. So the average number of nearest neighbors is

$$\frac12\cdot 4+\frac14\cdot 3+\frac14\cdot 2=\frac{13}{4}=3.25\tag{3}$$

Now we need to find the function $d(M,E_b)$. Again, for even $k$ (quadratic constellations) a general expression for $d(M,E_b)$ can be derived with relative ease (see e.g. here):

$$d(M,E_b)=\sqrt{\frac{6\log_2 (M) E_b}{M-1}}\tag{4}$$

For the cross-constellation with $M=32$ we first need to compute the average energy per symbol. There are 5 different amplitude values, each contributing to the average energy per symbol (the symbols are at locations $\{a,3a,5a\}+i\{a,3a,5a\}$):

$$E=\frac{4}{32}\cdot 2a^2 +\frac{8}{32}\cdot10a^2 +\frac{4}{32}\cdot18a^2 +\frac{8}{32}\cdot26a^2 +\frac{8}{32}\cdot34a^2=20a^2$$

Since the minimum distance between two symbols is $d=2a$ we get for $d$

$$d=2a=2\sqrt{\frac{E}{20}}=\sqrt{\frac{E}{5}}$$

With $E=\log_2M\cdot E_b=5E_b$ we finally get for the $M=32$ cross-constellation

$$d=\sqrt{E_b}\tag{5}$$

Combining (1), (3), and (5) we get the approximation for the symbol error rate of the $M=32$ cross constellation:

$$P_s\approx 3.25\cdot Q\left(\sqrt{\frac{E_b}{2N_0}}\right)$$

From the above formulas we get for 4-QAM and for 16-QAM

$$P_s\approx 2\cdot Q\left(\sqrt{\frac{2E_b}{N_0}}\right)\quad\textrm{(4-QAM)}\\
P_s\approx 3\cdot Q\left(\sqrt{\frac{4E_b}{5N_0}}\right)\quad\textrm{(16-QAM)}
$$

The probability of bit error can be obtained from the probability of symbol error quite easily if we assume Gray encoding (i.e. nearest neighbors differ by only one bit), and if we assume that errors are only caused by confusion with nearest neighbors (which is a good assumption at moderate to high SNRs). In this case a symbol error causes only one bit error and the bit error probability can be approximated by

$$P_b\approx\frac{1}{\log_2M}P_s$$


Spectral Analysis of a Time Series with Missing Data Points

Given $ \left\{ x \left[ n \right] \right\}_{n \in M} $ where $ M $ is the set of indices given for the samples of $ x \left[ n \right] $.
The trivial solution (Which it would be great to have a faster more efficient solution is what I'm looking for) would be:
$$ \arg \min_{y} \frac{1}{2} \left\| \hat{F}^{T} y - x \right\|_{2}^{2} $$
Where $ \hat{F} $ is formed by subset of columns of the DFT Matrix $ F $ matching the given indices of the samples, $ x $ is the vector of the given samples and $ y $ is the vector of the estimated DFT of the full data of $ x \left[ n \right] $.
The solution is then given by the Pseudo Inverse (Least Squares Solution):
$$ y = { ( \hat{F} \hat{F}^{T} ) }^{-1} \hat{F} x $$
In practice, the matrix will be very poorly conditioned hence solution must be generated using the LS Solution using the SVD.
A sample code is shared on GitHub Repository.
Result of the code:

Remark: Borrowed from Estimate the Discrete Fourier Transform / Series of a Signal with Missing Samples.


Definitions for amplitude spectrum and phase spectrum

Beni I think you want something like this:

samplingFrequency = 1000; % / Hz
signalFrequency = 50;     % / Hz
signalLength = 500;
signal = cos(2*pi*signalFrequency./ samplingFrequency .*[0:1:signalLength - 1]);
[hh ff] = freqz(signal, 1, 2048, samplingFrequency);
figure(10);
plot(ff,10*log10(abs(hh)));
ylabel('Amplitude in log scale');
xlabel('Frequency / Hz');
figure(20);
plot(ff, angle(hh));
ylabel('Phase');
xlabel('Frequency / Hz');

You do have MATLAB correct?... if you do try it and let me know how else I can help.

Well, a couple things. First off, remember, since this is a digital signal, it means of course that it was sampled at some frequency $samplingFrequency= 128 or 127$ in your case. The physical frequency of your signal in this case seems to be $signalFrequency = 21$.

The ratio of $\frac{signalFrequency}{samplingFrequency}$ is equal to the $f_0$ in your case. This is called the 'digital frequency'.

For your problem it sounds like he also wants you to fix the fft size, (change the 2048 in my code to 128) and the rest of the numbers as I described and run the code. Does he mention how long he wants your signal to be?


How to calculate the frequency of a signal without knowing the sampling frequency using Matlab

If you do not have the information about the sampling frequency $F_s$ of your digital data, the best option is to talk about dimensionless relative frequencies $f$, or reduced frequency. The frequencies you observe on periodograms will be $f = F/F_s$, where $F$ would be the true frequency (hoping you have no aliasing). This amounts to saying that your sampling period is $1$ (dimensionless) and that the maximum observable frequency in your signal is $1/2$.

However, it is likely that the actual sampling frequency can be obtained from the data file, the  experiment or sensor, a recorded phenomenon with known frequency (like the 50 or 60 Hz power) or the person who gave you that signal.


The fact that sampling in one domain corresponds to aliasing in the other domain is an important property of the Fourier transform. Of course, this can also be shown formally. If $X[k]$ is the length $N$ DFT of $x[n]$, we have

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{1}$$

Furthermore, $H(z)$ is related to the impulse response by

$$H(z)=\sum_{m=-\infty}^{\infty}h[m]z^{-m}\tag{2}$$

If we define

$$X[k]=H(e^{j2\pi k/N})=\sum_{m=-\infty}^{\infty}h[m]e^{-j2\pi mk/N}\tag{3}$$

we can determine $x[n]$ as follows. From $(1)$ and $(3)$ we have

$$\begin{align}x[n]&=\frac{1}{N}\sum_{k=0}^{N-1}\sum_{m=-\infty}^{\infty}h[m]e^{-j2\pi mk/N}e^{j2\pi nk/N}\\&=\frac{1}{N}\sum_{m=-\infty}^{\infty}h[m]\sum_{k=0}^{N-1}e^{j2\pi(n-m)k/N}\tag{4}\end{align}$$

The last sum in $(4)$ equals $N$ for $n-m=lN$ for integer $l$, and zero otherwise. Consequently, this last sum can be rewritten as

$$\sum_{k=0}^{N-1}e^{j2\pi(n-m)k/N}=N\sum_{l=-\infty}^{\infty}\delta[n-m-lN]$$

and $(4)$ becomes

$$\begin{align}x[n]&=\sum_{m=-\infty}^{\infty}h[m]\sum_{l=-\infty}^{\infty}\delta[n-m-lN]\\&=\sum_{l=-\infty}^{\infty}\sum_{m=-\infty}^{\infty}h[m]\delta[n-m-lN]\\&=\sum_{l=-\infty}^{\infty}h[n-lN]\tag{5}\end{align}$$

where I've used the fact that the sum over $m$ is the convolution of $h[n]$ with $\delta[n-lN]$, and consequently equals $h[n-lN]$.

Eq. $(5)$ establishes the result that the IDFT of an equidistantly sampled discrete-time Fourier transform (DTFT) equals an aliased version of the sequence corresponding to the given DTFT.


What kind of information can we obtain from curvature for exemplar-based inpainting?

I haven't read the paper before, but the rationale is state in III.A:

  The filling order is very crucial which directly affects the quality of image inpainting. However, during filling procedure we find the data terms may be zeros for many pixels, which make the priorities being close to zeros, although the confidence terms may be very big. This would result in incorrect filling order.[...] by introducing the curvature term,
  we effectively avoid the priority being close to zero while the
  data term is zero.

The foundation for the method used in this paper is in Object removal by exemplar-based
inpainting, so be sure to read that first.


In fact Welch was a good idea. End of post. Problem solved.

# Loop for FFT data
for dataset in [fft1]:
    dataset = np.asarray(dataset)
    freqs, psd = welch(dataset, fs=266336/300, window='hamming', nperseg=8192)
    plt.semilogy(freqs, psd/dataset.size**2, color='r')

for dataset2 in [fft2]:
    dataset2 = np.asarray(dataset2)
    freqs2, psd2 = welch(dataset2, fs=266336/300, window='hamming', nperseg=8192)
    plt.semilogy(freqs2, psd2/dataset2.size**2, color='b')


Trouble understanding an article about MFCC

So the output of the Discrete Cosine Transform gives a 13-vector. To this, they append the log mean energy. That gives a 14-vector.

They then subtract from this the previous 14-vector to yield another 14-vector, which they concatenate with the existing 14-vector.

That gives a 28-vector.

I suspect the log mean energy is just:
$$
{\tt LME} = \log\left(\frac{1}{N} \sum_{N} x_n^2 \right)
$$
where $N$ is the length of the frame used.


Dealing with spectral leakage between two frequency bins

Yes, concatenating the data will greatly reduce the spectral leakage problem.  There will still be spectral leakage of course, but there will be more bins in between your signals-of-interest, so they won't leak as much into each other.

And yes, windowing should help too.  You are already using a "rectangular" window.  I would suggest trying other windows such as Hamming.


expression for the FT of the frequency response of a system

Try writing your equation (2) as $Y(e^{j\omega})=W(e^{j(\omega + \pi)})=X(e^{j\omega})H_1(e^{j(\omega + \pi)})$, and now try solving for $H(e^{j\omega})$. Remember that phase shifting by $\pi$ and $-\pi$ gives the same result.


Calculating an image's fourier spectrum by hand?

Your particular case can be solved without a calculator. Fourier transforms are linear operations. You can decompose your input as:

$$3 \left( \begin{array}{cccc} 0 & 0 & 1 & 1 \\
0 & 0 & 1 & 1 \\
0 & 0 & 1 & 1 \\
0 & 0 & 1 & 1 \end{array} \right) + \left( \begin{array}{cccc} 0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \end{array} \right) $$

The first matrix is constant on the Y axis, so its 2D FFT is only the 1D-FFT of one row, times 4, and 0 everywhere else. The 1D FFT of $(0, 0, 1, 1)$ is $(2, -1+j, 0, -1-j)$, so this gives you, for the first term:

$$3 \times 4 \times \left( \begin{array}{cccc} 2 & -1+j & 0 & -1-j \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \end{array} \right)$$

As for the second term, there's only one value so there's only one component in the sum - it's just one of the basis functions, the checkerboard pattern:

$$\left( \begin{array}{cccc} 1 & -j & -1 & j \\
-j & -1 & j & 1 \\
-1 & j & 1 & -j \\
j & 1 & -j & -1 \end{array} \right)$$

Summing these two, you get:

$$\left( \begin{array}{cccc} 25 & -12+11j & -1 & -12-11j \\
-j & -1 & j & 1 \\
-1 & j & 1 & -j \\
j & 1 & -j & -1 \end{array} \right)$$

Note that this differs from Paul R's answer because of the normalization and sign conventions nobody agrees with... Divide by 4 and take the conjugate and you'll find the same result as his.


What is the difference between a controller and a compensator?

As the book's title is control theory; a controller is the general broad name given to any piece of hardware (or sowftware) that aims to control some processes, in your aerodynamics field probably the flight process; by menas of receiving signals from sensors and sending commands to propellers - actuators, and in between is the control hardware or software...

A compensator is a specific control-box-tool to correct (or compensate) some lacking aspect of a system property under concern, such as phase, lag, magnitude etc... It is generally useful as part of a bigger control system, rather than itself being a controller except some simple applications.

Obviously there are many different control system examples in addition to a compensator...