The official documentation of that block:

  An SNR estimator for M-PSK signals that uses 2nd (M2) and 4th
  (M4) order moments. This estimator uses knowledge of the
  kurtosis of the signal ($k_a)$ and noise ($k_w$) to make its
  estimation. We use Beaulieu's approximations here to M-PSK
  signals and AWGN channels such that $k_a=1$ and $k_w=2$. These
  approximations significantly reduce the complexity of the
  calculations (and computations) required.

  Reference:
  D. R. Pauluzzi and N. C. Beaulieu, "A comparison of SNR
  estimation techniques for the AWGN channel," IEEE
  Trans. Communications, Vol. 48, No. 10, pp. 1681-1691, 2000.

The documentation of the constructor vaguely hints:

  alpha   the update rate of internal running average calculations."

I.e. probably a single-tap IIR ("exponentially weighted moving average") with feeback coefficient alpha.

Verifying that with the source code:

for(int i = 0; i < noutput_items; i++) {
    double y1 = abs(input[i])*abs(input[i]);
    d_y1 = d_alpha*y1 + d_beta*d_y1;

    double y2 = abs(input[i])*abs(input[i])*abs(input[i])*abs(input[i]);
    d_y2 = d_alpha*y2 + d_beta*d_y2;
}

Yes, that alpha is just the coefficient of a single-tap IIR (further up you'll find d_beta = 1 - d_alpha).

So, this estimator should do exactly what it promises: for a zero-mean signal (where the variance is simply $\text{Var }X = \text E(X^2)$), it calculates the second and fourth central moment and then just

double
mpsk_snr_est_m2m4::snr()
{
  double y1_2 = d_y1*d_y1;
  d_signal = sqrt(2*y1_2 - d_y2);
  d_noise = d_y1 - sqrt(2*y1_2 - d_y2);
  return 10.0*log10(d_signal / d_noise);
}

So, it does the usual kurtosis-based estimate.

By the way, I'd fully agree the parameter documentation – to put it friendly – leaves a bit of room for improvement. Hit us / me up if you'd be willing to contribute a bit of a documentation improvement!


As you suspect, this paper says that zero-lag Butterworth filtering is obtained by passing the signal through the filter twice: once in the forwards direction and once in reverse.


I think you are having trouble constructing your frequency axis properly. Once you do this, you can do a simple peak pick. I have re-written the code for you:

Fs = 48000;
x = dataSet;            % 256 samples data set
fftLength = length(x);  % Always make sure to be at least as long as your data.
xdft = fft(x,fftLength);
maxAmp = max(abs(xdft));

freq = [0:fftLength-1].*(Fs/fftLength); % This is your total freq-axis
freqsYouCareAbout = freq(freq < Fs/2);  % You only care about either the pos or neg
                                        % frequencies, since they are redundant for
                                        % a real signal.

xdftYouCareAbout = abs(xdft(1:round(fftLength/2))); % Take the absolute magnitude.

[maxVal, index] = max(xdftYouCareAbout); % maxVal is your (un-normalized) maximum amplitude

maxFreq = freqsYouCareAbout(index); % This is the frequency of your dominant signal.


$N = kTB$ is the expression for the total thermal noise power within bandwidth $B$. $kT$ (where $k$ is Boltzmann's constant and $T$ is the temperature in Kelvin) is the noise density, meaning power/Hz. It is also a white noise process thus this single constant (assuming a fixed temperature) is the noise density at all frequencies. ($kT$ at room temperature is approximately -174 dBm/Hz, a number worthy of committing to memory).

If you want to compute the total power within a bandwidth $B$, you simply multiply by $B$ to get $kTB$.

The second part of your question is not clear; typically we would compute the noise power within the signal's bandwidth in the process of estimating the signal to noise ratio as the well designed receiver would have filtered out the (amplified) noise power beyond the signal's frequency bandwidth. So $B$ in this case should be the signal's frequency bandwidth and not an arbitrarily "large bandwidth".


Is sift a good way to extract features from an image?

What you are supposed to do when matching a template to an image using sift is to run sift against your template and then look for those sift features in that arrangement in your scene.

Rule of thumb: Compare like to like.

Sift(Template) Contained Within Sift(Image)

You cannot tune Sift to extract the "features you want" Sift *uses invariant measures to find features that it believes are most salient. (*Chris A)

If you want it to find particular features from your template make smaller templates of those features and search for them.

Just a suggestion. I do not know exactly what type of template you are trying to fit to what type of scene or image.


What does one call this method of grid fitting?

I don't think this has a name, but you could just apply a clustering algorithm to the X-coordinates of the red dots and separately to the Y-coordinates. For example, in Mathematica:

img = Import["http://i.stack.imgur.com/DfhTN.png"];
blobs = ComponentMeasurements[ColorNegate@Binarize[img], "Centroid"][[
   All, 2]];
{meanX, meanY} = Mean /@ FindClusters[#] & /@ Transpose[blobs];

This will find "clusters" in the X- and Y-coordinates of the red dots:

Show[img,
 Graphics[{Line[{{#, 0}, {#, 480}}] & /@ meanX,
   Line[{{0, #}, {640, #}}] & /@ meanY}]]

(If you want to implement the clustering yourself, k-means is a good starting point: In a nutshell, you start with random cluster centers; then assign each point to the nearest cluster center; then update the cluster centers to the mean of the values assigned to them. Repeat until convergence, split clusters if the assigned points spread too far.)

The grid you're looking for is just the lines in between these clusters:

inBetween = Mean /@ Partition[Sort[#], 2, 1] &;
Show[img,
 Graphics[{Line[{{#, 0}, {#, 480}}] & /@ inBetween[meanX],
   Line[{{0, #}, {640, #}}] & /@ inBetween[meanY]}]]


For Direct-Sequence Spread Spectrum, can we reduce the transmit power of the signal to maintain a fixed BER?

BER depends on your E_b/N0. You need to keep that constant if you want constant BER.
Now, you don't say how your comparing your DSSS system, but it would make sense to assume that your symbol rate constant, thereby increasing your bandwidth by the spreading factor, and keeping the duration of a symbol constant.
In that case, Energy is still the product of power and duration, and the latter stayed constant.


pyo: Retrieve amplitudes of audio signal frequenciy ranges

pyo is not really the issue here at all. The question still stands if you were to substitute it with a million other packages of similar functionality.

What is at the centre of this question is "How to build a classifier to distinguish between different types of hits on a drum". Once you have that sorted, you can then implement it with whatever package you like, including pyo.

One of the simplest classifiers you can use for this task is the Nearest Neighbour Classifier.

The end effect of this classifier is not too different than your current thinking: I will split the frequencies of the sound and then try to work with if-then-else constructs, the final outcome of which would be the classification. Every if-then-else, puts a binary decision boundary on the outcome. For example, if you monitor the audio volume then you can create a rule that tells you whether there is something playing or not (if the volume is above a threshold, then something is playing). The nearest neighbour classifier over spectral features (of the FFT), does exactly the same thing but puts $N_{DFT}$ if-then-else on a given frame (it is not exactly if-then-else, but the end-effect is the same).

To get it to work, you would need a few representative recordings of each drum shot. For example, 20 different rimshots, 20 different straight hits. If you pass each one through the Discrete Fourier Transform (DFT) (of which the FFT is a quick implementation of) this will give you $c_{k,i,N_{DFT}}$ coefficients. That would be one set of $N_{DFT}$ coefficients for the $i^{th}$ recording of the $k^{th}$ type of drum hit.

Find the centroids (means of each DFT coefficient) that correspond to each drum hit and this will give you $k$ points in feature space that represent what each drum hit tends to sound like (obviously, here $k=2$).

To figure out what sort of drum hit has just been recorded (a drum hit you have not encountered during the "training" phase), take the waveform of it, pass it through DFT, take its DFT coefficients and work out which centroid they are closer to.

Doing this in real time might be a bit more challenging as you would have to establish an audio "loop" that continuously records frames from the audio interface, transforms them and passes them over to the classifier to make a decision.

Eventually, you will probably also need a way to describe a "third class" for sounds that are not drum hits and in which case you would like your classifier to remain "silent".

Overall, the system will also require some tuning to bring it as close as possible to a real time identification of drum hits, presumably, from a clean audio source, but it would be a start.

Hope this helps.


Viterbi Decoder - Decoding Continuously Encoded Stream Partially

Since the decoder that you are working with is a black box, yes, adding prefixes and postfixes to the stream to "window-ize" it is appropriate.  The pre and postfixes need to be the actual data from the previous and post windows, and they need to be at least as long as the traceback depth.  If you do this and get rid of the result bits that correspond to the prefixes and postfixes, then you should do fine.

In a normal situation you would not need to add any postfixes or prefixes.  There are only two things you need to do to handle continuous streams of Viterbi-encoded data.

1) occasionally check the state metrics and lower them when all of them get "high".  This works because the only thing that matters is the difference between the state metrics, not their absolute value.  An easy way to do this is to check if all of them are equal to or above half the maximum state metric value.  If they are, subtract half the maximum value from all of them.

2) You have to do tracebacks while you are decoding.  This is do-able in both FPGA's and microprocessors (it's actually quite easy in microprocessors).  Every $N$ transitions (pick an $N$ that works for your architecture) you do a traceback and output $N$ bits.  You have to go back through $M$ bits ($M$ being equal to the traceback depth) before you get to the $N$ bits that you will output.


From what I read there are certain quantisation matrices for different applications, Adobe Photoshop has like 15 or something. The idea is that although there are error based mathematical calculations for these matrices the best way to deduce a good quantisation matrix is by simply using a human eye (in the case of JPEGs)I haven't seen any quantisation matrices that don't follow the general trend of having heavy penalties(High number) on the lower right frequencies, although some images have a large high frequency coefficients the lower frequencies are generally larger.


Why 8-PSK modulation has higher Eb/N0 than QPSK or BPSK?

example - at a distance of 10 m (transmitter to receiver) SNR for QPSK was 22.9 and for 8-PSK it was 26. Does it depend on transmitting more bits symbols per cycle (8-PSK transmits 3 while QPSK/BPSK only 2)?

This makes no sense, physically: The attenuation through the channel is the same for all these signals, since you'd put them all through the same pulse shape to have a fair comparison.
The fact that they don't transport the same number of bits per symbol doesn't matter to physics.
Maybe you're confusing $E_b/N_0$ with SNR?

Ah yes, you were confusing $E_b/N_0$ with SNR. Now, you didn't compare fairly. In BPSK, a bit gets the energy of a full symbol. That's the transmit power times the symbol period.
Since you need to compare fairly, you need to keep both the transmit power constant (obvious) and either the bit rate constant, or the symbol rate constant:

constant power, constant bit rate: Now your 8-PSK transmit signal has only 1/3 of the symbol rate of your BPSK signal. That means narrower bandwidth, 1/3 of BSPK. Which means 1/3 of the noise. You cheated here by reducing the bandwidth.
constant power, constant bandwidth: your system now divides the transmit power on three times as many bits as BPSK, so $E_B/N_0$ actually is reduced, since within the same bandwidth, the same original noise occurs.

Note that in (almost) any sensible comparison of transceiver systems, you keep the $E_b/N_0$ constant and vary the bandwidth or necessary transmit power – your receiver cares less about the SNR than it cares about how likely it is that a bit gets mis-decided, which is a function of the energy per bit. This becomes slightly more obvious when one considers what interleavers and channel decoders do: they distribute noise energy that were concentrated on individual bits or symbol-bit-groups over a larger duration, and then correct these less-than-perfect bits using the knowledge from all the other bits.


Is there an equivalent time-domain operation to the spectrum averaging in Welch's method?

Yes, you can average in time lag (autocorrelation) as well as in frequency. You can also combine the two domains, as long as it is reasonable to treat the data as being stationary.

A. H. Nuttall and G. C. Carter, "Spectral estimation using combined time and lag weighting," in Proceedings of the IEEE, vol. 70, no. 9, pp. 1115-1125, Sept. 1982.


What is the mean power of a complex random variable?

The correct expression for the mean power of a complex random variable $\alpha=x+jy$ is
$$
\begin{align}
\bar P &= \operatorname E\left[\alpha \alpha^*\right]\\
&= \operatorname E\left[x^2 + y^2\right]\\
&= \operatorname E\left[x^2\right] + \operatorname E\left[y^2\right]\\
&= \bar P_\mathrm{x} + \bar P_\mathrm{y}
\end{align}
$$

In other words, the mean power of a complex random variable is the sum of the mean powers of its real and imaginary part, respectively. In contrast, expression no. 1 from your question evaluates to
$$
\begin{align}
\operatorname E\left[\alpha^2\right] &= \operatorname E\left[x^2 + j2xy - y^2\right]\\
&= \operatorname E\left[x^2\right] - \operatorname E\left[y^2\right] + \operatorname E\left[j2xy\right]\\
&= \bar P_\mathrm{x} - \bar P_\mathrm{y} + j2\operatorname E\left[xy\right]
\end{align}
$$
Note that for the special case of $x$ and $y$ being uncorrelated $\operatorname E\left[xy\right]=E\left[x\right]E\left[y\right]$. If, in addition, $x$ or $y$ has zero mean $E\left[xy\right]=0$.


Detecting a curve with certain characteristics in a spectrogram

This answer might sound hard initially, but once understood it is intuitive and certainly deals well with the noise.

One smart way to model such constraints is either through HMM's or through MRFs. The design is similar, but for HMM I would give you some hint:

Prepare your transition matrix in such a way that Viterbi decoding would provide you the resulting curve. For example, the first state will be at 0 with certain variance, since it starts there. Then it is only allowed to increase, within a given tolerance. This tolerance would be obtained from a Gaussian prior. Say, the slope of the line (thus the variance of the Gaussian) would range from 1 to 5. Finally, you should only allow for a decrease or keeping constant. Similarly, for every other constraint you could do this. Here allowance means that the probability of transition from the constant state to an increasing state is 0. Otherwise you have a probability (might be due to any distribution, prior etc.)

You could also make use of Hierarchical HMMs to elaborate your model with such tracking parameters:

http://en.wikipedia.org/wiki/Hierarchical_hidden_Markov_model
http://www.cs.princeton.edu/courses/archive/spr06/cos598C/papers/FineSingerTishby1998.pdf


Wiener filtering/deconvolution for non-stationary noise

The Wiener Filter can written as a Least Squares problem (See How Is the Formula for the Wiener Deconvolution Derived?).
In you case, Since noise is basically with different STD per sample you may use the simple extension of Weighted Least Squares on top of the terms for the Wiener Filter.
Once you do that you'll have the optimal solution.


Finding autocorrelation matrix of an autoregressive process AR(1)

Hi: You're expression will be a 3 by 3 matrix with the 1's along the diagonal and the correlations in the respective places.

$\begin{matrix}  1 & \sigma_{12} & \sigma_{13}  \\
                 \sigma_{21} & 1 & \sigma_{23}  \\
                 \sigma_{31} & \sigma_{32} & 1
\end{matrix}$

The 12 and 21 elements are $\rho$   because the lag one autocorrelation of an AR(1) with parameter $\rho$ is $\rho$.

The 23 and 32 elements are also $\rho$  for the same reason as above:  The lag one autocorrelation of an AR(1) with parameter $\rho$ is $\rho$.

The 13 and 31 elements are $\rho^2 $  because the lag 2 autocorrelation of an AR(1) with parameter $\rho$ is $\rho^2$.

The transfer function implies that $\rho = 0.7$. This is because

$ y_t = \rho \times y_{t-1} + \epsilon_t$ has transfer function $H(z) =
\frac{1}{(1-\rho z^{-1})}$.


Can you count signal energy equivalent in frequency domain?

Parseval's theorem applied to the signal itself:

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

So, if you sum the magnitude squared of the frequency coefficients of the DFT, you get the energy of your time-domain vector.

If you want power (i.e. energy per sample), divide both sides by $N$.


How to calculate the energy in the dirac delta function signal?

\begin{align}
&\int_{\omega_1}^{\omega_2}|\delta (\omega-\omega^*)|^2.d\omega \\
&= \int_{\omega_1}^{\omega_2}|\delta (\omega-\omega^*)|.|\delta (\omega-\omega^*)|.d\omega \\
&= \int_{\omega_1}^{\omega_2}\delta (\omega-\omega^*).\delta (\omega-\omega^*).d\omega \\
&= \int_{\omega_1}^{\omega_2}\delta (\omega^*-\omega^*).d\omega \\
&= \int_{\omega_1}^{\omega_2}\delta (0).d\omega \\
&= \delta (0).\Delta \omega
\end{align}
Which is undefined! See this for better clarification: https://en.wikipedia.org/wiki/Distribution_(mathematics)#Adding_and_multiplying_distributions

What you have to understand is that $\delta$ is not a function, it is a distribution. What you learn in introductory DSP courses is just an "engineering trick" so that things will work the way they should (you will understand what that means once you study sampling and analog to digital conversion).

EDIT: Sorry I misinterpretation the question. (Thanks Jazzmaniac for pointing that out)

As far as I know the $\delta$ "function" in the continuous domain is not in $L^2$, it isn't a square-integrable function.

Think about your question in the time domain. If you have 2 symmetric $\delta$'s in the frequency domain, you will have a INFINITLY lasting sinusoid. It's energy is infinite, i.e, you can't associate a value to it (undefined, remember that $\infty$ is not a number). If you have a single $\delta(f)$, centered at zero, you will have a complex exponential in the time domain, and it's energy is also undefined (unbounded since $|e^x|^2=1$.)


If the cosine function is periodic, why does it have a Fourier Transform?

Indeed there are two things you have to know.

First, it can be shown that the continuous-time Fourier transform can be obtained from the continuous-time Fourier series by letting the period $T$ go to infinity.

Second, formally speaking the Fourier transform integral for periodic signals do not converge, hence do not exist. The solution is a generalisation of the Fourier transform by the use of Dirac impulse functions.

The result is an interpretation that the Fourier transform of periodic functions is a sum of scaled Dirac impulses at the Fourier harmonic frequencies, and the scale being the corresponding Fourier series coefficients.


In reverse order,

Yes the sampling frequency needs to be a multiple of $f_0$.
The number of cycles of the sinusoid in $N$ samples needs to be an
integer, not necessarily a prime number.
If you do 3. and 2., all the energy will be concentrated in two bins
(corresponding to $100$ Hz and $-100$ Hz.

Be careful, though. If you have a sinusoid $x(t)=\sin(2\pi (100) t)$ at 100 Hz
and you sample it at 1000 Hz, you will have $10$ samples per cycle of the
waveform. With $x[0] = 0$ and $x[9] \neq 0$, you can do an "FFT"
of length $10$ if you like, and all the energy will be in $X[1]$ and $X[9]$.
I put FFT in quotation marks because at such short lengths, a plain
vanilla brute-force DFT calculation might be almost as fast. Note also
that $x[10] = 0$ is not part of the DFT calculation.
Alternatively, you can use $N = 10n$ samples and do an FFT of length $1000$  (say) on samples $x[0]$ through $x[999]$.  In this case, the
energy will all be in $X[10]$ and $X[1014]$.
Regardless of what some other answers will likely tell you,

  DO NOT zero-pad the length $1000$ FFT to length $1024$ to get an
  "efficient" implementation.

If you do follow the advice on zero-padding, your length-$1024$ FFT
will NOT have the energy concentrated in two bins as you desire.


Why many DACs uses series of half band filters for interpolation instead of single one?

The reason cascaded half-band filters, for interpolation applications, should have different numbers of taps can be found at:
https://www.dsprelated.com/showarticle/903.php


Given a signal $x[n]$ with the nonzero range of $N_1 \leq n \leq N_2$ and the signal $h[n]$ with the nonzero range $M_1 \leq n \leq M_2$, then the nonzero range of their comvolution $y[n] = x[n] \star h[n]$ is $N_1+M_1 \leq n \leq N_2 + M_2$.

Applying this to your case yields the following:
For the signal $x[n]$ with $n=0$ at third position the nonzero range is $-2 \leq n \leq 2$ and for $h[n]$ with $n=0$ at the second position the nonzero range is $-1 \leq n \leq 2$ which therefore yields the nonzero range of their convolution $y[n]$ is:
$$ -3 \leq n \leq 4 $$

which means that the $n=0$ for $y[n]$ is at the fourth position.


How do I apply the Nyquist criterion on the following spectrum to see if it satisfies the criterion or not?

The Nyquist criterion for avoiding intersymbol interference (ISI) requires the pulse to have zeros at multiples of the symbol interval $T$, such that symbols do not interfere with each other. The time domain requirement is

$$h(t)=\begin{cases}1,&t=0\\0,&t=nT\end{cases}\tag{1}$$

In the frequency domain, criterion $(1)$ can be formulated as

$$\sum_{k=-\infty}^{\infty}H\left(f-\frac{k}{T}\right)=T\tag{2}$$

where $H(f)$ is the Fourier transform of $h(t)$. I.e., the shifted spectra (shifted by the symbol rate) must add up to a constant.

I leave it up to you to decide whether or not the given pulse satisfies $(1)$ and $(2)$.

I suggest you go through your lectures notes because it would be surprising if the Nyquist criterion has never been mentioned and you get an exercise like this.


Conceptual question on entropy and its relation to information

The basic idea behind maximum entropy models is that you want to make the least assumption about the data. This is considered equivalent to retaining as much unpredictability as possible, as quantified by entropy. For more information I refer you to this Wikipedia article.

Information and unpredictability are closely related. If you think of the signal as a stochastic process, its information content is determined by unpredictable it is. At one extreme, if the signal tends towards being deterministic, you can tell what it will be any an arbitrary time so there is no value in sampling the signal; it has no information content. Entropy formalizes this notion.

To find the density that maximizes the entropy we have to use the calculus of variations, and that involves determining several coefficients called Lagrange multipliers. This is where the moments come in: they're the constraints that determine the multipliers. We set the $n$'th absolute moment for n=1..N to constants and we get the corresponding maximum entropy distribution. For example, the maximum entropy distribution whose first two moments are $\mu$ and $\sigma^2$ is the Gaussian distribution $N(\mu, \sigma^2)$. For details of the derivation please refer to a textbook.

I haven't read the paper but I'm guessing the rationale is that if you're trying to fit a model, you want to capture as much information of the signal as possible. Equivalently, you want to leave as little residual information as possible on the table; i.e., in the error. For details see 3.4 Minimum error entropy algorithm in Information Theoretic Learning: Renyi's Entropy and Kernel Perspectives.


Destructive interference on the Autocorrelation of the time signal of periodic CPFSK-signals.

After thinking for a while i came up with my own proof. I hope it is helpful to you and correct!

Definition: Frame-Structure

The frame structure used in the following is one of one-periodic signals, with a constant header and a uniform distributed payload. The ratio between payload and total frame is given as $$ p_{\%} = \frac{|payload|}{|frame|} $$
the ratio of samples in the header vs samples per frame.

Example:
$$ |HHH|012|HHH|223|HHH|978|...$$
,here: $p_{\%}=0.5$ .

Proof:

The complex envelope of an FSK-signal is

\begin{array}{cc}
x\left(t\right)=e^{i\left[2\pi\Delta FT\left(\varphi\left(j,T\right)+\left(\frac{t}{T}-j\right)d\left(j\right)\right)\right]} & ,\text{}jT\leq t\leq\left(j+1\right)T\end{array}

with

$$\varphi\left(j,T\right):=2\pi\Delta FT\sum_{l=0}^{j-1}d\left(l\right)+\varphi_{0}$$

and the modulation index $$\eta:=2\Delta FT.$$

We introduce the model by having a look at the ACF of truely periodic, i.e. $p_{\%}\left(x\right)=1 $, discrete signals. In a second step we generalize to partly periodic signals with $p_{\%}\left(x\right)<1$. Let $x$ be a digital FSK signal of length N+1, that is, $x$ is a sampled version of the complex envelope above. Then the ACF of x is defined by

\begin{array}{ccc}
\rho\left[x\right]\left(\tau\right) & = & \sum_{k=0}^{N-\tau}x\left(k\right)x^{*}\left(k+\tau\right)\\
 & = & \sum_{k=0}^{N-\tau}e^{i\left[\pi\eta\left(\left(\frac{k}{T}-j_{k}\right)d\left(j_{k}\right)+\varphi\left(j_{k},T\right)\right)\right]}e^{-i\left[\pi\eta\left(\left(\frac{k+\tau}{T}-\left(j_{k+\tau}\right)\right)d\left(j_{k+\tau}\right)+\varphi\left(j_{k+\tau},T\right)\right)\right]}
\end{array}

Let $x$ be $\tau_{0}$-periodic, consequently $\tau_{0}$ is an integer-multiple of the symbol duration T measured in samples. Given our data model, we can expect to see a distinct peak in the ACF feature of such a signal. Let $j_{k}:=\left\lfloor \frac{k}{T}\right\rfloor$  be the index of the symbol being modulated at point of time $k\in\mathbb{N}$. While the signal $x$ is $\tau_{0}$-periodic, the series of symbols is $\frac{\tau_{0}}{T}$-periodic, i.e. $j_{k+\tau_{0}}=\left\lfloor \frac{k+\tau_{0}}{T}\right\rfloor =j_{k}+\frac{\tau_{0}}{T}$. Therefore the symbols of a signal are repeated periodically, i.e. $d\left(j_{k+\tau_{0}}\right)=d\left(j_{k}+\frac{\tau_{0}}{T}\right)=d\left(j_{k}\right)$. Hence we can simplify the ACF of x at lag $\tau_{0}$ to

\begin{array}{cc}
\rho\left[x\right]\left(\tau_{0}\right)\\
= & \sum_{k=0}^{N-\tau_{0}}e^{i\left[\pi\eta\left(\left(\frac{k}{T}-j_{k}-\frac{k+\tau_{0}}{T}+j_{k}+\frac{\tau_{0}}{T}\right)d\left(j_{k}\right)\right)\right]}e^{i\left[\pi\eta\left(\sum_{l=0}^{j_{k}-1}d\left(l\right)-\sum_{l=0}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)\right]}\\
= & \sum_{k=0}^{N-\tau_{0}}e^{i\left[\pi\eta\left(\left(-\frac{\tau_{0}}{T}+\frac{\tau_{0}}{T}\right)d\left(j_{k}\right)\right)\right]}e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}\\
= & \sum_{k=0}^{N-\tau_{0}}e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}.
\end{array}

The value $\rho\left[x\right]\left(\tau_{0}\right)$ is decisively influenced by the term $p_{m}:=-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)$. It contains exactly one period of the cyclic symbol stream and it shifts cyclically with time and, thus $p_{m}$ is constant. Accordingly we can simplify the ACF feature for periodical signals to

$$\left|\rho\left[x\right]\left(\tau_{0}\right)\right|=\left|\sum_{k=0}^{N-\tau_{0}}e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}\right|=\left(N-\tau_{0}+1\right).$$

The complex value of the ACF of a truely periodic signal is independent of modulation index and symbol stream.

In the second step we assume a partly-periodic signal x of ratio $p_{\%}$ with lag $\tau_{0}$. The ACF of a uniformly distributed series of samples converges to zero at all non-zero lags. Thus, according to our model the autocorrelation of the data part of the frames all but vanishes. To simplify the analysis, we introduce an indicator function b; indicating wether a sample is part of the periodic share of a signal:
$$\begin{array}{ccc}
b\left(k\right) & := & \begin{cases}
1 & k\in\text{header-part of x}\\
0 & \text{else}
\end{cases}\end{array}$$

We get

\begin{array}{cc}
\rho\left[x\right]\left(\tau_{0}\right)\\
= & \sum_{k=0}^{N-\tau_{0}}b\left(k\right)e^{i\left[\pi\eta\left(\left(\frac{k}{T}-j_{k}-\frac{k+\tau_{0}}{T}+j_{k}+\frac{\tau_{0}}{T}\right)d\left(j_{k}\right)\right)\right]}e^{i\left[\pi\eta\left(\sum_{l=0}^{j_{k}-1}d\left(l\right)-\sum_{l=0}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)\right]}\\
= & \sum_{k=0}^{N-\tau_{0}}b\left(k\right)e^{i\left[\pi\eta\left(\left(-\frac{\tau_{0}}{T}+\frac{\tau_{0}}{T}\right)d\left(j_{k}\right)\right)\right]}e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}\\
= & \sum_{k=0}^{N-\tau_{0}}b\left(k\right)e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}.
\end{array}

The term $p_{m}$ (see [eq:ACF_FSK_PhaseMemory]) is not constant like it was above. As the frames are not truely periodic, the data-part of the frames introduce symbols randomly choosen from a uniform distribution. We split the results by modulation index:

If the modulation index $\eta$ is integer, $\eta p_{m}\in\mathbb{Z}$ and we get $e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}=1.$

Thus,
$$
\left|\rho\left[x\right]\left(\tau_{0}\right)\right|=p_{\%}\left(N-\tau_{0}+1\right).$$

If $\eta$ is not integer,

$$\left|e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}\right|=1$$

but the argument of the summands in ${arg}\left(e^{\pi\eta i\left(-\sum_{l=j_{k}}^{j_{k}+\frac{\tau_{0}}{T}-1}d\left(l\right)\right)}\right)$

is not constant, as it depends on the random data frame, and thus  consists of a sum of complex exponentials of differing phase. The sum decreases compared to the magnitude of the sum of an integer modulation index signal. I conclude that a random walk is in order to determine the effects of the cancellation.


How can I prove stability of a biquad filter with non-zero initial conditions

Tricky problem. Not sure I can answer this but here are a few pointers:

Direct Form II is the worst biquad for audio processing. The transfer function between the input and the state is given just by the poles. The gain can be really large, I have seen gains of in excess of a 100 dB for reasonable audio filters.
This makes real time switching of Direct Form II very tricky. The magnitude of your state can change 10s of dB with seemingly small overall filter changes. This creates massive clicks and pops and occasionally unstable filters.
Much better choices are Direct Form I or transposed Form II. For these filter types the state variables are always in the same order of magnitude than input and output. Much less potential for loud pops and clicks or instability
There is not "sure fire" way of updating an IIR filter in real time without the potential for audible artifacts. Some potential methods are running two filters in parallel and slowly cross fading the output or updating pole and zero locations really slowly with very small steps at the time
If you must update in real time, make sure that you don't update one coefficient at a time but ideally both poles and zeros in one go. At least do the poles all at ones (i.e. always update a1 and a2 at the same time)
Update lower gain first. If the pole section of the new filter has lower gain, update the poles first and then the zeros. If the pole section of the filter has higher gain, update zeros first and then poles
If you need to stick with Direct Form II, it's probably best to zero out the state during update. This will still create an artifact but it's probably less than if you don't and at least it guarantees stability.


The sum part of the formula for a signals power

A more general expression states that for $ M \geq N$:

$$ \sum_{n= N}^{n = M} c =  (M-N+1) \cdot c $$

where the derivation simply relies on fact that the epxression has (M-N+1) terms :

$$ \sum_{n= N}^{n = M} c = \{ c + c + ... + c\} = (M-N+1) \cdot c $$

And when applied for your particular case (with $N = -M$) it becomes:
$$ \sum_{n= -M}^{n = M} c = (M-(-M)+1) \cdot c = (2M+1) \cdot c $$


Matlab normalizes in the inverse Fourier transform, so if you want to calculate the fft accurately, you have to normalize by the fft vector by the lenght of the fft itself.

This matlab script will demonstrate the normalization. There is a 1 seconds long 1000 Hz
pure sine wave with amplitude 1. We will sample it with fs = 6000 Hz, and calculate the spectrum to it in 60 points.

After plotting, we expect two pure spikes at 1000 and -1000 Hz with amplitude 0.5.

fs = 6e3;
T = 1;
f = 1000;
x = sin(2*pi*f*(0:1/fs:T-1/fs));
Nfft = 60;
X = fft(x,Nfft)./Nfft;
stem((0:Nfft-1)*fs/Nfft-fs/2,abs(fftshift(X)))
title('Normalized spectrum')
xlabel('frequency')


Slight crackle with overlap-add frequency-domain convolution

For FFT fast convolution, the complex multiplication in the frequency domain needs to be done over the full length of the FFT (fftSize), not just the length of the impulse response in the time domain.  Otherwise (for strictly real data and a strictly real impulse response) the IFFT input will likely no longer be conjugate symmetric.


Generate and Recognize Graphical Patterns

I think deep learning comes very handy to solve such problems in these days. Simply create a convolutional neural network (CNN) architecture. Your input are the images, which are resized to one common dimension. The output is the class (type) of the image. For a CNN, this problem is rather simple and if you don't have a lot of training examples, a simple architecture such as LeNet would suffice. Of course, if you have a lot of classes, maybe a larger network (such as AlexNet) might be more appropriate. In either case, you could benefit from a large set of libraries, which are mostly coded in Python. To name a few: Tiny-DNN, Caffe, Tensorflow and Torch.


How does OpenCV find chessboard corners?

The source code of OpenCV is available, so I'd recommend just taking some time going through the code. The relevant file for this particular function is:

calib3d/src/calibinit.cpp

I've not looked into it in detail, but it looks like

CV_IMPL
int cvFindChessboardCorners( const void* arr, CvSize pattern_size,
                             CvPoint2D32f* out_corners, int* out_corner_count,
                             int flags )

is the main implementation of this method. In here they

Use cvCheckChessboard to determine if a chessboard is in the image
Convert to binary (B&W) and dilate to split the corners apart
Use icvGenerateQuads to find the squares.

The code then seems to go though a set of checks to condense these quads to chessboard corners, including icvFindConnectedQuads, icvCleanFoundConnectedQuads to remove extra corners, icvCheckQuadGroup, and icvCheckBoardMonotony.

All of these functions are implemented in the same file, apart from cvCheckChessboard which is in calib3d/src/checkchessboard.cpp. Depending how well you wanted to understand the code, there seem to be a number of debugging lines, which can be included if you #define DEBUG_CHESSBOARD, that may help you see what's going on.


How can I construct a Band-pass filter from a low and a high-pass filter?

Yes you are correct. You apply them if series in they are linear.

One simple band-pass filter you could use is called the difference of Gaussian (DoG)

The procedure is:

Create a Gaussian filter with a small variance
Create a Gaussian filter with a large variance
Subtract the latter from the former to create a band-pass filter
Apply the filter to the image

There are a lot of different filters to choose from depending on the application. A good paper is "On the choice of band-pass quadrature filters" by Boukerroui, Noble and Brady. It has an analysis of the DoG filter as well.


the demonstration that states that the FT of the ACF function is the square of the DTFT of the signal

You are right, the derivation is full of typos. The first equation below Eq. $(8.39)$ should read

$$\int_{-\infty}^{\infty}x(t+\tau)e^{\color{red}{-}j\omega\tau}d\tau=X(\omega)e^{j\omega \color{red}{t}}\tag{1}$$

Substituting into $(8.39)$ gives

$$\begin{align}\mathcal{F}\big\{R(\tau)\big\}&=\int_{-\infty}^{\infty}x(t)X(\omega)e^{j\omega t}dt\\&=X(\omega)X(-\omega)=|X(\omega)|^2\end{align}\tag{2}$$

where the last equality is only true for real-valued $x(t)$. However, the overall result is also true for complex-valued $x(t)$ because in that case the ACF is defined differently:

$$\mathcal{F}\big\{R(\tau)\big\}=\int_{-\infty}^{\infty}x^*(t)x(t+\tau)dt\tag{3}$$


voice recognition vs speech recognition MFCC

This is the exact same extraction procedure. Parameters such as frame rate, number of filter banks, band spacing etc. might vary a bit from one system to another, but this is just a consequence of the tuning procedure (the parameters yielding the best performance are selected), more than a fundamental difference. The "sweet spot" for some parameters such as frame duration is the same for both tasks.

You have to keep in mind that MFCCs capture an overall timbre parameter which measures both what is said (phones) and the specifics of the speaker voice & gender, and it is hard to factor these two dimensions apart. The MFCC vector of Alice saying "o", Bob saying "o", Alice saying "a" and Bob saying "a" are all different - which would make it look that this feature would be hard to use for speech recognition (due to speaker variability) and for voice identification (due to phone variability). There are workarounds, though...

In voice recognition systems, you have to make sure that the training data contains enough occurrences of all phones for each speaker. If you train your Bob voice model on Bob's "o"'s, it won't work when you feed it Bob's "a"'s. Because of that, the training set of MFCC data for a given speaker will have to be by design quite heterogenous - and this is why gaussian mixture models have been so successful at this task.

The reverse problem appears in speech recognition system: if you train a recognition system trained on Bob's voice, and try to feed it with MFCC data computed on Alice's voice, it will not work. One solution to this problem (vocal tract length normalization) is to find a linear transform of the MFCC vector such that, when applied to Alice's MFCC data, make it similar to what the model has captured.


Depends on context.

For strictly real-valued signals it sometimes refers to creating an analytic signal (see: http://en.m.wikipedia.org/wiki/Analytic_signal ) from a signal with no imaginary component, and taking the magnitude (or complex with both magnitude and phase) envelope of that hypothetical now complex creation.


What is the meaning of $l_p$ norm in this model for sparse channel estimation?

The symbol ${||h||}_p$ refers to the p-norm of the vector. It is calculated as

$${||h||}_p=\sqrt[p]{\left|h_1\right|^p+\left|h_2\right|^p+...+\left|h_n\right|^p}$$

When they write ${||h||}^p_p$, it is just that multiplied by itself $p$ times, so that the root disappears. Mathematically:

$${||h||}^p_p=\sum\limits_{i=1}^n \left|h_i\right|^p$$


Odd-order Butterworth filter gain at cutoff frequency

Congratulations, you've found an error on the wikipedia page on Butterworth filters. The squared magnitude of the frequency response of an $n^{th}$-order Butterworth low pass filter is

$$|H(j\omega)|^2=\frac{1}{1+\left(\frac{\omega}{\omega_c}\right)^{2n}}\tag{1}$$

There should be no imaginary unit $j$ in that formula.


No problem at all if $B_s \ll B_c$.
However, $B_s = 10\text{Hz}$ is considerably small that makes synchronization very difficult especially for low-cost receivers.
For example, Sigfox Ultra Narrow Band uses $10\text{Hz}$ bandwidth (are you studying the standard?) that requires precise but expensive oscillators like TCXO. Even with that, there must be random access enhancement to deal with the frequency imprecision.
There is also some effort to design transceivers without TCXO for such narrow bandwidth. Check e.g. Foxy.


If you are using Matlab, the plots looks like Matlab, you can use indexing directly to get a threshold value:

% 'time' is the time series
% 'signal' is the signal itself
idx = signal > 0.2;
plot(time(idx), signal(idx));

You are not clear on what you're trying to achieve, but there could be some information found in the frequency spectrum, so that you could do for example a fourier transform on the signal.


Dimension in mathematics is defined as the number of independent components of a structure (with some generalizations, like fractal dimensions). Now which number of dimensions you assign to e.g a discrete image (as in picture) depends on what aspect of the image you are describing. As a vector space where each pixel is a component, you get as many dimensions as pixels. The same for a discrete signal, the number of samples is your dimension of the vector space.

However, sometimes you're rather interested in the manifold the object is a function on. The real line is a 1-dimensional real manifold, and a real function on it would be a signal. So you can call the signal "one dimensional", explicitly referring to this property. Sometimes your function may map to the complex numbers or to R^2 (still from the real line!) instead, and you might want to call that a "two dimensional signal", and that's also correct. However, a two dimensional signal may also just be a function on a two dimensional manifold, like a plane.

So there's a lot ambiguity in verbal descriptions, which is why to be clear you should use exact mathematical notation. Opposite to common belief it was not just invented to torture students.


What is the Sound of a Solution to the Wave Equation?

Since you already have the solution to the wave equation, what you need to do now is define a "pick-up" point somewhere along the length of the string (some $x$) and obtain its trace in time, i.e. look at what that specific $x$ value is "doing" across all the $t$ instances that you have. In other words, draw a straight line down your vibrating string diagram and take the vertical displacement of your string at each time point as the amplitude of your sound signal at that time point.

In reality, when you hear an electric instrument you hear exactly that, you hear what the vibrating string is doing exactly above the pick-up coil. Furthermore, electric guitars (and other instruments) have more than one pick-up points along the length of a string and a switch that is used to change the way these pick-up points are combined to offer different qualities of sound. For example, the coil closest to the neck (closest to the middle of the distance between the neck and bridge) produces sound with more low frequencies, while the coil closest to the bridge (the fixed point at the lower end of the instrument) produces sound with more high frequencies. The placement of the coils is a very important issue. Imagine for example what would happen if you placed the coil at a nodal point on the string (when it is excited to produce a standing wave).

Similarly, in acoustic instruments, the sound that reaches the ear is the projection of the vibrating instrument (notice now we are talking about the whole instrument) towards the direction of the ear or microphone. So, in effect, you are still applying this 'slice of vibrating element in time' concept. By the way, this is why mic placement is important when you record acoustic instruments (see for example: http://www.soundonsound.com/sos/1995_articles/sep95/mikinginstruments.html)

Hope this helps.


Can someone explain what a butterfly is in regards to DSP?

Wikipedia has this one covered:

  In the context of fast Fourier transform algorithms, a butterfly is a portion of the computation that combines the results of smaller discrete Fourier transforms (DFTs) into a larger DFT, or vice versa (breaking a larger DFT up into subtransforms). The name "butterfly" comes from the shape of the data-flow diagram in the radix-2 case, as described below.1 The earliest occurrence in print of the term is thought to be in a 1969 MIT technical report.[2][3] The same structure can also be found in the Viterbi algorithm, used for finding the most likely sequence of hidden states.

(Source)


How to find PSD (Power spectral density) of spatial data

A spectrum estimation method that works well in many cases is Welch's method. The MATLAB implementation is pretty straight forward to use, even when not knowing about the math behind it. Following Marcus' recommendations, you could try:

if norm(diff(X,2)) < sqrt(eps) % is X evenly spaced ?
    Fs = 1/mean(diff(X)); % find the step size of X
    [Pyy,F] = pwelch(Y,[],[],[],Fs,'onesided','PSD'); % compute the PSD
    plot(log10(F),log10(Pyy)), xlabel('log [\mum^{-1}]'), ylabel('log[PSD[\mum^2/\mum^{-1}]]')
    axis tight, grid on
end


Precise Centre frequency of an All-pole digital filter

For a filter consisting of a complex conjugate pair of poles, the $z$-domain transfer function is:
$$H(z) = \frac{a}{\left(1-r(\cos{\theta}-i\sin{\theta})z^{-1}\right)\left(1-r(\cos{\theta}+i\sin{\theta})z^{-1}\right)}\\
= \frac{a}{1 - 2r\cos(\theta)z^{-1} + r^2z^{-2}},$$

where $a$ is a constant by which the magnitude frequency response can be normalized, $r$ is the pole radius, and $\theta$ is the pole angle. The peak of the magnitude frequency response will not be exactly at frequency $\theta$, but at a nearby frequency $\omega_0$ at which the derivative

$$\frac{|a|2r\sin(\omega)\left(2r\cos\omega - (r^2 + 1)\cos\theta\right)}{\left(4r^2\cos^2\theta - 4r(r^2 + 1)\cos(\theta)\cos(\omega) + 4r^2\cos^2\omega + r^4 - 2r^2 + 1\right)^{3/2}}$$

of the magnitude frequency response $|H(e^{i\omega})|$ with respect to $\omega$ equals 0. That frequency is:

$$ω_0 = \operatorname{acos}\left(\frac{(r^2 + 1)\cos\theta}{2r}\right)\text{ if }-1 < (r^2 + 1)\cos(\theta)/(2r) < 1$$

Figure 1. The condition for that the magnitude frequency response peaks elsewhere than at $w_0=0$ or $w_0=\pi$ is that poles stay outside these droplets.

In the following it is assumed that the condition is satisfied as otherwise the peak would be trivially at frequency $0$ (if $\theta < \frac{\pi}{2}$) or $\pi$, whichever is closer to the pole. The normalized magnitude frequency response with $a$ chosen so that $|H(e^{i\omega_0})| = 1$ reads:

$$|H(e^{i\omega})| = \sqrt{\frac{(r^4 - 2r^2 + 1)\sin^2\theta}{4r^2\cos^2\theta - 4r(r^2 + 1)\cos(\theta)\cos(\omega) + 4r^2\cos^2\omega + r^4 - 2r^2 + 1}}.$$

Given $\omega_0$ and $r$:

$$\theta = \operatorname{acos}\left(\frac{2r\cos\omega_0}{r^2 + 1}\right)\tag{1}$$

which is something you can use directly to keep the peak of your all-pole bandpass filter where you want it:

$\omega_0$ pole loci">
Figure 2. Pole loci for a selection of constant $\omega_0$ traced over $0 < r < 1$.

Applying the substitution given by equation (1) gives:

$$H(z) = \frac{a}{1 - \frac{4r^2\cos(\omega_0)}{r^2 + 1}z^{-1} + r^2z^{-2}}$$

Just like $\theta$ does not exactly give the peak frequency, $r$ does not give the exact bandwidth. The exact -3 dB to -3 dB bandwidth $B_\text{-3 dB}$:

Figure 3. -3 dB bandwidth of the filter, with $\omega_0 = \frac{\pi}{4}$ and $r=0.9$.

is given by:

$$B_\text{-3 dB} = \operatorname{acos}\left(\frac{(r + 1/r)\cos\theta + (r - 1/r)\sin\theta}{2}\right)\\- \operatorname{acos}\left(\frac{(r + 1/r)\cos\theta + (1/r - r)\sin\theta}{2}\right)$$

This was found by solving the values of $\omega$ for which $|H(e^{i\omega})| = \sqrt{2}/2$. The first $\operatorname{acos}$ gives the solution found to the right from the peak and the latter $\operatorname{acos}$ the one to the left. More generally, for a bandwidth from magnitude frequency response value $g$ to $g$:

$$B_g = \operatorname{acos}\left(\frac{(r + 1/r)\cos\theta + \sqrt{1/g^2 - 1}(r - 1/r)\sin\theta}{2}\right)\\-\operatorname{acos}\left(\frac{(r + 1/r)\cos\theta + \sqrt{1/g^2 - 1}(1/r - r)\sin\theta}{2}\right)$$

This is needed when designing a composite filter, for example a gammatone filter, that is a cascade of identical two-pole filters. Usually we want to design the filter to have a given bandwidth, but it is very difficult to solve $r$ from the above. So we resort to approximation. Using the value $D_2$ of the second derivative of the magnitude frequency response at $\omega_0$:

$$D_2 = \frac{(r^2 + 1)^2\cos^2\theta - 4r^2}{(r^2 - 1)^2(1 - \cos^2\theta)},\tag{2}$$

we can approximate the magnitude frequency response with that of a complex one-pole filter that has at $\omega_0$ the magnitude frequency response and its first two derivatives matched to those of the two-pole filter:

Figure 4. One-pole approximation with pole angle $\omega_0$.

The magnitude frequency response of the one-pole filter is:

$$|\hat H(e^{i\omega})| = \frac{\sqrt{2}\left(\sqrt{1 - 4D_2} - 1\right)}{2\sqrt{\left(2D_2\cos(\omega - \omega_0) - 2D_2 + 1\right)\left(1 - \sqrt{1 - 4D_2} - 2D_2\right)}}.$$

It is symmetrical around $\omega_0$ and has $g$ to $g$ bandwidth:

$$\hat B_g = 2\operatorname{acos}\left(\frac{1/g^2 - 1}{2D_2} + 1\right),$$

or inversely and more usefully:

$$D_2 = \frac{1 - g}{2g^2\left(\cos(\hat B_g/2) - 1\right)}.$$

Equations $(1)$ and $(2)$ when combined give:

$$D_2 = \frac{2r^2(r^2 + 1)^2\left(1 - \cos(2\omega_0)\right)}{(r^2 - 1)^2(2r^2\cos(2\omega_0) - r^4 - 1)}.$$

Solving that for $r$ gives a cumbersome equation with $\omega_0$ and $D_2$ as arguments that can be used to set the approximate bandwidth of the two-pole filter (c language compatible plaintext with w_0 = $\omega_0$ and d_2 = $D_2$):

r=sqrt(-sqrt(2)*sqrt(sin(w_0))*sqrt((-sin(w_0)*(pow(d_2+1,2)*pow(cos(w_0),2)
+d_2-1)*sqrt(-pow(d_2+1,2)*pow(cos(w_0),2)+pow(d_2,2)-2*d_2+1)+pow(d_2+1,3)
*pow(cos(w_0),4)-(d_2+1)*(pow(d_2,2)-d_2+2)*pow(cos(w_0),2)+pow(d_2,2)-2*d_2
+1)/sqrt(-pow(d_2+1,2)*pow(cos(w_0),2)+pow(d_2,2)-2*d_2+1))+sin(w_0)*sqrt(-
pow(d_2+1,2)*pow(cos(w_0),2)+pow(d_2,2)-2*d_2+1)-(d_2+1)*pow(cos(w_0),2)+1)
/sqrt(-d_2)

or with alternative symbols:

r = √(- √2·√(SIN(ω_0))·√((- SIN(ω_0)·((d_2 + 1)^2·COS(ω_0)^2 + d_2 - 1)·√(-
(d_2 + 1)^2·COS(ω_0)^2 + d_2^2 - 2·d_2 + 1) + (d_2 + 1)^3·COS(ω_0)^4 - (d_2
+ 1)·(d_2^2 - d_2 + 2)·COS(ω_0)^2 + d_2^2 - 2·d_2 + 1)/√(- (d_2 + 1)^2
·COS(ω_0)^2 + d_2^2 - 2·d_2 + 1)) + SIN(ω_0)·√(- (d_2 + 1)^2·COS(ω_0)^2 +
d_2^2 - 2·d_2 + 1) - (d_2 + 1)·COS(ω_0)^2 + 1)/√(-d_2)

When $\hat B_\text{-3 dB} = \hat B_{\sqrt{2}/2}$ is set proportional to $\omega_0$ by a constant $c$ stepped from 0.1 to 1.0, the upper half-plane pole resides on tracks like these:

Figure 5. Pole loci for bandwidth proportional to peak frequency with proportionality constant $c$.


Calculate cut-off frequency of lowpass IIR-filter

This will work but only if the it's actually a low pass filter. Your formula describes any possible IIR filter so it could really be anything.


Extraction of non-sinusiodal repetition rates

You are right that the repetition is around 650 by how exactly do I compute that automatically? Seems like a peak-picking problem to me? Or is there some other methods that can be used?

Yes, it's just peak-picking.  Your period is the x value of the first strong peak:

Your peaks are all similar in height, probably because you're doing the autocorrelation using the FFT?  So it's a circular autocorrelation.  You can

make it non-circular by zero-padding before doing the FFT, or
could skew the plot you currently have by adding a linear function that emphasizes the peaks close to 0 and reduce the height of the farther ones.

It depends on what you're specifically looking for.  In either case you then have to pick the highest peak that's not at 0 lag.

I adapted my script from here and it (surprisingly) worked without any tweaking:

raw1 x = 709.23

raw2 x = 710.77

from pylab import *

# Import matlab files
import scipy.io
#sig = scipy.io.loadmat('raw1.mat')['raw1'][0]
sig = scipy.io.loadmat('raw2.mat')['raw2'][0]

# Calculate autocorrelation (same thing as convolution, but with
# one input reversed in time)
corr = fftconvolve(sig, sig[::-1], mode='full')

# throw away the negative lags, so x = 0 means 0 lag
corr = corr[len(corr)/2:]

# Find the first low point from derivative
d = diff(corr)
start = find(d > 0)[0]

# Starting at the first low point, find the highest peak
peak = argmax(corr[start:]) + start

# Fit parabola to estimate more precise location of peak
px, py = parabolic(corr, peak)

plot(corr)
axvline(px, color='r')

print px

The variable peak is the integer x value at the peak (709).  Then the parabolic function fits a parabola to the peak and tries to get a more precise estimate (709.232...), which may or may not be relevant to what you're doing.

The python fftconvolve() function does zero-padding internally to produce a linear convolution, like convolve() does


I do not think that there is a standard approach to this problem, because as far as I know it has not yet been solved satisfactorily. I believe that you can find a solution to a (very) simplified version of the problem by restricting it for example to the recognition of a fixed set of animal sounds. The way to go would probably be to use methods from speech recognition (Hidden Markov Models, HMM), because the problem is similar to voice control systems, where commands from a fixed set of predefined commands are identified. With some good engineering, these systems can work relatively well.

In the case of animal sounds you would need to define a fixed set of animal sounds that you want to be able to recognize, collect a large database of such sounds, and use the data to train HMM models, which can then be used to score unknown sounds. Of course, such a system does not even come close to human capability of distinguishing sounds, but for such a restricted application it may work well enough.

So in that sense I do not agree that voice control systems are "very inaccurate" if there is no context available, because that's the usual scenario for simple voice control systems, which can in fact be made to work in practice.


Relation between bandwidth and data rate

A pure never changing sinusoidal carrier of a given frequency carries no information.  Changes or modulations to that carrier over time (as in: per second) can carry information.  The maximum rate of information is related to the bandwidth and S/N ratio, not the number of unmodulated carrier cycles.

Of course, it may be easier to implement certain modulation schemes with a higher carrier to bandwidth ratio.


Is there a common name for the first order IIR averaging filter?

This is often called a leaky integrator, a special case of a first-order IIR lowpass filter. They are discussed in more detail in several previous questions:

Is a leaky integrator the same thing as a low pass filter?
Is there a technical term for this simple method of smoothing out a signal?
Designing a practical integrator
"Low-pass filter" in non-EE, software API contexts


The first thing you must change for understanding why IFFT is used for OFDM is to forget about time and frequency domains.

As the name states OFDM is a technique about frequency domain multiplexing. When you think in FDM you'll have several channels modulated with a guard interval between then to prevent interference. If the carriers are orthogonal you'll not have interference, but without the waste of band guard.

But how you can have this?

One way is to use several modulators, as in FDM, but with the correct frequencies. This is hard to accomplish.

If you will implement digitally the filter bank,and this is the correct way to view IFFT on OFDM systems: a transmultiplexer(TMUX), is equivalent to an IDFT. This is why you can use IFFT as a fast implementation for OFDM. You can read one of the books about multi-carrier communication for a better explanation. This paper
http://wsl.stanford.edu/~ee359/ofdm_bingham.pdf is also one that I like about multi-carrier modulation.


A cubic interpolation function: folkloric copypasta or clever trade-off?

A reasonable, analytically derived cubic interpolation method will give a straight line if the inputs are on a straight line. With $a = 0$, $b = 1$, $c = 2$, $d = 3$ you get from cubic_bourke $A = d - c - a + b = 3 - 2 - 0 + 1 = 2$ as the coefficient for $x^3$ whereas for a straight line you would have $A = 0$ (and $B = 0$ as the coefficient for $x^2$). So the method doesn't give a straight line for straight-line input. Instead it will, as you say, wobble. In the frequency response that will manifest itself as a low multiplicity of the roots at multiples of the sampling frequency.
The impulse response is not symmetrical. Normally splines have a symmetrical impulse response. The asymmetry will result in phase non-linearity:
Figure 1. Impulse response of cubic_bourke.
The magnitude frequency response shows a rather nice treble response for audio purposes, but gives a lot of aliasing, worse than linear interpolation:
Figure 2. Magnitude frequency rseponse of cubic_bourke (red), linear interpolation (blue) and cubic Hermite spline (green). Half the input sampling frequency is at $\omega = \pi$.
Cubic Hermite (Catmull-Rom) will give better treble response up to 15 kHz for a 44.1 kHz input sampling frequency, and lower aliasing, than cubic_bourke.
So I say scrap it.


How to demodulate an AFSK signal in software

I think that you could get the best performance in terms of demodulator bit-error rate (BER) with a phase-locked loop.  You need it to be fast, though.  I think your best bet for a fast algorithm that still performs reasonably well is zero crossing.

On a side note, I would like to suggest that you change the 2200 Hz to 2400 Hz.  A naive implementation of the 1200/2200 Hz scheme would yield discontinuities, as seen about two thirds into the plot below, where the 2200 Hz transitions to 1200 Hz.

In order to minimize the bandwidth that you're using and avoid discontinuities that will distort the signal you'll need to make the phase continuous.  Even if you make the transmitter phase continuous, though, there will still be the issue that the 2200 Hz symbols will not always have the same number of zero crossings due to the different phases.   Usually they will have four zero crossings, but sometimes they will have three.  The 1200 Hz symbols on the other hand, will always have two zero crossings because the baud rate divides evenly into the FSK frequency.

You can solve both of these problems by changing the 2200 Hz to 2400 Hz.  Then the symbols will always start and end at 0 degrees (thus automatically making them phase continuous), and they will always have the same number of zero crossings- two and four.


Short answer

Normally a chromagram constits of chroma feature vectors per time frame. If you want to have a musical time axis in beats instead of a physical time axis in seconds, you have to get beat information (e.g. by a beat tracker) and average all chroma features that belong to the time span of a beat.

Long answer

Since the comments asked for a longer answer, here it is: A chromagram presents the energy that is contained in the twelve pitch classes over time. Starting from the STFT, you could compute a chromagram version by summing up all frequency bins that correspond to a particular pitch class, individually for each time frame. Since the STFT uses a constant hopsize, you end up with a chromagram which is constant in physical time. Since chromagrams are used for music applications, it is sometimes desirable to have a musical time axis, e.g. in beats. Since the tempo of music can change, the physical time span of an beat can vary. For converting physical to musical time, you could annotate the beat positions within the music or apply an automated procedure, kown as beat tracking. Then you could simply average the chroma features that belong to the time span of a beat.

In more technical wordings: $\mathcal{C} \in \mathbb{R}^{N \times 12}$ is your normal chromagram matrix with $N \in \mathbb{N}$ time frames. $\mathcal{B} : \mathbb{N} \rightarrow \{s| s \subset \{1,\dots,N\}\}$ is a function that gets the index of a beat as input and outputs a set of all time frame indices belonging to this beat index. You want to compute the beat-synchronous chromagram $\mathcal{C}^\text{beat} \in \mathbb{R}^{B \times 12}$ for $B \in \mathbb{N}$ beats. Usually $N \gg B$. You can do this simply like this: $\mathcal{C}^\text{beat}_{b, c} = \frac{1}{|\mathcal{B}(b)|} \sum_{n \in \mathcal{B}(b)}  \mathcal{C}_{n, c}$ with beat index $b \in \{1,\dots,B\}$ and pitch class index $c \in  \{1,\dots,12\}$.

This is not the place for comprehensive literature, but only two references for general chromagrams and beat-synchronous ones:

Meinard Müller: Fundamentals of Music Processing, Springer 2015. (Chapter 3.1)

Daniel P. W. Ellis and Graham E. Poliner: Identifying Cover Songs with Chroma Features and Dynamic Programming Beat Tracking, IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Honolulu, 2007, pp. IV-1429–IV-1432.


Find the length of the impulse response of a Linear Phase Type 4 FIR filter

You are missing a couple of zeros.

First of all, you must also include the reciprocal of the one located at $z=-2$, that is one at $z=-0.5$
You should also have one more zero coming from the fact that types 3 and 4 are anti-symmetric. This zero must be located at $z=1$, and is responsible for the minus sign in the anti-mirror image polynomial equation:

$H(z)=-z^{-L} H(\frac{1}{z})$

Given that you now have an odd number of zeros, you don't need to add any other one. If it were not the case, if you had an even number of zeros (type 3 filter), then you should add one in $z=-1$.


Scalar Value from the output of each mel filter bank is accumulated scaled magnitudes of frequencies present in mel filters range/bandwidth.


A filter with low phase shift or no phase shift at all? Which one?

If your filter is $h[n] \leftrightarrow H(e^{j\omega})$, then the DTFT of the time reversed filter $h[-n]$ is $H^*(e^{j\omega})$. Consider the result of passing your signal through the combined filter $h[n]*h[-n]$, which has a DTFT of $H(e^{j\omega})H^*(e^{j\omega})=|H(e^{j\omega})|^2$ and has zero-phase. In this setup, I've assumed that $h[n]$ is causal so that makes $h[-n]$ anticausal, and this makes the overall filter $h[n]*h[-n]$ noncausal and not applicable for real time applications but can be used for offline processing.

This is the idea used in the MATLAB command filtfilt. You pass filtfilt the filter coefficients and it runs the input signal through the filter and then the time reversed filter. Of course, you can also just as well do that by hand.

Here is a dumb example of a noisy sine wave and a moving average filter just to illustrate.


I think it should be:

$$ \boldsymbol{y}_{i} = {H}_{i} \boldsymbol{x} + \boldsymbol{n}_{i} $$

For the $ i $ -th antenna and $ {H}_{i} $ being the convolution matrix of this specific channel.
You could write this in a Matrix Form where $ H $ becomes a tensor.
In case $ \forall i, j, \; {H}_{i} = {H}_{j} $ then it can be made simpler.

You could also work on the samples:

$$ \boldsymbol{y} \left[ m \right] = H \left[ m \right] \boldsymbol{x} \left[ m \right] + n $$

Where in this case $ {y}_{i} \left[ n \right] $ is the $ i $ -th antenna at time sample $ m $ and the $ i $ -th row of $ H \left[ m \right] $ is the response of the channel $ i $ at time $ m $.
This model allows time variant channels.


Opinion about this aproach on Machine Learning

This is a fine idea and pretty standard. However, your problem will be the augmentation of the data. If you render the image from many viewpoints/rotations etc. and just use a plain background, the network would include that information into the training stage e.g. Like a shape cue. So, you have to use different backgrounds. This is of course scene specific. However, you might be able to narrow down your scenes so that you could randomly generate/perturb different background structures.

Also note that, you might be able to scale up to thousands of classes, but getting to millions could become a headache, since training this network would be very tough.

You might as well go with standard indexing / hashing approaches, which already perform pretty well.


Opencv color segmentation of traffic signs: HSV ranges

To those looking at this in the future, I realized I was being silly and forgot how HSV worked! To pick up the yellow more so than the green I needed to heighten the threshold of the saturation value within the range!

As a side note, I also was able to "normalize" the colors of the image a bit more by converting to Lab color space and applying CLAHE (for localized histogram equalization) and as well on the saturation values within the HSV color space. This allowed me to get the signs in ideal lighting conditions almost 100% of the time (roughly speaking of course) with very little false positives, but didn't really alleviate my problems on images with a redish/orangeish hue as seen in the second example image. But, as was mentioned by harshkn, this approach is very simple so I suppose I can only ask for so much success! Thanks to all who replied with suggestions!


Why do I have to divide by $\pi$ in order to scale amplitude of sine frequency spectrum?

The FFT should be divided by the number of elements in the vector.  You specified the vector length to be duration = int(2*pi*fs) which is a really weird way to defining the length, btw.  Thus, you need to divide by int(2*pi*fs).  Yes, you need to divide by a 2 as well, because the signal energy is distributed evenly between positive and negative frequencies.


Data smoothing in temperature measurements from MLX90614 BCC IR sensor

In any measurement application the VERY FIRST thing you should to do is to calibrate the hardware and characterize it's noise floor and data variance by doing measurements on a known device.
Start with measuring something that's guaranteed to be at a constant and known temperature and take histograms. These SHOULD have a more or less a gaussian shape with the mean matching the temperature of the device under test and the standard deviation being the capability of your measurement setup. If that's not the case or if the standard deviation is larger than you need it to be, you need to start debugging the hardware.
You histograms do not make a lot of sense. Typically you would expect a gaussian shape with a continuous distribution around the mean. The large gaps you have in your histograms indicate that there is something not working as it should .


Efficient calculation of exact frequency from a broad frequency range (1 Hz - 1000 Hz)

Thanks for citing my blog article on frequency calculation.  Based on your emphasis on efficiency I am going to recommend a different approach which should get you equivalent or superior results with way fewer calculations.

The first thing you should do is smooth your signal heavily with exponential smoothing.  This will squash your harmonics and mitigate noise.  You can carry your results forward buffer by buffer so it is extremely efficient.  While you are doing the smoothing, you can also check for zero crossings.  The signal need not be a "nearly perfect sine curve" for zero crossings to give you a good estimate of the frequency.  After heavy smoothing, the zero crossing should correspond to your fundamental frequency.

If your signal is sufficiently sine like, you can calculate the frequency quite accurately by using the three point (simplest case) formula with the appropriate spacing that you can find here:  Exact Near Instantaneous Frequency Formulas Best at Peaks (Part 1).  You will want to select the center point at the peak (or trough) and the other two points at the nearest zero crossings.  This formula will give you as good results as if you did interpolation at the zero crossings and calculated the frequency from there.  You can then calculate the frequency at every peak and trough, that is every half cycle, and generate a new sequence of values.  You can then apply a best polynomial fit on intervals of this sequence to give you very good calculations of the frequency and how it is varying.

A more accurate, but slightly more computationally expensive approach would be to calculate two bins of the DFT and use the formula found in my article: Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT.  For the best result the frequency should fall halfway between the bins so you want to select a frame interval of a whole number of cycles and a half.  Two and a half would be a good compromise between getting an accurate reading with a varying frequency and getting a better reading with multiple cycles.  Point spacing is key for saving calculations.  Four points per cycle will give the best results, so you want to be close to that.  For two and a half cycles, that means ten points, which is not many.  Since you are calculating the bin values separately, there is no advantage to a power of two.  You can build lookup tables for the ten points for the two bin and three bin ahead of time to make it even more efficient.  Pick a uniform spacing parameter that gets you as close to two and a half cycles as you can.  Your calculations result in a cycles per frame value, then use your spacing value to convert to cycles per whatever your measure is.

Feel free to add any followup questions in your question and I will respond with an edit on my own.

I'm glad to hear you find my blog articles useful.

"How do you "heavy smooth" second harmonic when the base frequency could be 100x higher? "

Exponential smoothing has the property that the higher the frequency the more the tone will be attenuated.  If you are still getting multiple zero crossing then you still have higher harmonics present with significant amplitude.  Increase your smoothing parameter even more.  If you are working in floating point, you don't have to worry about destroying your fundamental frequency.  It will always be less dampened than the higher frequencies.

Since you have such a high sampling rate, you may also want to downsample your signal some as well.  There are many ways to do this, but since this is not an area of expertise for me, I will hold off on recommendations.  I can say the lowest calculation method is to just sample every nth point, but this does nothing to mitigate alias frequencies.  Another simple method is average intervals to get a new sequence.

The nice thing about picking a whole plus a half number of cycles for your frame is that the second harmonic (twice the fundamental frequency), and all even harmonics above it, are a whole (or very close to) number of cycles in the frame so they will be orthogonal to your basis vectors.  If there are a lot of harmonics present, you may want to increase the number of points above four per cycle to push the Nyquist frequency up for your frame length.

You can find my email address on my profile page.  Feel free to contact me there.

P.S.  Thanks to the upvoters.


How to set STFT-parameters to visualize woodpecker / Sawtooth Signal?

Your first problem is that you have your axes mixed up.

Your "17 Hz" picture shows this quite clearly.  The red bars are the impulses from the signal.  The bars should indicate an impulse (wide range of frequencies) but your axes indicate they are across time.

This is what it should look like:

Notice that the bars are along the frequency axis rather than time as in yours.

This is your woodpecker recording:

It contains bars, too.  But, they are not as wide band and are much closer (in intensity) to the background noise.

The bars repeat at just above 20Hz - you can see that the frequency is around 23Hz in your last plot (the red line plot below the blue oscilloscope style plot.)

Your overlap of 0.95 amounts to a time resolution of about  just about 40Hz.  You need a higher overlap of maybe 0.99 to separate the bars.

You could downsample the recording down to like 11025 Hz sampling rate.  That would get you better frequency resolution and make calculating the spectrogram faster.  You can do that because there's not really much of interest in your signal above about 4500Hz.


Segmenting and labeling chessboards through computer vision

I have once done something similar. You should give a try to OpenCV's checkerboard detection. There are also works on Sudoku square detection, which might give you a hint.

If these don't work for you, then I would suggest the following: Forget about detection the checkerboard. Initially, assume that every piece is in place and take the snapshot of the board (you can even manually place an initial board, doesn't matter). Then, whenever a piece is moved, grab a new image and apply a differencing of two images (this is somewhat like a background subtraction). Note that this way, you easily get the location of the moved piece. And since you know the piece and what was in the previous location beforehand, you get the complete information of the move.

Use a downward looking camera, just like you showed.


Understanding what channel gain is and why it is important side information

If you transmit signal $s(t)$ and receive $$r(t)=gs(t)+n(t),$$ where $n(t)$ is noise, then $g$ is the channel gain. Note that, when the signals are low-pass equivalents (that is, complex envelopes), then in general $g$ is complex.

The channel gain can also be written for a discrete channel. If you matched-filter $r(t)$ and sample the result at the symbol times, then you'll obtain a sequence of numbers $$r_i = g s_i +n_i,$$ where $s_i$ are the transmitted pulse amplitudes and again $g$ is the channel gain. Note that the channel gain has two main effects:

It reduces the SNR, because when $|g|<1$, the received signal loses power compared to the noise.
It changes the phase of the signal.

In other words, it scales and rotates the transmitted constellation. It is common to assume that the receiver knows $g$, because it can be learned by using a training sequence. If the receiver knows $g$, then it can compensate the constellation rotation. In the case of a wireless channel with diversity, the receiver can align and add coherently the different received signals (see maximal ratio combining).

It is less common for the transmitter to know the channel gain, but if it does, it can do beamforming and increase the SNR. Basically the idea is to transmit a signal that is "aligned" with the channel's phase $\angle g$.

You're correct that the BER (along with maybe other channel statistics) is useful side information. Wireless systems such as Wi-Fi include a control channel where the receiver can provide these statistics to the transmitter, so that it can select an appropriate modulation, code rate, and symbol rate.


How to obtain a meaningful group delay graph using matlab

Does the impulse response of the FIR filter you've shown need to be reversed along the sample axis (i.e. horizontal axis)?


Yes, this is correct. Other values would cause a DC offset.


Computing main lobe width of Hamming window in Matlab by plotting

The problem is the way you shift the transform of the rectangular window to compute the transform of the Hamming window. You can't just shift on the discrete frequency grid, but you have to evaluate the (continuously) shifted functions on the grid. If $W_r(\omega)$ is the transform of the rectangular window

$$W_r(\omega)=\frac{1}{N}e^{-j\omega (N-1)/2}\frac{\sin(N\omega /2)}{\sin(\omega/2)}\tag{1}$$

then the transform of the Hamming window can be written as

$$W_h(\omega)=\alpha W(\omega)-\frac{\beta}{2}\left[W_r(\omega-\omega_0)+W_r(\omega+\omega_0)\right]$$

with $\alpha=0.54$, $\beta=0.46$, and $\omega_0=2\pi/(N-1)$.

In Matlab/Octave you get:

w = -pi:0.01:pi;
N = 20;
w0 = 2*pi/(N-1);
Wr = 1/N*exp(-j*w*(N-1)/2).*sin(N*w/2)./sin(w/2);
W1 = 1/N*exp(-j*(w-w0)*(N-1)/2).*sin(N*(w-w0)/2)./sin((w-w0)/2);
W2 = 1/N*exp(-j*(w+w0)*(N-1)/2).*sin(N*(w+w0)/2)./sin((w+w0)/2);
Wh = .54*Wr -.23*(W1+W2);

plot(w/pi,20*log10(abs(Wr)),w/pi,20*log10(abs(Wh)/max(abs(Wh))))
axis([-1,1,-60,0]), grid

EDIT:

Another error I've mentioned in the comments is the sign of the shifted windows in your formula. The positive sign that you use is valid for a zero-phase Hamming window, i.e. a window that is centered at $n=0$. However, your rectangular window starts at $n=0$, so it's not a zero-phase window. Consequently, combining the two windows as if they were both zero-phase doesn't work. That's where the negative sign in my formula comes from (both windows start at $n=0$). If you change the sign in your code the result is much closer to a Hamming window than before, but it's still off due to the shift on the grid I've mentioned above. In the figure below you can see the difference (blue: exact Hamming; red: your solution with corrected sign):


Phase difference between signals sampled at different frequencies

SOLUTION

Bottom Line

$$(\theta_2-\theta_1) = 2\pi f(T_2-T_1)n -(\phi_2[n]-\phi_1[n]) \tag{1}$$

$f$: frequency in Hz of two tones of the same frequency and fixed phase offset

$(\theta_2-\theta_1)$: phase difference in radians of tones being sampled

$T_1$: period of sampling clock 1 with sampling rate $f_{s1}$ in seconds

$T_2$: period of sampling clock 2 with sampling rate $f_{s1}$ in seconds

$\phi_1[n]$: phase result from sampling tone with $f_{s1}$ in radians/sample

$\phi_2[n]$: phase result from sampling tone with $f_{s2}$ in radians/sample

This shows how any standard approach of finding the phase between two tones of the same frequency that are sampled with the same sampling rate (common phase detectors approaches including multiplication, correlation etc) can be extended to handle the case when the two sampling rates are different.

Simpler explanation first:

Consider the exponential frequency form of equation (1):

$$e^{j(\theta_2-\theta_1)} = e^{j2\pi f(T_2-T_1)n}e^{-j(\phi_2[n]-\phi_1[n])} \tag{2}$$

The term $e^{j2\pi f(T_2-T_1)n}$ is the predicted difference in frequency between the two tones that would result from sampling a single tone with two different sampling rates (when observing both on the same normalized frequency scale).

The observed difference in frequency between the two tones would be $e^{j(\phi_2[n]-\phi_1[n])} $.

Both terms have the same frequency with a fixed phase offset. This phase offset is to the actual difference in phase between the two continuous-time tones. By conjugate multiplication we subtract the two, removing the phase slope and the fixed phase difference results.

Derivation

The approach is to carefully work with all units with a time axis of samples. The frequency domain is thus in units of normalized frequency: cycles/sample or radians/sample corresponding to cycles/sec or radians/sec when the time axis is seconds. Therefore our sampling rate, regardless of what it is in time given in seconds, will be always equal to $1$ cycle/sample (or $2\pi$ radians/sample if working in normalized radian frequency).

The two signals with the same analog frequency once sampled each with a different rate in the time domain, will be two signals each with a different normalized frequency.

This simplifies the problem to gives us the following result:

Given our original signals as normalized sinusoidal tones at the same frequency with different phase offsets:

$$x_1(t) = \cos(2\pi f t + \theta_1) \tag{3}$$
$$x_1(t) = \cos(2\pi f t + \theta_2) \tag{4}$$

Once sampled but still with time in seconds:
$$x_1(nT_1) = \cos(2\pi f n T_1 + \theta_1) \tag{5} $$
$$x_2(nT_2) = \cos(2\pi f n T_2 + \theta_2) \tag{6}$$

Equation (5) and Equation (6) expressed time in units of samples is:

$$x_1[n] = \cos(2\pi f T_1 n+ \theta_1) \tag{7}$$
$$x_2[n] = \cos(2\pi f T_2 n+ \theta_2) \tag{8}$$

Convert to complex exponential form so that we can easily extract the phase terms using complex conjugate multiplication, (for a single tone we just need to split the input signal into quadrature components; $\cos(\phi) \rightarrow [\cos(\phi),\sin(\phi)]\rightarrow \cos(\phi)+j\sin(\phi) = e^{j\phi}$, this is described using the Hilbert Transform as $h\{\}$)

$$h\{x_1[n]\} =e^{-j(\phi_1[n])} =  e^{2\pi f T_1 n+ \theta_1} = e^{2\pi f T_1 n}e^{\theta_1} \tag{9}$$
$$h\{x_2[n]\} = e^{-j(\phi_2[n])} =e^{2\pi f T_2 n+ \theta_2} =e^{2\pi f T_2 n}e^{\theta_2} \tag{10}$$

The complex conjugate multiplication gives us the difference phase term we seek and its relation to our measured results:

$$e^{-j(\phi_2[n]-\phi_1[n])} = e^{2\pi f T_2 n}e^{\theta_2}e^{-2\pi f T_1 n}e^{-\theta_1} \tag{11}$$

Resulting in

$$e^{j(\theta_2-\theta_1)} = e^{j2\pi f(T_2-T_1)n}e^{-j(\phi_2[n]-\phi_1[n])} \tag{12}$$

Note that $e^{-j(\phi_2[n]-\phi_1[n])}$ represents the measurement which for single tones will result in a frequency and this frequency is predicted to be $\omega = 2\pi f(T_2-T_1)n$, given by the  $e^{j2\pi f(T_2-T_1)n}$ term. If we remove the frequency offset (by the multiplication above), the result is the phase difference of the original signal.

Taking the natural log of both sides reveals the result in units of phase (radians):

$$(\theta_2-\theta_1) = 2\pi f(T_2-T_1)n-(\phi_2[n]-\phi_1[n]) \tag{13}$$

So in summary, $\phi_1[n]$, $\phi_2[n]$ come from our measurements given as $cos(\phi_1[n])$, $cos(\phi_2[n])$ and we establish the difference that we need to get our answer through the complex conjugate multiplication of the Hilbert Transform of those measurements.

Demonstration

I demonstrate this with the script below similar to the OP's configuration with the results plotted below, which now includes the decimation and was tested for both f = 25 MHz and f = 400 MHz (undersampled) with similar results This shows each step to demonstrate the process above, and the operations can be further combined. The Hilbert Transform in implementation would be any approach of choice to delay the sampled tones 90° (A fractional delay all-pass filter is a reasonable choice).

Len = 10000;
phase_diff_in = 45;
f=400e6; % Incoming signal frequency
D = 13
Fs1 = 157e6*D;
Fs2 = 121e6*D;
t1 = [0:Len-1]/Fs1;  % Time vector channel 1
t2 = [0:Len-1]/Fs2;  % Time vector channel 2
phi1 = 2*pi*f*t1;
phi2 = 2*pi*f*t2 + deg2rad(phase_diff_in);
sign1 = cos(phi1);
sign2 = cos(phi2);

% emulation of perfect Hilbert Transform for each tone:
c1_in = 2*(sign1 - 0.5*exp(j*phi1));
c2_in = 2*(sign2 - 0.5*exp(j*phi2));

% create expected phase slope to remove
n = [0:Len-1];
comp_in = exp(-j*2*pi*f*(1/Fs2-1/Fs1)*n);

% decimation
c1 = c1_in(1:D:end);
c2 = c2_in(1:D:end);
comp = comp_in(1:D:end);
pdout = c1.*conj(c2);
result = pdout.*comp;

%determine phase_diff
phase_out = rad2deg(unwrap(angle(result)));
mean_phase = mean(phase_out);

Below is the result for two test cases, 0° as the OP was trying in his example and then a 45° phase shift.

Below shows the result for the copies of the input signal at frequency $f$ sampled by $f_{s1}$ as sig1 and $f_{s2}$ as sig2 for the case of zero degree phase between them. The real of the complex conjugate product pdout is the bold red sinusoid, and we note that it has zero phase offset.

To confirm the calculations, the plot below compares it directly to the real of the compensation term $cos(2\pi f(T_2-T_1)) to see that they are the same frequency consistent with the equation.

And repeating with $\theta_2-\theta_1 = 45°$

The result of the raw phase data for every sample shows that each sample individually has extremely low noise (limited by numerical precision, so the result can be determined with very few samples!). Such performance will depend on the actual quality of the Hilbert transform to accurately delay the input tone by 90° to create a qaudrature copy. Under conditions of noise the result can be averaged to the extent of waveform stability for a very robust solution.

Extended testing of performance with the undersampling case shows excellent results (f = 400e6):

Testing every difference angle in 1 degree steps:

RMS error of 10,000 samples (Note vertical axis is in increments of 0.5e-11)

Result of a greatly extended frequency sweep of the input frequency from 1e6 to 4000e6 in steps of 1e6 with a 45 degree phase shift with 10,000 points measured at each frequency showed a consistent result for phase determination at all frequencies (oversampling and undersampling). This is with the OP's configuration with the two frequencies including the decimate by 13. (Thus the sampling rate of each of the input tones after decimation for this test was at fs = 157e6 and 121e6, thus the far right of this graph with the frequency of the tone being sampled being 4e9 is significantly under-sampled. The RMS error is proportional to the frequency of the tone as shown, but even under this extreme condition, the error is still less than 5e-10 degrees. (8.7e-12 radians or -221 dB).

Practical Limitations

The accuracy of the above result is limited by knowledge of the exact frequencies and phase relationship given by $f_{s1}$ and $f_{s2}$, and knowledge of the frequency $f$ of the tone being sampled.

(As written the solution also assumes that the two sampling clocks both start at time $t=0$, but the sampling offset can be added starting with equation (8) with a similar result; bottom line is the starting phase relationship between the two sampling clocks must be known or measured as it will introduce an additional offset).

The reality is that no two free-running clocks will stay in perfect synchronization; there will be an inevitable drift in the actual frequency and phase difference between the sampling clocks that are not locked to a common reference (see Segal's Law https://en.wikipedia.org/wiki/Segal%27s_law). One of the clocks must be declared our reference of time (and our measurement will be to the accuracy of that one clock). If the clocks are not physically co-located, two-way time transfer techniques (see https://tf.nist.gov/time/twoway.htm) can be used to measure one clock versus the other. If they are physically co-located, then the simple thing to do would be to sample one clock with the other.

Below I show how this approach can completely eliminate one of the sampling clocks from the equation for our solution:  (I haven't yet tested this so may contain math errors)

Consider sampling $f_{s1}=\frac{1}{T_1}$ with $f_{s2}=\frac{1}{T_2}$. This will ultimately remove $f_{s2}$ from the equation entirely by using $f_{s1}$ as the common reference (we essentially measured $f_{s2}$ with $f_{s1}$ by sampling $f_{s1}$ with $f_{s2}$ allowing us to put the samples of $f_{s2}$ in units of $f_{s1}$ counts.):

$f_{s1}$ as a cosine:

$$x_{s1}(t) = cos(2\pi f_{s1}t) \tag{14}$$

When sampled with $f_{s2}$ given the constraint they both start at t=0 becomes:

$$x_{s_1}(nT_2) = cos(2\pi f_{s1}nT_2) =  cos(2\pi nT_2/T_1) \tag{15}$$

Which in units of samples is:

$$x_{s_1}[n] =cos(2\pi T_2/T_1 n) \tag{16}$$

Resulting in a third phase measurement in units of samples that we can get by sampling $f_{s1}$ with $f_{s2}$ (importantly to be done at the same time $x_1(t)$ and $x_2(t)$ are sampled!):

$$\phi_3[n] = 2\pi T_2/T_1 n \tag{17}$$

Thus if we don't know $T_2$ but have $\phi_3$ we can substitute the above  equation  to get:

$$T_2 = \frac{T_1 \phi_3[n]}{2\pi n} \tag{18}$$

substituting into (4):

$$
\phi_2[n]- \phi_1[n] = 2\pi nf\bigg(\frac{T_1 \phi_3[n]}{2\pi n}-T_1\bigg) + (\theta_2-\theta_1)  \tag{19}
$$

Resulting in the following solution for the desires phase difference of the original input signals:

$$
\theta_2-\theta_1=  2\pi f\bigg(\frac{T_1 \phi_3[n]}{2\pi n}-T_1\bigg)n - (\phi_2[n]- \phi_1[n])  \tag{20}

$$

Where

$f$: frequency of tone being sampled

$T_1$: period of sampling clock 1 with sampling rate $f_{s1}$

$\phi_1[n]$: result from sampling tone with $f_{s1}$, values will be $cos(\phi_1[n]) $

$\phi_2[n]$: result from sampling tone with $f_{s2}$, values will be $cos(\phi_2[n]) $

$\phi_3[n]$: result of sampling $f_{s1}$ with $f_{s2}$, values will be $cos(\phi_3[n]) $

Thus by only knowing $T_1$ which is $1/f_{s1}$, we can measure $f$ from the samples of $x_1(t)$ directly, measure $\phi_1[n]$ by sampling $x_1(t)$ with $f_{s1}$, measure $\phi_2[n]$ by sampling $x_2(t)$ with $f_{s_2}$ and measure $\phi_3[n]$ by sampling  $f_{s1}$ with $f_{s2}$ and from those measurements resolve $\theta_2-\theta_1$.

Similarly if your application is for a phase offset that would not be changing, then you can measure $f_{s2}$ error using the slope of the result without having to sample $f_{s1}$ with $f_{s2}$.

The true results will depend on the actual clock accuracy of $f_{s1}$ but we have completely removed $f_{s2}$ from the equation. If you can consider $f_{s1}$ to be your true reference of time, meaning it is accurate enough for the precision and accuracy of your measurement, then the result will be the phase difference of your two waveforms being sampled. This means that ultimately you need something to be your common reference of time.


Unlikely.  The top 6 FFT magnitudes would more probably contain several overtones or harmonics that are stronger than the 6 note pitches in a chord.  Furthermore, some of the spectral frequency peaks can be between FFT bins, depending on the length of the FFT.  So the frequency you've calculated by just multiplying by a bin index could be different from any pitch frequency in the audio.  Depending on any window used, you might also find windowing artifacts in the FFT magnitude spectrum.


why do we need the particle filter?

When your model is non-linear there is no analytical solution to obtain the PDF of the state. The Kalman filter/smoother can't be applied in the sense that it is not optimal. The particle filter is one method to obtain an estimate of state PDF in a sequential setup.

  can't we just do ML estimation

Dependent on the model you may find that even an ML estimator can't be found without going through some approximations. But the particle filter is also used in a Bayesian contexts where other estimators are relevant.


I think I know the answer now. If I made an error, please correct me.

Assume there are $m$ sensors and $d$ sources. As the eigenvectors in $S$ and $G$ of the Hermitian matrix $R$ is selected to be orthogonal and normalized,
$$UU^*=\begin{bmatrix}S&G\end{bmatrix}\begin{bmatrix}S^*\\G^*\end{bmatrix}=SS^*+GG^*=UU^{-1}=I_m;$$
$$S^*S=I_d;$$
$$\mathrm{span}\{G\}=\mathrm{span}\{S\}^{\bot}.$$

We know the noise subspace $\mathrm{span}\{G\}$ and $\mathrm{span}\{A\}$ are orthogonal as in the MUSIC algorithm:
$$\mathrm{span}\{G\}=\mathrm{span}\{A\}^{\bot}.$$
So $\mathrm{span}\{A\}=\mathrm{span}\{S\}$, which implies their projection matrices into this subspace are equal:
$$A(A^*A)^{-1}A^*=S(S^*S)^{-1}S^*.$$
Using $S^*S=I_d$, we get
$$A(A^*A)^{-1}A^*=SS^*;$$
with $SS^*+GG^*=I_m$, we get
$$GG^*=I_m-SS^*=I_m-A(A^*A)^{-1}A^*.$$


Pick up an apple and raise it one meter. You have just spent (approximately) one joule of energy. You had to spend energy because the action (raising the apple) goes against a force (Earth's gravity).

The same is true for electrons. An electrical signal is nothing but electrons moving in an electric field. Causing these electrons to move takes energy, exactly like moving a train or an apple.

In the formula you posted, $x(t)$ is a voltage causing electrons to move in a resistance of $1\,\Omega$, and $E_s$ is the energy required to make those electrons move around in the resistor, causing a current to develop.


What is the right way to interpolate a 2D grid?

What a correct solution constitutes depends on the application. What you have proposed sounds more like an inference or estimation problem (that needs to be fleshed out). It would make sense if you thought your values were noisy.

If that were the case, try looking at Kriging.

https://en.m.wikipedia.org/wiki/Kriging

If you just want to interpolate MATLAB has interp2

Which (not limited to the  2 types mentioned above) is correct depends on what your ultimate purpose is.


Let me answer your two questions in turn:

For your first question, generally, yes that is correct; if you filter a 10KHz sinewave that has a DC offset with a filter that has a cutoff frequency below the frequency of the sinewave, then the sinewave would be rejected. The amount of rejection specifically depends on the performance of the filter, but given you said you have a 1KHz cutoff frequency, the sinewave is significantly higher and therefore sufficiently rejected.

I see in your plot that the order of the filter is 5, which for a Butterworth filter as also shown would have a rejection of 20dB/decade *5 (where 5 is the order of your filter), or 100 dB per decade. Depending on other factors such as your digital dynamic range, this suggests that you would be able to filter your 10KHz sine wave up to 100 dB (10KHz is a decade above the cutoff frequency).  The DC signal, which is below the cutoff frequency would pass through to the output, unless something in your system blocked DC or introduced other DC -offsets (which is possible).  Also the filter itself can have gain or loss, so the actual DC output level if it did pass through can be modified by this gain or loss accordingly.

For your second question, sampling frequency is the sampling rate for the signals passing through this digital filter implementation.  From the figure, you are using a sampling rate of 200KHz, and yes this would be the sampling rate of the sinewave that is created.

I hope this helped to clear up some of your questions.


Calculate the time delay introduced by group delay for IIR-Filters

If you want to compute the delay experienced by a sine wave at a given frequency $\omega_0$, then it's not the group delay but the phase delay that you need to consider. In the case of linear phase FIR filters the phase delay and the group delay are identical, but in the general case they aren't.

The phase delay is defined by

$$\tau_p(\omega)=-\frac{\phi(\omega)}{\omega}\tag{1}$$

where $\phi(\omega)$ is the frequency dependent (unwrapped) phase of the system's frequency response. Clearly, if $\phi(\omega)=c\cdot \omega$ with some constant $c$, then the phase delay $(1)$ equals the negative derivative of the phase, i.e., group delay and phase delay are equal.

Of course, the (de-)normalization by the sampling frequency is always done in the same way. If $\omega$ is normalized, then you get the phase delay in samples, otherwise in seconds.

EDIT:

The response of an LTI system with frequency response $H(\omega)=M(\omega)e^{j\phi(\omega)}$ to a sinusoidal input signal $x(t)=A\sin(\omega_0t+\theta)$ is given by

$$\begin{align}y(t)&=AM(\omega_0)\sin(\omega_0t+\theta+\phi(\omega_0))\\&=AM(\omega_0)\sin\left[\omega_0\left(t+\frac{\phi(\omega_0)}{\omega_0}\right)+\theta\right]\end{align}\tag{2}$$

So the delay experienced by $x(t)$ is given by $\tau_p(\omega_0)$ as defined in $(1)$.


Is the filter described by specific difference equation time independent

From the given difference equation, nothing can be said about linearity or time-invariance. Note that the difference equation (DE) doesn't even uniquely specify the output for a given input $x[n]$. To see this, assume that $y_1[n]$ is a solution to the DE. Let $y_h[n]$ be a solution to the corresponding homogeneous equation

$$y[n]+2y[n-1]+y[n-2]=0\tag{1}$$

Then any sequence $y_2[n]=y_1[n]+y_h[n]$ must also be a solution to the original DE. Note that if $y_h[n]\neq 0$, the corresponding system is neither linear nor time-invariant because that part of the output signal does not depend on the input $x[n]$. So no matter if you shift the input signal, or if you multiply it with a constant, that term remains unchanged.

In order to uniquely specify $y[n]$ for a given input $x[n]$, we need auxiliary conditions. There are several options, but the auxiliary conditions that make sure that the system described by the DE is linear and time-invariant (and causal) are initial rest conditions. If we assume that $x[n]=0$ for $n<0$, then the initial rest conditions are $y[-1]=y[-2]=\ldots=0$. This makes sure that $y_h[n]=0$ for $n\ge 0$, and consequently, the output $y[n]$ contains no term that is independent of the input $x[n]$.

A system described by a linear DE with constant coefficients (such as the one in your question) describes a linear and time-invariant (and causal) system only in combination with the aforementioned initial rest conditions.


If you're averaging multiple (two, in your case) spectral magnitude sequences your result will be valid, correct, as long as the sequences are equal in length.  Now whether the "averaged spectral magnitude" result gives you useful information, that you will have to decide on your own.


Don't understand what is meant by signal dimension

In digital communication systems, the dimension of a modulation scheme refers to the number of basis function, i.e. independent/orthogonal signals (in this case $\sin$ and $\cos$ function) used to represent the symbols. In these systems $k=\log_2 M$ binary digits are mapped into analog waveforms below
$$\left\{s_m(t), m = 1, 2, \ldots, M \ \big\vert \ M = 2^k\right\}$$
Take a PAM signal for instance, this has the form
\begin{align}
s_m(t)&=\Re\left[A_mg(t)\exp\left(j2\pi f_c t\right)\right]\\
&= A_mg(t)\cos\left(2\pi f_c t\right)
\end{align}
One basis function (one axis), the in-phase component (here $x$-axis) is used to represent the signal. The signal is one dimensional. Now consider PSK modulation, the signal waveforms are represented as
\begin{align}
s_m(t)&=\Re\bigg[g(t)\exp\left(j2\pi\frac{m - 1}{M}\right)\exp\left(j2\pi f_c t\right)\bigg]\\
&= g(t)\cos\left(2\pi f_c t + 2\pi\frac{m - 1}{M}\right)
\end{align}
Which can be decomposed down into sine and cosine components. This is a two-dimensional signal. The same goes for QAM signals.

A common misconception is taking the dimension for the number of bits per symbol. In brief, the dimension is the number of orthogonal basis functions used in the linear combination representing the waveforms $s_m (t)$. Note that despite the digital part of it, what is not always aparent is that the modulator takes in charge the mapping of the binary info into the symbols in the analog domain ready for transmission.

Having said that, there is also a notion of multi-dimensional signals where the time and frequency intervals are subdivided into short intervals for transmission. More on this and my explanation above see $[1]$.

$[1]$, J. G. Proakis, Digital Communications, 4th edition, McGraw Hill, chapter 4.