Subsampling property of wavelet transform

I am not completely confident with the full theory of wavelet frames, however I will share some pointers and my beliefs, do not take any of  the following for a truth. The property you show is (partly) a characterization of an orthogonal system. In Hernandez & Weiss, A first course on wavelets, 1996, you find Theorem 2.12:

If $\psi\in L^2(\mathbb{R})$ is a band-limited function such that $\hat{\psi}$ is zero in a neighborhood of the origin and $\{ \psi_{j,k}: j,k\in \mathbb{Z}\}$ is an orthonormal system, then this system is coomplete if and only if
$$\sum_{j\in \mathbb{Z}} |\hat{\psi}(2^j \xi)|^2=1 \quad \text{a.e. on }\mathbb{R}-\{0\}$$
and
$$\sum_{j=0}^{\infty} |\hat{\psi}(2^j \xi)\hat{\psi}(2^j (\xi+2k\pi))| =0\quad \text{a.e. on }\mathbb{R},\quad k\in 2\mathbb{Z}+1$$
with the system defined as:
$$ \psi_{j,k}(x)= 2^{\frac{j}{2}}\psi(2^j x-k)\,.$$
So the conditions are a little bit stricter for an orthonormal system.
So it seems rather to me a characterization in 1D (f band-limited wavelets), than can help define radial 2D wavelets, than a characterization of radial wavelet per se. So if I understand your first question well, this imply the possibility of discrete wavelets and downsampling.

Since it is a characterization, I (very unsure) would say that you can downsample, but you have no guarantee about the completeness. I do not understand your condition  $$h(\omega)=1, \forall \omega>\pi/2, h(\omega)=1, \forall\omega>\pi/2,$$
but I can add that in Chapter 3.4 of the above book, you have characterization of wavelets with support in $[-\frac{8}{3}\pi,-\frac{2}{3}\pi] \cup [\frac{2}{3}\pi,\frac{8}{3}\pi] $, so I would say yes, as $\frac{8}{3}> \frac{1}{2}$. I would say that the limit on $\pi/2$ is not an issue because of the dilations of the mother wavelet.

For the last question, without further information on the $h_i$, it is unlikely that you generate a wavelet system (possibly because of refinement equations), but as long as you cover the whole spectrum, who really cares? You get a decomposition that is energy preserving, but I am not sure yet it suffices to define a tight frames.

Additional lectures:

Affine Density, Frame Bounds, and the Admissibility Condition for Wavelet Frames, 2007, Gitta Kutyniok, with results in frames with generalized samplings  $ \psi_{j,k}(x)= a^{\frac{j}{2}}\psi(a^j x-b)$
A first course on wavelets, 1996, Eugenio Hernandez and Guido L. Weiss
Ten Lectures on Wavelets, 1992, I. Daubechies, esp. Chapter 3
Tight Wavelet Frame Construction and its application for Image Processing, 2005, Kyunglim Nam
Frames and bases. An introductory course, 2008, Christensen, O.


Is a good practice to use one PID in rise time and another PID for steady time in order to eliminate oscillations during steady time?

Switching between 2 PID controller is called gain-scheduling. There are various ways to implement gain scheduling. But basically, the idea is that the gain should change smoothly. You should not change suddenly. You could have one set of gains if T < T1, another set if T > T2 and use interpolation between the 2 sets if T1 < T  < T2.

You can also use 2 PID controller with the unused PID set in tracking mode while the other controls the actuator also. However, it will be harder to have a smooth transition between the 2 controllers.

Secondly, before implementing gain scheduling, make sure your issue is really due to non-optimal gains. Typically, there are 2 problems with control loop for temperatures

1 - Actuator saturation which causes integrator windup. This happens often in temperature control when the difference between the setpoint and the measured temperature is significant. This can be fixed with an appropriate anti-windup scheme.

2 - Asymetry between heating and cooling. Typically, if you have a heater, you can only "heat" you cannot cool your system. This makes the choice of gain a bit tricky as the cooldown rate is not controllable by your controller. If you have some kind of Peltier cooler, you can both cool down and heat up, but the cooling rate is usually slower.


how were DTMF frequencies determined exactly and how can I extend them?

From page 293 of Engineering and Operations in the Bell System (1982-1983) by AT&T Bell Labs:

  The tones have been carefully selected to minimize harmonic interference and the probability that a pair of high and low tones will be simulated by the human voice, thus protecting network control signaling.

P.S. Yes, I know this answer is 6 years late, I just happened to be researching for the answer and I thought this could help.


The 2D version is just the same as the 1D version:
$$
x(n,m) \leftrightarrow X(k,l) = \sum_{n=0}^{N-1} \sum_{m=0}^{M-1} x(n,m) e^{-j2\pi nk/N - j2\pi ml/M}\\
x(n-a,m-b) \leftrightarrow
\begin{array}{ccl}
&&\sum_{n=0}^{N-1} \sum_{m=0}^{M-1} x(n-a,m-b) e^{-j2\pi nk/N - j2\pi ml/M} \\
&=& \sum_{n'=0}^{N-1} \sum_{m'=0}^{M-1} x(n',m') e^{-j2\pi (n'+a)k/N - j2\pi (m'+b)l/M}\\
 &=& X(k,l) e^{-j (2\pi ka/N + 2\pi lb/M)}
\end{array}
$$


a-weighting signal without digital filter

The FFT has a processing gain of 10Log(N) that would apply to your measurements of the standard deviation, where N is the length of the signal, and the Hamming Window has a coherent gain of 0.54 (which is a loss). Refer to fred harris' famous paper On the Use of Windowing that details these gains and losses.

Let's see if this works out. You have 44100 samples, so the net processing gain is $10Log10(44100\times.54) = 43.8 \text{ dB}$

Your peak sinewave is 32768 so the rms should be $20log_{10}(32768/\sqrt{2})=87.3 \text { dB}$ which is consistent with your measurement.

After taking the window and the FFT, you saw the RMS signal increase to 130 dB, while you measured 87 dB before the FFT. 130-87 = 43 dB.

The difference of 43 dB is consistent with that predicted by the formula for the net processing gain given above.

The following is more details that may further help to understand how this occurs:

The coherent gain of your input signal is N. Your input signal in time is

$$A sin(\omega t) = \frac{A}{2}\frac{e^{j\omega t}}{j}-\frac{A}{2}\frac{e^{j\omega t}}{j}$$

The important item to see from Euler's expansion above is the sine is composed of two exponentials each with an amplitude of half of the sine wave itself. After the FFT, each of these is increased by a factor of N due to the coherent processing gain in the FFT.

Thus each of the two tones in the FFT would have a magnitude of $\frac{AN}{2}$ (if you sampled coherently at an exact integer number of samples this would be exact, otherwise we will see it slightly reduced from this due to "scalloping loss").

If sampled coherently, every other bin besides the two tones described will be zero (we could further show through Parseval's theorem the general case where all the energy in the signal is conserved, but to keep my explanation simple I will describe just the case of coherent sampling, meaning the sampling clock is an integer multiple of the sine wave).

Therefore when computing the rms, since there are only 2 non-zero values but we are using N samples, the rms calculation of the FFT would reduce to:

$$\sigma = \sqrt{\frac{\left(\frac{AN}{2}\right)^2+\left(\frac{AN}{2}\right)^2}{N}}$$

Which reduces to

$$\sigma = \sqrt{N}\frac{A}{\sqrt{2}}$$

Which we recognize as the $\sqrt{N}$ times the rms of our input sinusoid.  When we window the signal before taking the FFT, the coherent amplitude (that previously increased by N) would now be reduced by the coherent gain of the window, which for the case of the Hamming window was 0.54.

Thus the standard deviation of the FFT is expected to be $\sqrt{(N\times0.54)}$ times larger than the standard deviation of the input sinusoid for the case of a Hamming Window being used. So in dB:

$$20log_{10}(\sigma) = 10log_{10}(N\times0.54) $$


here is code for a "rolling median" copied from the reference of the reference i made above.  i believe the computational cost of the alg (per sample) is $O(\log_2(N))$.

//Copyright (c) 2011 ashelly.myopenid.com under <http://w...content-available-to-author-only...e.org/licenses/mit-license>

#include <stdlib.h>

//Customize for your data Item type
typedef int Item;
#define ItemLess(a,b)  ((a)<(b))
#define ItemMean(a,b)  (((a)+(b))/2)

typedef struct Mediator_t
{
   Item* data;  //circular queue of values
   int*  pos;   //index into `heap` for each value
   int*  heap;  //max/median/min heap holding indexes into `data`.
   int   N;     //allocated size.
   int   idx;   //position in circular queue
   int   ct;    //count of items in queue
} Mediator;

/*--- Helper Functions ---*/

#define minCt(m) (((m)->ct-1)/2) //count of items in minheap
#define maxCt(m) (((m)->ct)/2)   //count of items in maxheap

//returns 1 if heap[i] < heap[j]
int mmless(Mediator* m, int i, int j)
{
   return ItemLess(m->data[m->heap[i]],m->data[m->heap[j]]);
}

//swaps items i&j in heap, maintains indexes
int mmexchange(Mediator* m, int i, int j)
{
   int t = m->heap[i];
   m->heap[i]=m->heap[j];
   m->heap[j]=t;
   m->pos[m->heap[i]]=i;
   m->pos[m->heap[j]]=j;
   return 1;
}

//swaps items i&j if i<j;  returns true if swapped
int mmCmpExch(Mediator* m, int i, int j)
{
   return (mmless(m,i,j) && mmexchange(m,i,j));
}

//maintains minheap property for all items below i/2.
void minSortDown(Mediator* m, int i)
{
   for (; i <= minCt(m); i*=2)
   {  if (i>1 && i < minCt(m) && mmless(m, i+1, i)) { ++i; }
      if (!mmCmpExch(m,i,i/2)) { break; }
   }
}

//maintains maxheap property for all items below i/2. (negative indexes)
void maxSortDown(Mediator* m, int i)
{
   for (; i >= -maxCt(m); i*=2)
   {  if (i<-1 && i > -maxCt(m) && mmless(m, i, i-1)) { --i; }
      if (!mmCmpExch(m,i/2,i)) { break; }
   }
}

//maintains minheap property for all items above i, including median
//returns true if median changed
int minSortUp(Mediator* m, int i)
{
   while (i>0 && mmCmpExch(m,i,i/2)) i/=2;
   return (i==0);
}

//maintains maxheap property for all items above i, including median
//returns true if median changed
int maxSortUp(Mediator* m, int i)
{
   while (i<0 && mmCmpExch(m,i/2,i))  i/=2;
   return (i==0);
}

/*--- Public Interface ---*/

//creates new Mediator: to calculate `nItems` running median.
//mallocs single block of memory, caller must free.
Mediator* MediatorNew(int nItems)
{
   int size = sizeof(Mediator)+nItems*(sizeof(Item)+sizeof(int)*2);
   Mediator* m=  malloc(size);
   m->data= (Item*)(m+1);
   m->pos = (int*) (m->data+nItems);
   m->heap = m->pos+nItems + (nItems/2); //points to middle of storage.
   m->N=nItems;
   m->ct = m->idx = 0;
   while (nItems--)  //set up initial heap fill pattern: median,max,min,max,...
   {  m->pos[nItems]= ((nItems+1)/2) * ((nItems&1)?-1:1);
      m->heap[m->pos[nItems]]=nItems;
   }
   return m;
}

//Inserts item, maintains median in O(lg nItems)
void MediatorInsert(Mediator* m, Item v)
{
   int isNew=(m->ct<m->N);
   int p = m->pos[m->idx];
   Item old = m->data[m->idx];
   m->data[m->idx]=v;
   m->idx = (m->idx+1) % m->N;
   m->ct+=isNew;
   if (p>0)         //new item is in minHeap
   {  if (!isNew && ItemLess(old,v)) { minSortDown(m,p*2);  }
      else if (minSortUp(m,p)) { maxSortDown(m,-1); }
   }
   else if (p<0)   //new item is in maxheap
   {  if (!isNew && ItemLess(v,old)) { maxSortDown(m,p*2); }
      else if (maxSortUp(m,p)) { minSortDown(m, 1); }
   }
   else            //new item is at median
   {  if (maxCt(m)) { maxSortDown(m,-1); }
      if (minCt(m)) { minSortDown(m, 1); }
   }
}

//returns median item (or average of 2 when item count is even)
Item MediatorMedian(Mediator* m)
{
   Item v= m->data[m->heap[0]];
   if ((m->ct&1)==0) { v= ItemMean(v,m->data[m->heap[-1]]); }
   return v;
}

/*--- Test Code ---*/
#include <stdio.h>
void PrintMaxHeap(Mediator* m)
{
   int i;
   if(maxCt(m))
      printf("Max: %3d",m->data[m->heap[-1]]);
   for (i=2;i<=maxCt(m);++i)
   {
      printf("|%3d ",m->data[m->heap[-i]]);
      if(++i<=maxCt(m)) printf("%3d",m->data[m->heap[-i]]);
   }
   printf("\n");
}
void PrintMinHeap(Mediator* m)
{
   int i;
   if(minCt(m))
      printf("Min: %3d",m->data[m->heap[1]]);
   for (i=2;i<=minCt(m);++i)
   {
      printf("|%3d ",m->data[m->heap[i]]);
      if(++i<=minCt(m)) printf("%3d",m->data[m->heap[i]]);
   }
   printf("\n");
}

void ShowTree(Mediator* m)
{
   PrintMaxHeap(m);
   printf("Mid: %3d\n",m->data[m->heap[0]]);
   PrintMinHeap(m);
   printf("\n");
}

int main(int argc, char* argv[])
{
   int i,v;
   Mediator* m = MediatorNew(15);

   for (i=0;i<30;i++)
   {
      v = rand()&127;
      printf("Inserting %3d \n",v);
      MediatorInsert(m,v);
      v=MediatorMedian(m);
      printf("Median = %3d.\n\n",v);
      ShowTree(m);
   }
}


How do I interpret this block diagram correctly?

It is better to "parse" these networks from the output back towards the input, calling their input some general $x$ and performing substitutions and/or compositions.

So, let's call these networks $U$pper and $L$ower.

From the upper diagram:

$$UH_2[n] = x[n] + x[n-1] \cdot -a_1$$

and

$$UH_1[n] = x[n] \cdot b_0+x[n-1] \cdot b_1$$

Now, the output of $U$pper is given by the composition of the two:

$$UH_y[n] = UH_2[UH_1[n]]$$

This is because, the output of one, becomes the input to the other.

Or...

$$UH_y[n] = UH_1[n] + UH_1[n-1] \cdot -a_1$$

...and if you substitute...

$$UH_y[n] = x[n] \cdot b_0 + x[n-1] \cdot b_1 + (x[n-1] \cdot b_0 + x[n-2] \cdot b_1) \cdot -a_1 \Rightarrow \\ x[n] \cdot b_0 + x[n-1] \cdot b_1 + x[n-1] \cdot b_0 \cdot -a_1 + x[n-2] \cdot b_1 \cdot -a_1$$

Notice here that if you are trying to get to the $-1$ of the $x[n-1]$, then that would be the $x[n-2]$.

And this concludes the $U$pper part.

For the $L$ower part, we are not going to go through the whole thing, because of two reasons:

If you notice, the $H2, H1$ are reversed (The $L$ower network calls $H2$ what the $U$pper network calls $H1$).
We have already mapped $H2, H1$.

So, the lower network's response is:

$$LH_y[n] = LH_1[LH_2[n]]$$

Or...

$$LH_y[n] = LH_2[n] \cdot b_0+LH_2[n-1] \cdot b_1$$

And if you substitute:

$$LH_y[n] = (x[n]+ -a_1 \cdot x[n-1]) \cdot b_0 + (x[n-1]+ -a_1 \cdot x[n-2]) \cdot b_1 \Rightarrow \\
   x[n] \cdot b_0 + -a_1 \cdot b_0 \cdot x[n-1] + x[n-1] \cdot b_1 + b_1 \cdot -a_1 \cdot x[n-2]$$

Which, after you re-arrange, looks exactly the same.

Hope this helps.


Help calculating / understanding the MFCCs: Mel-Frequency Cepstrum Coefficients

Step by step...

1. & 2. This is correct. Note that the frames are usually overlapping, for example, frame 0 is samples 0 to 440; frame 1 is samples 220 to 660; frame 2 is samples 440 to 880 and so on... Note also that a window function is applied to the samples in the frame.

3. The Fourier transform is done for each frame. The motivation behind this is simple: a speech signal varies over time, but is stationary over short segments. You want to analyze each short segment individually - because on this segments the signal is simple enough to be described efficiently by few coefficients. Think of someone saying "hello". You don't want to see all phonemes collapsed into one single spectrum (FFT collapses temporal information) by analyzing all the sound at once. You want to see "hhhhheeeeeeeeeeelloooooooooo" to recognize the word stage by stage, so it has to be broken down into short segments.

4. "Mapping to the Mel scale" is misleading and that's probably why you're getting confused. A better description for this step would be: "Compute the signal energy through a bank of filters tuned to mel-scaled frequencies". Here is how this is done. We consider $N$ frequencies (a commonly used value is $N = 40$) equally spaced according to the mel scale, between 20 Hz (bottom of the hearing range) and the Nyquist frequency. Practical example: The signal is sampled at 8kHz and we want 40 bins. Since 4kHz (Nyquist) is 2250 mel, the filterbank center frequencies will be:
0 mel, 2250 / 39 mel, 2 x 2250 / 39 mel .. 2250 mel.

Once these frequencies have been defined, we compute a weighted sum of the FFT magnitudes (or energies) around each of these frequencies.

Look at the following picture, representing a filter bank with 12 bins:

The 8th bin has a center frequency of around 2kHz. The energy in the 8th bin is obtained by summing weighted FFT energies in the 1600 to 2800 Hz approximately - with the weight peaking at around 2kHz.

Implementation note: This bunch of weighted sums can be done in a single operation - a matrix multiplication of a "filterbank matrix" by the FFT energies vector.

So at this stage we have "summarized" the FFT spectrum into a set of 40 (12 in the illustration) energy values, each corresponding to a different range of frequencies. We take the log of these values.

5. The next step consists in taking the DCT of this sequence of 40 log-energies. This will yield 40 values. The first $K$ coefficients are the MFCC (Usually, $K = 13$). Actually, the very first DCT coefficient is the sum of all the log-energies computed at the previous step - so it is an overall measure of signal loudness and is not very informative as to the actual spectral content of the signal - it is often discarded for speech recognition or speaker id applications where the system has to be robust to loudness variations.


When you quantize a signal, you introduce and error which can be defined as
$$q[n] = x_q[n]-x[n]$$
where $q[n]$ is the quantization error, $x[n]$ the original signal, and $x_q[n]$ of the quantized signal. The maximum quantization error is simply $max(\left | q \right |)$, the absolute maximum of this error function. Dx in this definition seems to be the range of the input signal so we could rewrite this as
$$Q = \frac{max(x)-min(x)}{2^{N+1}}$$

Let's look at a quick example. Let's assume you have a signal that's uniformly distributed between -1 and +1 and you want to quantize this with 3 bits. You have a total 8 of quantizaton steps which would map to [-1 -.75 -.5 -25 0 .25 .5 .75]. The difference between steps is 0.25. If you round during quantization the maximum error will be half of that (i.e. 0.125).

Now let's try the formula:
$$Q = \frac{max(x)-min(x)}{2^{N+1}} = \frac{1-(-1)}{2^{3+1}}= \frac{2}{16} = 0.125$$


Interpolation of a gray scale image in matlab

Try the imfill function, see http://www.mathworks.com/help/toolbox/images/ref/imfill.html


You could normalize the average amplitude, i.e.

YData = YData / mean(abs(YData));

Or you could normalize the signal power to one, i.e.

YData = YData / sqrt(mean(abs(YData).^2));

If just the peaks are bothering you, you could use dynamic range compression, but that would introduce nonlinear distortions. As Phonon hinted, please tell us why you are not happy with the results and why you want to normalize your data in the first place. That would allow for a more precise answer.


Does zero-forcing equalizer need known channel impulse response?

If the received signal can be written as

$$\mathbf{y} = \mathbf{H}\,\mathbf{x} + \mathbf{n}$$

where $\mathbf{H}$ is the channel matrix, $\mathbf{x}$ is the transmitted vector, and $\mathbf{n}$ is the AWGN of the channel, then a zero forcing equalizer is simply (assuming that the channel matrix is square, and it's estimated perfectly at the receiver)

$$\mathbf{H}^{-1}\mathbf{y} = \mathbf{x}+\mathbf{H}^{-1}\mathbf{n}$$

Obviously, you need the channel impulse response, which is captured in the channel matrix. This channel matrix is estimated in practice using any channel estimation technique, but the estimation is usually not perfect, and thus the aforementioned ZF equalizer serves as a theoretical limit.


White noise generation for different bandwidth signals in MATLAB

One technicality to be aware of is that white Gaussian noise has infinite bandwidth. What you are asking for is filtered white noise, which has a flat PSD in a given bandwidth, and zero power outside of it. It is common to call this filtered noise "white noise", but you should be aware of the difference.
With that out of the way, consider the channel between the ground station and the satellite. The noise in this channel is modeled as white Gaussian, with $\text{PSD} = N/2$. Assuming a filter with bandwidth $B$, the power of the filtered noise is, as you say, $NB$ (because the power is the integral of the PSD from minus infinity to infinity, which includes negative frequencies).
Now you need to define the receiver filter. It is very unusual, in my experience, to transmit a signal with two different bandwidths. You have two choices for the receiver design: a single filter with bandwidth 180 kHz, or two filters, one of 180 kHz and another of 15 kHz.
The single-filter design is obviously cheaper and much simpler to implement. The drawback is that it is quite suboptimal, as the data signal will be subject to much more noise than it would with a 15 kHz filter.
The two-filter design will be much more complex (I have never seen a real system implemented like this). The benefit is that it will be optimal.
In Matlab, you can generate filtered white noise of variance 1 with randn. After that, it's just a matter of scaling the noise appropriately to get the variance you want.


Any (melody) musical instrument note has a fundamental frequency (think "a sine wave") which our ears tend to interpret as a "pitch", and a number of harmonics (extra sine waves) which are integer multiples of that fundamental frequency. These extra sine waves with double, triple, ... times the fundamental frequency add up to something that our ears tend to interpret as a timbre. The reason for why you find this special relation between fundamental frequency and harmonics is a consequence of the physics of how real instruments make sound.

The difference between an "A4" note on a piano and an "A4" note on a violin is in which harmonics are present in the sound, and how strongly they are present. One way to examine this is by comparing their respective spectrograms. Several audio tools allow you to draw a spectrogram for a given sound (e.g. sonic visualizer or audacity).

For percussive instruments there's no easy relation between fundamental frequency and harmonics and the spectrum typically looks much noisier (many more high frequency components). Also the spectrum at the start of a melody instrument note (the so-called "attack") will typically look noisier.


Unexpected frequency components found in Periodogram

First of all, if you use a logarithmic graph, it is customary to use dB. In order to convert a FFT in decibels, we use 20*log10, not log10.

Second of all, the "unwanted" frequencies are 10^-30 times smaller than the wanted frequencies, i.e. they are insignificant. . The unwanted frequencies are most likely due to the limited precision of floating-point arithmetic and not due to the sampling process.

https://en.wikipedia.org/wiki/Floating-point_arithmetic

That should answer your second point too...


A non-linear phase response is equivalent to group delay variation, and similar to multi-path fading can cause inter-symbol interference if it is severe enough, but either can be corrected with channel equalization. In this case it would be a static correction so would not be overly complicated to implement in the transmitter, receiver or both. Factoring the equalization into the cascade of two equal filters and then implementing it in both the transmitter and receiver would also provide the benefit of SNR optimization (against white noise) as a matched filter (but the actual performance benefit of doing that may be negligible compared to the added complication, so just throwing out the thought).

The metric to analyze would be group delay variation since that would show it directly. A linear phase component is only a delay and does not provide any distortion, so the group delay which is the derivative of the phase with respect to frequency would convert any linear phase to a constant value and highlight the variation which is the destructive component. One approach is to simulate the filter response and analyze the eye diagram to assess the signal quality and if further equalization and or changes to the analog filter design would be needed.

The graphic below depicts this for the case of multipath, but multipath is simply a filter with phase (and amplitude) variation across the signal bandwidth.


Denoising short, non-uniformly spaced, bandlmited sequences

Not a lot of references that include nonuniform, band limited, and filtering

D. Bonacci and B. Lacaze, "Lowpass/bandpass signal reconstruction and digital filtering from nonuniform samples," 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), South Brisbane, QLD, 2015, pp. 3626-3630.
doi: 10.1109/ICASSP.2015.7178647
Abstract: This paper considers the problem of non uniform sampling in the case of finite energy functions and random processes, not necessarily approaching to zero as time goes to infinity. The proposed method allows to perform exact signal reconstruction, spectral estimation or linear filtering directly from the non-uniform samples. The method can be applied to either lowpass, or bandpass signals.
keywords: {band-pass filters;digital filters;linear phase filters;low-pass filters;random processes;signal reconstruction;signal sampling;band-pass signal reconstruction;digital filter;finite energy function;linear filter;low-pass signal reconstruction;nonuniform sampling problem;random process;spectral estimation;Baseband;Estimation;Fourier series;Interpolation;Nonuniform sampling;Random processes;Time-frequency analysis;Nonuniform filtering;Periodic nonuniform sampling;Sampling theory;Signal reconstruction;Spectral estimation},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=7178647&isnumber=7177909

H. Johansson and P. Lowenborg, "Reconstruction of nonuniform sampled bandlimited signals using digital fractional filters," ISCAS 2001. The 2001 IEEE International Symposium on Circuits and Systems (Cat. No.01CH37196), Sydney, NSW, 2001, pp. 593-596 vol. 2.
doi: 10.1109/ISCAS.2001.921140
Abstract: This paper considers the problem of reconstructing nonuniformly sampled bandlimited signals using a synthesis system composed of digital fractional delay filters. The overall system can be viewed as a generalization of time-interleaved ADC systems. By generalizing these systems, it is possible to eliminate the errors that are introduced in practice due to time-skew errors
keywords: {analogue-digital conversion;bandlimited signals;digital filters;filtering theory;signal reconstruction;digital fractional filters;nonuniform sampled bandlimited signals;signal reconstruction;synthesis system;time-interleaved ADC systems;time-skew errors;Analog-digital conversion;Delay;Digital filters;Electronic mail;Filter bank;Frequency;Interleaved codes;Nonuniform sampling;Sampling methods;Signal synthesis},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=921140&isnumber=19923

Y. C. Eldar and A. V. Oppenheim, "Filter bank interpolation and reconstruction from generalized and recurrent nonuniform samples," 2000 IEEE International Conference on Acoustics, Speech, and Signal Processing. Proceedings (Cat. No.00CH37100), Istanbul, 2000, pp. 324-327 vol.1.
doi: 10.1109/ICASSP.2000.861961
Abstract: This paper introduces a filter bank interpretation of various sampling strategies, which leads to efficient interpolation and reconstruction methods. An identity, referred to as the interpolation identity, is used to obtain particularly efficient discrete-time systems for interpolation to uniform Nyquist samples, either for further processing in that form or for conversion to continuous-time. The interpolation identity also leads to a new class of sampling theorems including an extension of Papoulis' (1997) generalized sampling expansion
keywords: {channel bank filters;discrete time filters;filtering theory;interpolation;signal reconstruction;signal sampling;discrete-time systems;filter bank interpolation;generalized sampling expansion;interpolation identity;recurrent nonuniform samples;sampling theorems;signal reconstruction;uniform Nyquist samples;Data acquisition;Digital signal processing;Filter bank;Fourier transforms;Interpolation;Iterative algorithms;Laboratories;Sampling methods;Signal processing;Signal sampling},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=861961&isnumber=18684


Why doesn't this complex multiplication in the frequency domain produce my expected phase shift?

If you want a phase-shifted strictly-real result, then you have to make sure the data you feed to the IFFT is conjugate symmetric.  So make sure to reflect the complex conjugate of any phase changes you make.  e.g. Make the upper or negative half of your IFFT input vector mirror the complex conjugate of the bottom or positive half of the input vector.

Added: (Also, "mirror" usually means the array index goes in the opposite direction.)


How to Calculate the number of DFTs taken in a STFT

I will first produce a lower number based on the covering of all the non-zero sample set.

With a window of size $W$, counting samples from $0$ to $S-1$, the first window containing at least  one sample covers $[-W+1,0]$. If the discrete overlap range is a number of samples $O<W$, the next window will be $$[1-O,W-O]=[-W+1,0]+(W-O)$$ and the others $k$th hops are $$[-W+1,0]+k(W-O)$$
Those intervals exceed the last sample $S-1$ when $0+k(W-O)\ge S-1$ or

$$k\ge \left\lceil\frac{S-1}{W-O}\right\rceil$$

so the covering takes $ \left\lceil\frac{S-1}{W-O}\right\rceil+1$ windows, hence DFTs.

Does this suffice in general, as we did not use any window property? For instance, if the window is non-null at one point only (or vanishes), the $O$ cannot be chosen arbitrarily.

The question in general seems open to me. For randomly chosen complex filterbanks, the smallest oversampling seems to suffice for recovery. For other tasks, who knows?


Most common modern day Image Segmentation techniques

I do not know about many segmentation techniques, but I have been dealing with structures that offer a "choice" of segmentation pieces that can be further examined to produce a satisfying segmentation.

Hopefully somebody else can write about some different state-of-the-art segmentation method that I don't know much about.

A small introduction as to why it is nice to have choices, or levels, of segmentation for a different image: segmentation is an ill defined problem. The ground truth does not exist: the desired results always depend on the user requirements and specifications. A nice quote:

  Even for a fixed image, there may be more than one "best" segmentation because the criteria defining the quality of a segmentation are application dependent. This motivated us to focus our research on image-partitioning techniques delivering "puzzle pieces" that can be used (...) to generate a segmentation satisfying specific user requirements

(P. Soille: Constrained Connectivity for Hierarchical Image Partitioning and Simplification (2008))

There are hierarchical structures, hierarchical image decompositions that propose image partitions with varying complexities. These structures are most simply represented as tree structures, where each node represents a region in the image. The idea with the structures:

the leaves of the structure are regions of a fine partition, or over-segmentation of an image (e.g. pixels, flat zones - connected regions of same intensity, watershed)
the links in the tree represent merging or union of (neighboring) regions and forming a more complex region, and are constructed in a way most likely to produce regions corresponding to objects (hopefully :) )
the complexity of the regions increases along each branch from nodes towards the root of the tree
every level of the tree (also, every cut) is a partitioning/segmentation of an image (coarser the closer to the root it is)
the root of the tree covers the whole image domain

The segmentation than consists of examining the regions and their proposed unions, to determine the regions in the tree or the cut of the tree corresponding to the required precision, or some known properties about object of interest, or other pre-defined user specifications.

Trees (that is, hierarchical image decompositions) with such properties are:

$\alpha$-trees
$(\omega)$-trees
Binary partition trees (aimed at segmentation, and object detection)

In addition to links to papers already provided, some more specific and less specific current papers on these kind of segmentation techniques:

Cousty, J., Najman, L.: Incremental Algorithm for Hierarchical Minimum Spanning Forests and Saliency of Watershed Cuts (2011)
Soille, P., Grazzini, J.: Constrained Connectivity and Transition Regions (2009)

(more practical):

Serna, A., Marcotegui, B.: Attribute Controlled Reconstruction and Adaptive Mathematical Morphology (2013)
Merciol, F., Lefèvre, S.: Fast image and video segmentation based on alpha-tree multiscale representation (2012)


This is a bug, either in your code, or in your version of the Matlab function grpdelay. As you correctly pointed out, making the last denominator coefficient imaginary just rotates the pole-zero diagram. This means that for the given example, the zeros stay where they were (in the origin), and the pole angles are rotated by $\pi/8$; or course, the pole radii remain unchanged.

For the group delay function this means that it only gets shifted (along the frequency axis), but its maximum and minimum values (and its overall shape) must remain the same.

Using Octave I get the expected result:

g1 = grpdelay(1,[1 0 0 0 .1],1000,'whole');
g2 = grpdelay(1,[1 0 0 0 .1j],1000,'whole');
plot([g1(:),g2(:)]), axis([0,1000,-.4,.5]), grid on


This is a linear filter because the signal values ($y_i$ in your notation) enter the equation linearly. It can be implemented as an FIR filter, and then it's called "Slope Filtering". The method was described in the IEEE Magazine on Signal Processing ("DSP Tips & Tricks"), Nov. 2008. You can read the article here.


Zero-padding data for a longer FFT is equivalent to interpolation by a (periodic) Sinc kernel.  Interpolation by a (periodic) Sinc kernel can reconstruct points between samples of a signal that was strictly bandlimited (to below the Nyquist frequency) prior to sampling.  See: https://ccrma.stanford.edu/~jos/resample/Theory_Ideal_Bandlimited_Interpolation.html .  This is a property of bandlimited signals (e.g. a linear-phase "brick-wall" bandlimiting FIR filter looks a lot like a Sinc function).

A lot of real-world data can be made close enough to bandlimited to meet existing or needed S/N criteria, and thus a zero padded FFT is "close enough" to "work" as needed.

It is not necessary to implicitly assume the repetition of a window of FFT data for all uses of an FFT.  The phase vocoder and fast convolution algorithms require the assumption that the data is windowed and not necessarily periodic in FFT aperture width, otherwise one could end up with a contradiction regarding the data samples across overlapped or adjacent windows.


Given an FIR filter coefficients, how can I find coefficients of its equalizer filter?

Assuming that you also want to equalize the filter's phase response (not only its magnitude response), you need an equalizer with a transfer function $E(z)$ that is the inverse of the FIR filter's transfer function $H(z)$:

$$E(z)=\frac{1}{H(z)}=\frac{1}{c_0+c_1z^{-1}+\ldots+c_{N-1}z^{-(N-1)}}\tag{1}$$

This is an all-pole filter which can be implemented by the following recursion:

$$y[n]=\frac{1}{c_0}\left(x[n]-c_1y[n-1]-\ldots-c_{N-1}y[n-N+1]\right)\tag{2}$$

where $x[n]$ is the input sequence, and $y[n]$ is the output sequence.

The problem with implementing Eq. $(2)$ is that the corresponding filter can be unstable. This depends on the locations of the zeros of the FIR filter's transfer function $H(z)$. If $H(z)$ has zeros on or outside the unit circle of the complex plane (i.e., zeros with magnitudes equal to or greater than $1$), then $E(z)$ has poles at those locations, and poles on or outside the unit circle cause instability of the equalizer. This means that if $H(z)$ has zeros with magnitude equal to or greater than $1$ (i.e. if $H(z)$ is not strictly minimum phase), then there exists no causal and stable equalizer.

In practice you can design an approximate equalization filter by choosing the filter coefficients of the equalizer such that the error

$$\epsilon=\sum_i\left|H(e^{j\omega_i})E(e^{j\omega_i})-e^{-j\omega_i D}\right|^2$$

is minimized, where $D$ is some appropriately chosen delay, and the sum is computed on a dense frequency grid. If the equalizer $E(e^{j\omega})$ is chosen to be an FIR filter, then the resulting optimization problem can be solved by solving a set of linear equations.


Count the number of times the phone rings before its answered

If you can identify one phone ring, then you can extract it and use it to construct a matched filter. The beginning of the wiki article is math heavy, but check out the example section. Additionally, Matlab has a page describing the matched filter from one of their toolboxes. You should really give it a try; it works like magic.

In the time domain, the matched filter is a cross-correlation of the signal against a known pattern (which is in this case a recording of the ringing). As a cross correlation it can efficiently be implemented using a DFT. The output of the matched filter is a time series.

The time series output will have a series of peaks whose locations correspond to the beginning of each ring. All you need to do is count the number of peaks which exceed some threshold.

If the rings could be different but you are guaranteed that they're selected from a finite set of rings (imagine different phone carriers or something) then you could use a bank of matched filters.


I agree with A_A. Just take deeper in Butterworth filter design and this article may help you to understand better.


Designing a 2x interpolation filter for audio

@Keith: Here's the situation as I see it. Your original signal's spectrum is shown in panel (a) below. The solid lines are spectral energy and the dotted lines are spectral replications. After the zero-stuffing (upsampling by two), the stuffed sequence's spectrum is shown in panel (b). Notice the spectral image centered at 44.1 kHz that needs to be removed (attenuated). So after zero-stuffing you need to implement a digital filter whose passband cutoff freq is 0.45*22.05 kHz and whose stopband begins at 0.55*22.05 kHz. (I show that filter as the dashed line in panel (b).) But keep in mind, the filter’s Nyquist rate is 44.1 kHz (88.2/2 kHz) and not 22.05 kHz. The final spectrum you want is shown in panel (c).


What exactly does compression say about correlation of data?

What might this indicate viz a viz any possible correlation within the data?  Is it at all possible to convert 2.7 to a lag figure of n bits, even loosely?

The compression ratio is directly proportional to the entropy ($H$) of the source. For example, Shannon's Entropy tells you how many bits on average would be required to encode a source, given its statistical distribution. The expression for Shannon Entropy is:

$$ H(S) = -\sum_{i=0}^{|S|-1} P(s_i) \log_2 P(s_i)$$

Where $|\cdot |$ indicates the size of an ordered set, $S$ is a memoryless source having $|S|$ symbols, each one appearing with probability $P(s_i)$ in a stream.

If the source was totally predictable ($1111111111111111111111111111 \ldots$) then $S_p=[1]$, $P(S_p) = [1]$ and to describe it you would need $H(S_p)=0$ bits.

If the source had to alternate between two states but in a totally unpredictable way then $S_u=[0,1], P(S_u) = [0.5,0.5]$ and to describe it you would need $H(S_u) = -(0.5 \cdot \log_2 0.5) \cdot 2 = 1$ bit.

BUT, if the same source started showing a preference for one of the states, it suddenly became just a bit more predictable. So, let's say $S_b=[0,1], P(S_b)=[0.7,0.3]$ and $H(S) = -(0.7 \cdot \log_2 0.7 + 0.3 \cdot \log_2 0.3 \approx 0.88129...$.

Remember that $H$ is expressed in bits here. If I need $H(S)$ bits per symbol in $S$ and produce a message $M$ that is $N$ symbols long, then I need $|M| = H(S) \cdot N$ bits to describe the message.

The unpredictable source $S_u$ needs 1 bit per symbol, so, we need as many bits as symbols in $M$ and $|M_u| = N_u \cdot H(S_u)$. But the "biased" source $S_b$ needs $\approx 0.88129$ per symbol, so $|M_b| = N_b \cdot H(S_b)$ and theoretically at least, we could fit it in less space.

The question now is, how do you do that? How do you re-encode a sequence to occupy less space? Because that is very, very important (for the interpretation of metrics such as this one).

Without wanting to go into too much detail (now), let's say that we do have a way of encoding a message in such a way that the compression ratio approximates its entropy as the length of the message tends to infinity. This is what you call $C$ (and it is encouraging that you have selected zip, bz2, etc).

Then, to form a ratio of the lengths of the shuffled and unshuffled compressed sequences is equivalent to forming the ratio of $\frac{|M_{u \leftarrow b}|}{|M_b|}$.

In other words, you have a file that is coherent because it bears content. Because of this reason, the bytes appear in some predictable way, following some kind of pattern, kind of like $S_b$ but this time around your $S_b$ contains 256 symbols and your $P(S_b)$ contains 256 probabilities (this is not exactly accurate here because a file is not memoryless but I will come back to this later).

When you shuffle the file, you basically construct an equivalent white noise sequence from a source $S_u$ where $P(S_u) = \frac{1}{256} \forall s \in S_u$ (Where $\forall$ means "for all"). This is what I mean up there with that $M_u \leftarrow M_b$.

Earlier, $S_u,S_b$ were unrelated, but here $S_u$ is the result of $S_b$ after you shuffled it.

So, having said all this, what you basically have expressed there is a ratio that is equivalent to the inverse of the compression ratio of $C$ on your data $d$.

Provided that the shuffling is adequate $|C(ds)|$ is equivalent to the length of $ds$ because it is impossible to compress it. And, the length of $|C(d)|$ is approximated by the entropy of $d$ times its length.

What does this tell you? It tells you how much bigger is a white noise sequence using the same symbols in $d$ compared to the coherent version of $d$.

What does this say about the auto-correlation of the data? Under very strict conditions it can tell you the average size of the longest redundant sequence of bits. If the data compression ratio is 0.5, on average, you are substituting 2 bits for 1.

And this, is where things can go really bad, in terms of choice of $C$ and interpretation.

The compressor $C$, doesn't care what is in $d$. The only thing it sees is a stream of bytes. If you have a file of doubles in a common format, the compressor will compress each and every one of these 8 bytes separately, not as one number. What your ratio will return, is information on the efficiency of compressing the stream, not your "data".

  Is it at all possible to convert 2.7 to a lag figure of n bits, even loosely?

This was answered above, but this discussion so far brings us to the wonderful subject of LZ Complexity.

The LZW algorithm is a dictionary based algorithm. It compresses a sequence by keeping a running dictionary of fragments it has encountered before.

Consider 0101001010010. We start at the first 0, have we seen it before? No, put it in the dictionary, output 0, move to the next. Have we seen 1 before?, No, put it in the dictionary, output 1, move to the next. Have we seen 0 before? Yes, carry on to the next. Have we seen 01 before? No, put it in the dictionary, output the index of the known sequence (the 0) and the new bit (1) and move to the next....and so on until the end of the file.

You see what's happening here, the algorithm compresses when it encounters a sequence it has seen before and it substitutes it for its index in the dictionary. In these early stages it does not seem to compress much but if the source is predictable then the algorithm will build a set of increasingly long repeating sequences very quickly and it will be able to start substituting repetitions. Have we seen 001100110011001100110011 before? Yes, it is entry 242. So, you substitute a long string of bytes for an integer.

The reason I am talking about this is because if you wanted to make such strong statements about the lag of the autocorrelation, it is not enough to examine the compression ratio. Ideally, you would have to examine the patterns within the dictionaries.

If you are lucky, your data is inherently digital (categorical) and you might be able to do this. But if your data is real numbers, then you have to quantise them and the way you do this makes it very difficult to apply these methods to real valued data.

Hope this helps.


Phase delay between two signal using FFT

Do an fftshift before the FFT (rotate your sample vector by N/2), this relocates the phase reference to the center of your sample window (instead of at a potential circular discontinuity around the ends).  This also often requires less unwrapping, as well as making phase interpolation (between FFT result bins) and phase differences more intuitive.


Measure how Straight/Smooth the Borders are Rendered in an Image

I'd try a very "tinkery" approach here:

Erode the image, so that the black area is shrunk by a fixed radius of pixels from its border (say, 5px).
Dilate the resulting image by the same amount
measure the amount of difference between original and processed image.

The idea is that something that is a locally convex border doesn't suffer through erosion (it's only shrunk) significantly, and that this erosion can be reverted by dilatation.


Suppose one signal is $a_1(t) = \cos(\omega t)$ and the other is
$a_2(t) = \cos(\omega t - \theta)$. Note the
phase shift $\theta$ between the two signals, and note also that we can write
$$a_2(t) = \cos(\omega t - \theta)
= \cos\left(\omega\left(t-\frac{\theta}{\omega}\right)\right)
= \cos(\omega (t - t_0)) = a_1(t-t_0)$$
where $t_0 = \frac{\theta}{\omega}$ showing that $a_2(t)$ is a delayed version of $a_1(t)$.  Suppose that $a_1(t)$ and $a_2(t)$
are transmitted from the two separate antennas simultaneously, that is,
no additional time delay between them,
and reach the receiver over paths of delays $\tau_1$ and $\tau_2$
respectively and gains $A_1$ and $A_2$ respectively.
Then the received signals are
$$\begin{align}
A_1a_1(t-\tau_1) &= A_1\cos(\omega (t -\tau_1)) = A_1\cos(\omega \hat{t})\\
A_2a_2(t-\tau_2) &= A_2\cos(\omega (t -t_0 - \tau_2))\\
&= A_2\cos(\omega (\hat{t} +\tau_1 -t_0-\tau_2))\\
&= A_2\cos(\omega (\hat{t} - (t_0 + \tau_2-\tau_1))\\
&= A_2\cos(\omega \hat{t} -\hat{\theta})
\end{align}$$
where $\hat{t} = t-\tau_1$ is the time as measured by the receiver clock.
Thus we see that the second signal is delayed
with respect to the former by $t_0+\tau_2-\tau_1$, that is,
the delay introduced at the transmitter plus the difference in the
delays between the two paths, and this can be expressed as a phase
shift of $\hat{\theta}= \frac{t_0+\tau_2-\tau_1}{\omega}$.

If $\tau_1 = \tau_2$ so that the
two paths are of identical lengths, the second signal
looks like a multipath signal with respect to the first with a delay
of $t_0$ (equivalently the same phase shift $\theta$ as existed at
the transmitter).
If $t_0 = 0$ so that effectively only one signal is
transmitted, but is received over paths of different lengths,
the second signal is a multipath signal with respect to the first
with a delay of $\tau_2-\tau_1$ and it too can be expressed
as a phase shifted version of the first signal.


I think all these answers fail to provide a descent methodology.

You should consider "the closeness" of the images as the availability of the initial pose and cast this as a registration problem. The solution can typically be obtained through a least-squares adjustment, i.e., by minimizing the distances of the model points (from first image) to their corresponding image points. In case of a perspective distortion, you could optimize the 9 parameters of the homography matrix. Methods like Lukas Kanade, ESM or hyperplane tracking are some ways to do that. In essence, the procedure is similar to ICP, where an iterative procedure is applied.

Here is one work, extending ESM and comparing a few methods. If your camera is calibrated, you could benefit from some alternative methods.

Note that, to begin with, one could always resort to a simple approach, where you try different combinations of rigid motion $[-2,2]$ degrees for rotation and $[-5,5]$ pixels for translation. Then transform your template with all possible poses (in this case ~$40-50$). Then interpolate the sum of squared differences or normalized cross correlation scores in the pose space (6DOF freedom or 9DOF homography) to get a subpixel accurate refinement. This should give you a rough yet a better alignment. Keep in mind that you could also use this scheme iteratively.

If you have absolutely no clue on how to do any of these, check out this link.


Frequency spectrum of a set of images in Matlab

If I'm interpreting the question correctly, you're looking for a way to find a distribution of the frequencies in each images. It sounds like you might need some sort of 3D scatterplot.

If you have 20 images that are 256x256, how I imagine this working is,

Take your 2D FFT or 2D DCT (or blockwise FFT/DCT) for each image
Normalize all the values so you can have a good comparison between each image
Concatenate all the frequency domain images into a 256 x 256 x 20 matrix
Filter out the values that are below a certain threshold because they do not represent a significant component of your image. For example, if the magnitude of the values in the corners that represent high frequency information are very low, then removing them won't make a huge difference to your image because most of the information is elsewhere. But removing those data points create holes in your 3D plot so you can see the information that you're interested in
Plot the data in 3D. If you're using Matlab, you can use scatter3(). The plot will be easier to read if you apply a heatmap.

Otherwise, if you're looking for the frequency information of a single image, there's several ways to represent the data depending on your problem. If you clarify the question you have and it's application, I can edit this answer.

Histogram

If you take the histogram of a grayscale image, you'll see a distribution of the intensity values of each pixel.

This easiest way to see this is to open up the image in a photo-editing software, like Photoshop, GIMP, Picasa and select histogram from the menu.

If you want to use Matlab, the hist() function plots a histogram similar to the one below. You can directly pass your image variable into the hist() function. But it's important to pick a good bin size, because different bin sizes can dramatically change the representation of the data. There are ways to calculate optimal bin sizes, but you can always guess and check.

2D FFT

Follow the Wikipedia article for the theoretical explanation, but basically the 2D FFT gives you an image of the intensity of low frequency to high frequency information. The center of the image is low frequency information (such as flat surface, sky, skin, walls) and the corners are high frequency information (edges, grainy noise, intricate patterns).

If you're using Matlab, you need to shift the 2D FFT to put the low frequency information in the center of the image. It's also a good idea to normalize the 2D FFT so you can see the details better; it's similar to contrast enhancing a photograph to see the details.

2D DCT

Follow the Wikipedia article for the theoretical explanation. The 2D DCT is very similar to the 2D FFT; It also gives you an image of the intensity of low frequency to high frequency information, but the low frequency information is in the top-left corner and the high frequency information is in the bottom-right corner. This makes the information easier to manipulate (at least for me).

Generally, the 2D DCT is used on blocks of an image (8x8 pixels). This gives you frequency information that is still related to its location in the spatial image.

The 2D DCT is useful for encoding images. A variant of JPEG encoding uses the 2D DCT transform.


Cross-Correlation Signal Delay Estimation Variance

One warning first: this is going to be a bit literature-heavy, because I think that's what you're more or less explicitly ask for. You want something that you can read, understand and cite for your publication.

  Firstly, my understanding is that a matched filter is the most common method for estimating signal-delay.

No. A matched filter can be used to effectively correlate a receive signal against a transmit signal shape. It's can but be a puzzle piece in a TDOA system.

  Thus, what I'm looking for is how to calculate the variance in the signal-delay estimate using a matched filter as an estimator.

Again, the matched filter is not the delay estimator; think of it as a signal preconditioner!

Now, what you probably mean is that you use matched filters for the TX pulse shape and the RX pulse shape filter, and then, based on the post-filter signal, do a delay estimate.

I'd like to cite a paper,[1] :

which doesn't actually deal with matched filtering directly, as the TDoA system isn't able to define the pulse shape the transmitter uses; instead, it's about cross-correlation of the received signals at different positions.

The interesting point is that the defining entity (aside from the obvious $\frac{c}{f_\text{Nyquist}}$ limit) is the SNR of the received signal -- that's pretty obvious, if you think about it: The "certainty" with which the correlation of RX and TX signal has a maximum at the "right" position is reduced with increased SNR.

In previous work of the author, she uses in [2] that with AWGN assumption and cross-correlation (i.e. matched filter with a known TX signal)

$$\begin{align*}
\sigma^2_{\hat\tau}&=\frac{3}{8\pi^2}\,\frac{1+2\, \text{SNR}}{\text{SNR}^2}\, \frac {1}{T_o B^3}\text,\,\text{with}\\
\sigma^2_{\hat\tau}&: \text{Variance of the timing estimate}\\
T_o&: \text{Duration of observation, and}\\
B&: \text{bandwidth of the system.}
\end{align*}$$

Interesting here, aside from the $B^{-3}$ relationship of variance and bandwidth, is the cited paper (if you find it in fulltext available online for free; I have IEEExplore access, but I don't know about you...) [9], which pretty much should give you plenty of background to work with.

[1]Noha El Gemayel, Holger Jäkel, Friedrich K. Jondral: Error Analysis of a Low Cost TDoA Sensor Network, IEEE/ION Position, Location and Navigation Symposium 2014, Monterey, CA, USA, May 5-8, 2014

[2]Noha El Gemayel, Sebastian Koslowski, Friedrich K. Jondral:
A low cost TDOA Localization System: Setup, Challenges and Results,
,
10th Workshop on Positioning, Navigation and Communication 2013 (WPNC'13), Dresden, Germany, March 20-21, 2013

[9] G. Carter, Coherence and time delay estimation, Proceedings of the
IEEE, vol. 75, no. 2, pp. 236 – 255, Feb. 1987


spectrogram(x, N, N/2, N, fs, ‘yaxis’) is my default one liner.

To get rid of the non-SI time units, I think you need to do
[S, F, T] = spectrogram(x, N, N/2, N, fs);
pcolor(S, F, 20*log10(abs(S)));
cb = colorbar()
cb.Label.String = «some unit»;
cl = get(gca, ‘CLim’);
set(gca, ‘CLim’, cl(2)+[-40 0];

From memory. You might want to check that S is «amplitude» rather than power, and if pcolor takes S and F in this order or backwards. Also consider shading flat (if dataset is dense) or using imagesc() instead of pcolor.


Design narrow bandpass filter for signal with high sampling rate

So, first, to put things into perspective: 4kHz is not a high sampling rate these days (add 5 orders of magnitudes, and things become hard).

Your 15 kHz passband doesn't say anything about the complexity of the filter; what counts is the transition width, ie. the distance between pass- and stopband, as well as the attenuation of the stopband. For a slightly specialized answer, see the answers to http://dsp.stackexchange.com/questions/31066/how-many-taps-does-an-fir-filter-need/31077.

  What is the recommended way to perform tasks like this?

Do a filter design, apply filter.

  What bandpass filter is suitable?

Any bandpass filter that suits your needs – which you haven't specified.

What's missing is:

transition width
stop-band attenuation
acceptable ripple
length constraints

What we know is:

You want a linear phase filter (because linear phase means constant group delay, and that means un-broken phase relationships)
Thus, you probably want a FIR that is symmetric in time

  Furthermore, the filter step response should have close to zero overshoot.

Gibb's phenomenon is non-negotiable :) so yeah, use a longish filter with a nice rolloff.

  The settling time could therefore be a little larger.

That will be the effect, yes. You could use a filter design tool that allows you to design with a windowing method and use a window that suits your application well – which I don't know.

  I am more interested in a general solution that delivers solid results using standard IIR or FIR filters (e. g. as available in Python's SciPy library)

Well, yeah, that's what I'd recommend, but you say just shortly above that you don't want to do a filter design? I'm a bit conflicted. Anyway, there's a lot of functions that will give you a proper design.

I'd recommend the following:

Design a low-pass FIR filter
convert it to a bandpass filter, by multiplying with a cosine (real-valued data, symmetric filter) or $e^{j2\pi\frac{f_{center}}{f_s}\cdot n}$ (complex-valued data, one-sided filter) as needed.

Designing is easy, something along (untested, straight from the back of my head)

from scipy.signal import fir_filter_design
from math import cos, pi

f_center=ToBeDefined!!!
f_cut = 15
f_s = 4e3
f_rel = f_s/2/f_cut
transition_width = ToBeDefined!!! # e.g. 10
trans_rel = f_s/2/transition_width
attenuation = ToBeDefined!!! # e.g. 10.0**-5

num_of_taps = estimate_by_reading_my_answer_above(f_rel, trans_rel, attenuation)

lowpass_taps = fir_filter_design.firwin2(numtaps=num_of_taps, cutoff=f_rel, window = "hamming"|"blackmanharris"|"hann"|"chebwin"|…)

bandpass_taps = [lp_tap * cos(2*pi*f_center/f_s*n) for n, lptap in enumerate(lowpass_taps)]

  Do I need to downsample the signal before filtering. If yes, what type of additional low pass filter should I use to avoid aliasing?

Well, considering you won't need most of your signal, yes, that seems advantageous. Typically, you'd just filter away (and decimate in the same step) an 1/N of your initial sampling rate. For example, if your bandpass lies about 100 Hz - 15 / 2 Hz to 100 Hz + 15 / 2 Hz, you won't need anything above 200 Hz – so just decimate to 1/20 of your input sampling rate, with a low pass filter that has 1/40 input sampling rate transition width.

But again, 4kHz is laughably little, and a couple thousand taps won't stress your PC at all. If in doubt,

from scipy import signal
signal.convolve(your_input_signal, filter_taps, method="fft")

does a fast convolution. (you can also just omit the method, or set it to auto, because scipy will just pick fast convolution for long filters, automatically). For a bit of an impression on what's possible with halfway-optimized code on a PC, see my answer to how to implement signal generation on a GPU (and why not); there, I showed that you can do a 107-tap filter at 160 MS/s on my PC back then, so you should be able to do a let's say 700 tap filter at easily 20 MS/s (just a wild guess). That's only 5 thousand times as fast as your sampling.


One does need to (not a necessary condition)! But if you do, you can recover (at least theoretically) the continuous signal from the discrete samples (sufficient condition).

The bandwidth most usually  denotes the span between the zeroth and the maximum frequency (in the analog signal). Sometimes, for a band-pass signal,  it denotes the span between the minimum (positive) and the maximum frequency, and there are similar theorems: twice the bandwidth (sometimes a little more depending on the location of the minimum or the maximum) is enough.

If you sample below, aliasing might occur, but not every-time, this depends on the structure of the spectrum.

You care only about 1 Hz, it is dangerous to simply resample: it is advised to filter a little above 1 Hz (say 1.5) and then resample consistently with 1.5 and the precision of your filter. Could be 4 Hz (above $2*1.5$ Hz). Those figures should not be taken for granted, just an idea of the whole method.

Why what it set to 100 kHz initially? Several possibilities, for instance:

people did not know initially in which part of the spectrum the information was located.
the acquisition scheme  was designed in a "the more data you have/the more information" mood


First, it seems to me that you are doing the analysis incorrectly. 128 linearly spaced channels is both too high (it won't capture a smoothed spectral envelope - vocoders usually have at most ~30 channels) and too small (it has a very poor resolution in the lowest frequencies and won't discriminate different formants). You really need to replace that by a bank of constant-Q band-pass filters, 2 or 3 bands / octave to get something that is coarse enough to capture the spectral envelope without picking individual frequency peaks, yet discriminates the frequencies of interest.

Your analysis and synthesis filters can be implemented with biquads. Running a dozen of those in parallel is not cost-prohibitive compared to a FFT.

As for the carrier wave, you can start with a sawtooth or narrow pulse - avoid waveforms like triangle, square, or sine whose spectrum have "gaps". Then stack a few of them to play a chord. You can thicken the carrier by adding slightly detuned version of each component or octaves. Note that in order to achieve the best sounding results, band-limited synthesis (for example with minBLEPs) has to be used to synthesize the square or sawtooth waves.


All-pass filter usage in automotive room correction

Allpass filters are very helpful in automotive mainly for getting the imaging dialed in. You often need to time align the contributions of different driver groups (left, center , right) at a specific listening position. Due to the complicated acoustics of the cabin. the delays are often dependent on frequency and so all passes can be very helpful in correcting this.
Allpasses in cross overs are typically not that useful, since the range of overlapping frequencies is pretty narrow and you can typically get by with a single delay


Assuming that you denote the output by $y(t)$ and the input by $x(t)$, the initial conditions are given by the values of $y(t)$ and its derivatives at a certain time $t_0$ (choose $t_0=0^-$ if you like). A value $x(t_0)$ has nothing to do with an initial condition.

For the zero-input response (ZIR) you just set $x(t)=0$ and compute the output caused by possibly non-zero initial conditions. For the zero-state response (ZSR), you set $y^{(k)}(t_0)=0$, $k=0,1,\ldots$, and compute the output for a given input $x(t)$. The total response is then the sum of the ZIR and the ZSR.


What would be the best way to estimate velocity from a video?

You can probably address that by tracking debris across frames, assuming that they move at the speed of water and aren't too far apart in successive frames.

Maybe start with the old stuff: Shi, J.; Tomasi, C., "Good features to track," in Computer Vision and Pattern Recognition, 1994. Proceedings CVPR '94., 1994 IEEE Computer Society Conference on , vol., no., pp.593-600, 21-23 Jun 1994.


Is uniform noise present in some real processes?

Typically, when a real value $x$ from a sensor is discretized ($x_Q$) with a uniform quantizer, the error $x-x_Q$ is likely to be uniformly distributed.

Caveat: as commented by @bone, this can be stated in precise terms, see for instance:

A necessary and sufficient condition for quantization errors to be uniform and white,A. Sripad and D. Snyder, 1977

  In this paper, a necessary and sufficient condition is given to model
  the output of a quantizer as an infinite-precision input and an
  additive, uniform, white noise. The statistical properties of the
  quantization error are studied, and a detailed analysis for Gaussian
  distributed inputs is given.

Statistical Theory of Quantization, Bernard Widrow, Istvan Kollar,  Ming-Chang Liu, 1996

  The effect of uniform quantization can often be modeled by an additive
  noise that is uniformly distributed, uncorrelated with the input
  signal, and has a white spectrum. This paper surveys the theory behind
  this model, and discusses the conditions of its validity. The
  application of the model to floating-point quantization is
  demonstrated.


Necessary image processing operations before template matching

The sample image you posted has relatively strong perspective (it is not imaged straight from the direction of surface normal) which can cause problems with template matching techniques witch use block processing. I assume that you have to take the image with strong perspective so first thing we want to do is estimate image transformation which will remove the perspective projection and will result rectified image ("straight image").

Correcting perspective projection when shape is known

The ultimate goal of the preprocessing is to find corners of that gray rectangle which is the background the color chart. I started by simply thresholding with constant value. Literature/web has lots of information on how to do thresholding when the image has uneven lighting, but if you can, try to correct the imaging setup such that the lighting is relatively even.

I assumed that the color chart takes quite large portion of the screen, so the it is probably biggest continuous region of the image. There exist many black and white labeling functions (MATLABs bwlabel, IPPs LabelMarkers, OpenCV has external library cvblob) which will give each distinct connected region its own index. After the labeling you can calculate the areas of the connected components using a histogram and then choose the component with largest area.

It is usually good idea to fill the components such that component does not have holes, this reduces the number of edge pixels.

Now we can use edge finding operators (I used gradient based method, but you can also trace the edge by following white pixels that are connected to black pixels) to get pixels belonging to the boarder.

As the object we are interested is rectangle, I fitted four lines to the edge pixels using robust line fitting (RANSAC). After finding new line, I removed all the points that were close to the found line, this forces subsequent fitting operations to return other edges of the object. Finally when all the edges were found, I computed all possible intersections between the set of lines. From the set of resulting points I removed all points outside of the image.

Now we have the corners of the color chart and as we know the dimensions of the chart, we can calculate the transformation matrix $H$ (homography, projection between two planes. use direct linear transformation (DLT) to solve $H$) between the millimeter coordinates and the pixel coordinates.

$\textbf{x} = H\textbf{X}$

This transformation matrix can be then used to transform the original image such that color chart is imaged straight from the direction of the surface normal.

As can be seen only radial distortion is left. When we used transformation $H$ for image transformation (2D interpolation), we chose the interpolation coordinates such, that the we know the precise resolution of the image, meaning that we now know for example the size of the letters.

If we still want to do some template matching, we can use any decent template matching method. The method does not have to be rotation/scale invariant as the image is already rectified.


Covariance matrix of an adaptive filter input

The Covariance Matrix is commonly defined as

$$\mathbf Q = E\left[ (\mathbf x -\mathbf\mu_{x})(\mathbf x -\mathbf\mu_{x})^*\right]$$
with $\mu$ denoting the mean value, i.e. $\mu_{x}=E\left[\mathbf x\right]$, and $\mathbf x$ being column vectors. The fact that you define the covariance matrix as

$$\mathbf{R}_i = E\left[\textbf{u}_i^*\textbf{u}_i \right]$$

indicates that your $\mathbf u_i$ have zero mean and are row vectors; that might be a result of the signal model you use, but it's pretty non-standard.

Furthermore, the product of a column vector with its hermitian is inherently hermitian (symmetric if real).

Now, one way to estimate an expectation value from multiple observations. Wikipedia has an article on estimating covariance matrices, so here's just the gist: the most intuitive estimator is the sample covariance matrix; adapted to your formula, that'd be

$$\hat{\mathbf R} = \frac1{N-1}\sum\limits_{i=1}^N \mathbf u_i^*\mathbf u_i$$


The computation procedure problem of calculating the BER

This is pure mathematics. However, these are some hints.

1) In 2-D coordinates $t$ and $u$ the multiple integral scans the region limited by $t = f_u(u) = 0$ and $t = g_u(u) = \sqrt{2u\textrm{SNR}}$ for $t > 0$ and $u >0$. It is the same region limited by two functions $u = f_t(t) = 0$ and $u = g_t(t) = t^2/2\textrm{SNR}$ for $t > 0$ and $u >0$.

2) because $\int_a^b e^{-u} du = e^{-a} - e^{-b}$

3) use 2) with the identity Q(0) = 1/2.


I would guess that this is the last point in the frequency response with a passband gain that is still inside the ripple spec. Meaning if you use the 3-dB bandwidth, you will have an attenuation of 3 dB (and therefore a "ripple" of at least 3 dB assuming the ripple spec for the digital filter is less than 3 dB). On the other hand, if you use the ripple limit bandwidth, your attenuation at that point will be within whatever the specified ripple is.

I couldn't find the spec for the ripple in the data sheet, but consider the following example. If the ripple were 0.1 dB and the passband gain were 0 dB then the attenuation at the ripple limit point would be 0.05 dB and from there the gain would decrease until reaching the 3 dB point where the attenuation would be 3 dB.

There is a chance I am wrong about this, since it isn't a standard term I'm familiar with, but I'm posting the answer because it seems very likely to be correct. Maybe someone else with more experience can provide a better response or point out my error if there is one.

The reason they probably include both is that a digital filter designer would only be interested in the ripple limit bandwidth (which we would generally refer to as the passband edge when doing an equiripple design) while many analog filters are specified with their 3-dB (or 6-dB) bandwidth. Probably there was some kind of internal discussion about which one to use, and they wound up with both. But that's all speculation on my part.


Cross power spectral density (CPSD) comparison between signals

I've ended up at a workable solution.

To deal with the noise, average the CPSD results of many realizations of the noise.

Pxy_AvB_mat = zeros(512/2+1,nRealizations); %preallocate
Pxy_BvC_mat = zeros(512/2+1,nRealizations); %preallocate
for i=1:nRealizations
    B = 1*(rand(length(t),1)-0.5); %new realization
    [Cxy_AvB,freqs_AvB] = cpsd(A,B,window,100,512,Fs);
    Pxy_AvB_mat(:,i) = abs(Cxy_AvB);
    [Cxy_BvC,freqs_BvC] = cpsd(B,C,window,100,512,Fs);
    Pxy_BvC_mat(:,i) = abs(Cxy_BvC);
end
Pxy_AvB = mean(Pxy_AvB_mat,2);
Pxy_BvC = mean(Pxy_BvC_mat,2);

Use a spectral analysis tool like periodogram, pwelch, or roll your own with FFTs to verify if a signal is pure noise or not. Averaging multiple realizations of noise's spectral powers will yield a flat line - showing there are no major noise constituents. Same idea as the code above. A single realization can be misleading for spectral analysis or CPSD. Many-realizations-averaged-CPSD will still show some minor peaks with respect to the noise signal, but the spectral analysis in combination will help confirm if there are any actual frequency components (noise on its own is flat).


Break a signal into segments, add them together and then perform an FFT on the result

That's not true (unless the signal is periodic, and those pieces are the base periods of it).

If you divide a signal $x_L[n]$ of a length $L$ into smaller pieces $x_i[n]$ of length $N$ each,

$$ x_L[n] = x_1[n] + x_2[n-N] + x_3[n-2N] + x_4[n-3N]$$

and sum them into a shorter signal $x_N[n]$ of length N

$$ x_N[n] = x_1[n] + x_2[n] + x_3[n] + x_4[n]$$

and take $N$-point DFT/FFT of their resultant $x_N[n]$, you won't be able to reconstruct the original long singnal $x_L[n]$ back.

Indeed it's not even about FFT. You would not be able to retrieve the components $x_i[n]$ from their sum $x_N[n]$. So since there is no way to reconstruct back $x_L[n]$ in this manner.


Suppose that you could write $x_.$ as $x_.(t) = a_1x_1(t-n_1)+a_2x_1(t-n_2)$, then by LTI hypothesis, you can derive $y_.(t) = a_1y_1(t-n_1)+a_2y_1(t-n_2)$ with very little computation.
You can get a first insight to that method  by looking at the shape of $x_2$ which is self-evident, and a little more imagination gives you the structure for  $x_3$.
I find that this exercise is very pedagogic in the progression of difficulties. Some other exercises of the book might have the same gradation, don't hesitate to look back before jumping into computations at first.


Linearity, Causality and Stability of a System

According to the interpratation (that belongs to you) that output $y[n]$ is given by the recursion relation $$y[n] = y[n-1] + u[n]$$ no matter what the input $x[n]$ is, then the following conclusions can be made about the linearity, time-invariance, stability and causality of the system.

First let's solve the output for $y[n]$

i.For $n<0$, since $u[n]=0$, we have $$y[n] = y[n-1]$$ which implies that $y[-1]=y[-2]=...y[-N]...= K$ for all $n < 0$

ii.For $n\geq0$, since $u[n]=1$, then we have
$$y[n] = y[n-1] + 1$$ which can be recursively solved to get

\begin{align}
y[0] &= y[-1] + 1 = K +1 \\
y[1] &= y[0] + 1 = (K+1) +1 = K+2 \\
y[2] &= y[1] + 1 = K +3 \\
&...\\
y[n] &= y[n-1] + 1 = K +(n+1) \\
\end{align}

Combining them yields: $$\boxed{y[n] = K + (n+1) u[n] }$$ for all $n$ as th fixed-unique output of the system for all inputs $x[n]$

Then we analyse the system for its properties:

The system is non-linear:
It's so evident since the system has a fixed, unique output for every input
it's a non-linear system as it cannot satisfy $$\mathcal{T}\{a x_1[n]+ b x_2[n] \} = a \mathcal{T}\{x_1[n] \} + b\mathcal{T}\{x_2[n]\} $$
The system is time-varying:
Again shifting the input have no effect on the output (which is fixed) so it cannot satisfy $$ y[n-d] = \mathcal{T}\{x[n-d] \}$$
NOTE: Does this make any sense? The system is already a fixed one. It has a fixed output for every input. Hence effectively nothing varies in its output. Nevertheless I called it as time-varying, because evidently it cannot saitsfy the time-invariance constraint.
The system is causal: it only depends on the past value of the output and it does not depend on any value of the input so it does not violate causality.
The system is unstable: as for a bounded input $x[n] = \delta[n]$ the output is $y[n] = K + (n+1) u[n]$ which grows unbounded as $n$ goes to infinity.
The system is not invertible.
The system can be implemented with and without any memory requirement.


Phase (or carrier) Recovery for BPSK can be done over the entire sequence using the information from every sample.  Here are common approaches to doing Carrier Recovery:

Frequency Doubling (squaring): If you square a BPSK signal (multiply it by itself) a strong tone will be created at twice your carrier frequency. The Squaring operation strips the data modulation; since the data is modulating the phase back and forth 0 to 180°, and since squaring doubles the frequency, and doubling the frequency doubles the phase, then the phase after doubling will modulate back and forth 0 to 360°, hence the modulation is removed. A PLL is used to track the new unmodulated signal at twice the carrier frequency, which is effectively a self tuning narrow band filter, and creates a copy of the signal with better noise performance. This is followed by a frequency divide by 2 and the carrier is recovered.

Cross product detectors: Regarding other approaches to carrier recovery, that are applicable to BPSK, QPSK and QAM in general, see the block diagram below showing the typical components of a carrier recovery loop:

This loop works in conjuction with a timing recovery loop, and most timing recovery loops can establish timing with a carrier offset, so assume that the timing recovery function has established proper timing and the carrier recovery is working off of the pre-decision best estimate of each symbol. For more info on timing Loops see  Location of Matched Filter and Gardner Timing Recovery for Repeated Sybmols.

Without going into deeper details of proper digital control loop design (which is necessary to understand), I will provide here details on implementations of the phase detector for use in carrier recovery.

First, importantly, the phase between two complex samples is directly proportional to the imaginary portion of the complex conjugate multiplication of those two samples. (Specifically it is proportional to the sine of the phase, but for small angles $\textrm{sin}(\theta ) \approx \theta$ ).  That is a mouthful! Luckily there is a simple and intuitive explanation that I have explained in detail in this post: How to correct the phase offset for QPSK I-Q data, but here is the direct mathematical explanation with reference to the upper plot below:

$Ae^{j(\omega t + \phi)} [Ae^{(j\omega t)}]^*= Ae^{j(\omega t + \phi)} Ae^{(-j\omega t)} = A^2e^{j\phi}= A^2cos\phi + jA^2sin\phi$

With the imaginary part equal to our output of interest: $A^2 sin\phi$

As described in the link, for two complex samples $I_1+jQ_1$ and $I_2+jQ_2$, and treating the $A^2$ as a proportionality constant, this becomes the cross product:

$A^2 sin\phi = I_2Q_1-I_1Q_2$

For small angles, phase is proportional to $I_2Q_1-I_1Q_2$, otherwise known as the cross product. Wow, it doesn't get much simpler than that!  (And for larger angles, all the way from $+\pi$ to $-\pi$ we can still use this phase detector, just be aware of the changing gain in the discriminator when designing the loop; mostly affecting acquisition time.)

See the upper plot below showing this, and first make note that we get a fairly linear phase discriminator output for small angles centered about 0° phase; further this discriminator will provide non-ambiguous phase information (although non-linear due to sine function) for angles from $+\pi$ to $-\pi$.

The lower plot shows the similar relationship in using a multiplier as a phase detector for real signals. In this case the output is proportional to the cosine of an angle, and if used in a carrier recovery loop, the loop will lock to a 90° phase offset since the 0 crossing of the discriminator occurs at 90°. In this post I am focusing on the complex implementation that would be more prevalent in an all digital design.

The question then in the specifics is what are the two signals that you use to compare phase? Here there are several possible approaches, but a common theme is you need to remove the data that is causing the phase to switch back and forth 180°- this can be done with the frequency doubling approach, comparing the doubled (squared) modulated signal with a local reference at twice the carrier, or a decision directed approach where you compare the received signal at the proper timing location to the closest decision of where that symbol should be. Since in this case it is comparing the signal to itself, the modulation effect is naturally removed.

There are also implementations that end up being a Frequency Lock Loop (so will lock in phase but to an arbitrary phase offset) that work by comparing the phase of the current symbol to the phase of the previous symbol (with data removed), or even the phase between samples within a symbol. This measures a change in phase versus a change in time, which is by definition frequency. This approach uses the imaginary of the complex conjugate method, and is often referred to as a cross product frequency discriminator.

Please refer to How to correct the phase offset for QPSK I-Q data for more details on using this concept for phase detector implementations.

This decision directed approach is basically the Costas Loop implementation; see below the implementation that works for BPSK and QPSK (adding more decision boundaries extends this to QAM). Notice that this is simply the implementation of $I_2Q_1-I_1Q_2$ described above, with one of the I, Q samples being the pre-decision sample, and the other being the post-decision sample:

And for context I show below a complete all-digital carrier recovery loop implementation for BPSK, QPSK and QAM, without the design details; the complexity is more in the Loop Filter and how to set the coefficients, other than that the block diagram shown is rather straight-forward: the accumulator followed by the two look-up tables is simply a complex NCO (numerically controlled oscillator), and the four multipliers and adders is a full complex multiplier implementation. The assumption is that earlier in the system the real signal is down-converted from RF and sampled with I and Q components (in either the analog or digital domain).  It is also feasible that the input be real (all on I and Q set to 0, and the digital multiplier in its down-conversion process by multiplying with a complex NCO will create the quadrature I and Q channels needed).


Affine Transformation for matching 2D point sets

There are two ways by which you can approach this problem:

Match the point clouds directly with RANSAC.
Consider matching the resulting graphs (closely related to 1, still)

What you have done already is closest to #1. To complete it, you need to repeatedly obtain subsamples, find the transformation and then check its "error". The "error" is how many projected points are within a distance threshold (which you are going to have to estimate). For more information, please see this thesis. It does a great job at defining the problem in the beginning (which is very important) and then talking about the approaches.

But, since you are already dealing with pose estimation, you can incorporate its constraints into the model and go for estimation via graph matching. In this case, you would be looking at extracting a graph from your images and then compare it with the pose graph, OR, compare the extent to which two graphs coincide (the graph obtained from the left image to the graph obtained from the right image). A key part of this is the derivation of image skeletons. Gonzalez' book has a very good section on skeletons (chapter 11) and their applications.

The existence of a graph imposes certain constraints on the relationships between the points. For example, the hand "leafs" are expected to be no more than 5-7 "hops" from the root node (e.g. the head node) and this can improve the detection rate. For more information please see this paper. Or more generally look for "graph based pose estimation" and "graph recognition". Be aware that "graph matching" as a term is associated with something entirely different in mathematics.

Hope this helps.


Understanding Integration Scale in Multi-Scale Harris Corner Detector

I think the first smoothing, by $\sigma_D$, is only done to get more stable derivatives whereas in the second step the convolution by a Gaussian with $\sigma_I$ is done to establish the 'scale-space' in which the operator is applied.

Ignoring $\sigma_D$, this looks like this in Matlab:

dx = [-1 0 1; -1 0 1; -1 0 1];  % Simple mask for derivative
Ix = conv2(im, dx, 'same');     % Convolve against image
Iy = conv2(im, dx', 'same');    % Again for y-direction with transposed mask

% Gaussian filter, as defined by Kovesi
% www.csse.uwa.edu.au/~pk/Research/MatlabFns/index.html
sigma = 1;
g = fspecial('gaussian', max(1, fix(6*sigma)), sigma);

% Determine smoothed squared image derivatives via convolution
Ix2 = conv2(Ix.^2, g, 'same');
Iy2 = conv2(Iy.^2, g, 'same');
Ixy = conv2(Ix.*Iy, g, 'same');

% Metric to define the corner score
cim = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2;    % Original measure
cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps);    % Kovesi measure


I think that is because your exponentials are wrong.
There is no variable named w in your code, although you refer to it, so I believe it should be something like:

per2(i)=sum(x(k)*exp(-2j*pi*k*w(i))));

But consider using fft instead:

pe2 = abs(fft(x)).^2;

or

pe2 = abs(fft(x(1:K))).^2;

if you insist on using only $K$ samples for some reason (why generate 10000 sample then??).

Perhaps you should consider the Welch method, with length 256

pe2 = pwelch(x,256,128,256);

BTW - you can generate AR process using filter:

x = filter(1,a,wn);

where a is the vector of AR coefficients.


What are the characteristics of a "good" smoothing convolution kernel?

Yes, in general, your #2 is correct.  That being said, both of the filters stink (with your triangle filter being a little better).

No, f3 does not have the same frequency response as f4.  To get an idea of why that is so, you generally have to zero-pad the impulse response before DFT'ing it to get a reasonable idea of what its frequency response looks like.


The other method would be to capture more data in time out to the length desired. This approach would improve the actual frequency resolution achieved, where here frequency resolution means the ability to distinguish in frequency the existence of two closely spaced tones. Zero padding decreasing the bin spacing, meaning how many Hz each bin represents but does NOT increase the resolution: two closely spaced tones would still appear as one tone (spread over multiple bins) if the frequency resolution was not sufficient and the length was only increased by zero-padding.

In general, without further windowing, the frequency resolution in Hz is related to the time duration of the signal according to:

$$f_{BW} \approx 1/T$$

Where $f_{BW}$ is the equivalent noise bandwidth in Hz (which is the same power that would result from a brickwall filter of that bandwidth if white noise was at the input) and $T$ is the duration of the data in seconds.

This relationship applies universally to continuous and discrete systems given the relationship between a rectangular window and its Fourier Transform in both the CTFT and DFT (A Sinc for the CTFT and essentially an aliased Sinc function in the DFT, approaching the Sinc as the number of samples increases).


NRZ vs. Manchester pulse shaping filter (Modulator)

Manchester coding effectively creates a bit sequence that is twice as long as the given bit sequence and then applies NRZ pulses to the resulting bit stream.
There are two types of Manchester coding. In one form, each $0$ in the original
sequence is replaced by $01$ and each $1$ by $10$ while the other form complements
these bit patterns.  Using the first form,
$$101011 \to 10\,01\,10\,01\,10\,10$$ while the second form gives the complementary
bit sequence $01\,10\,01\,10\,0101$. If each data bit is of duration $T$, the NRZ
waveform with Half Sine pulses has sines of period $2T$ so that each bit
pulse has duration half a period. Thus, that leading $1$ gets
modulated to $\sin(\pi t/T)$ lasting from $t = 0$ till $t = T$.
With Manchester coding, each channel bit
(of duration $T/2$) gets a Half Sine of a sine of period $T$. Put another
way,

with plain NRZ, that first $1$ in your bit sequence will have
pulse $\sin(\pi t/T)$ lasting from $t = 0$ to $t = T$. This is a positive
going pulse rising from $0$ at $t=0$, peaking at $1$ at $t=T/2$, and
falling to $0$ at $t = T$
with Manchester coding, that data $1$ will become channel $10$ giving a positive pulse $\sin(2\pi t/T)$ from $0$ to $T/2$. The next channel bit
$0$ has the delayed half-sine pulse $\sin(2\pi (t-T/2)/T$ from $t=T/2$
to $t = T$, but with a negative amplitude (the channel bit
during this interval is $0$ instead of a $1$) which
works out to be
$$\begin{align*}
-\sin(2\pi(t-T/2)/T) &= -\sin(2\pi t/T - \pi)\\
&= -\sin(2\pi t/T)\cos(\pi) + \cos(2\pi t/T)\sin(\pi)\\
&= \sin(2\pi t/T)
\end{align*}$$
lasting from $t = T/2$ to $t = T$. In short, during the entire data
bit interval $[0,T]$ the signal can be expressed as $\sin(2\pi t/T)$
and is one cycle of the sine wave of period T. If the data
bit were a $0$ instead of a $1$, the channel bits would be $01$
and we would get $-\sin(2\pi t/T)$ from $t=0$ to $t=T$. In short,
Manchester coding with half-sine pulses
effectively creates a phase-modulated signal
from the data sequence with a carrier frequency of $1/T$ Hz and
each data bit interval contains exactly one period of the sinusoid.
If you look at the figure for Thomas's Manchester coding
in Wikipedia,
what we have is the $\pm$ rectangular pulses comprising each
data bit interval, which look like a hard-limited sinusoid
of period $T$, being replaced by the sinusoid.

  The OP says, "My thoughts are that with a Half Sine ....  pulse shaping
  filter, I will see the transitions between the 1's and 0's smooth out,..."
  but this does not happen, there will be abrupt reversals of carrier
  phase at those data bit boundaries where there is a transition
  from $0$ to $1$ or from $1$ to $0$, pretty much as happens
  in plain PSK.


What is the inverse Laplace transform of squared denominator term?

This is quite straightforward to solve. Either just use this table where you can directly find the corresponding result, or "derive" it yourself with very basic knowledge of the Laplace transform.

You should know that

$$\mathcal{L}\{u(t)\}=\frac{1}{s}\tag{1}$$

Integration is equivalent to multiplication with $1/s$, so integrating $(1)$ gives

$$\mathcal{L}\{t\cdot u(t)\}=\frac{1}{s^2}\tag{2}$$

Replacing $s$ by $s+a$ is equivalent to multiplication with $e^{-at}$. Consequently, we have

$$\mathcal{L}^{-1}\left\{\frac{2}{(s+2)^2}\right\}=2\,t\,e^{-2t} u(t)\tag{3}$$


High audio frequencies in television broadcast

In digital television, I think they use mostly a 48 kHz sampling frequency, so it is up to whether the audio encoder allocates any bits to the ultrasonic frequencies:

Figure 1. Before encoding: Left) almost full scale white noise, Center) almost full scale pink noise, Right) full scale 20-24 kHz band pass white noise.

Figure 2. After encoding in MPEG 1.0 layer II, 224 kbit/s, 48000 Hz stereo (known to be used in DVB-T) and decoding, using mp2enc for encoding and mpg123 for decoding.

With white noise, in the 16-20 kHz range only some bursts of frequency survive, but above that nothing. Starting with "pink" noise, everything above 15 kHz gets stripped away, because more bits are allocated to the higher-intensity lower frequencies. Starting with content only >20 kHz, some ultrasonic frequencies survive, as no bits need to be allocated to lower frequencies.

Also other encoding methods than MPEG 1.0 layer II are in use in broadcast television.


What you might want to investigate is FFT fast convolution using overlap-add or overlap-save methods. You will need to expand the length of each FFT by the length of the impulse of your desired filter. This is because (1) FFT/IFFT convolution is circular, and (2) each index in the FFT array result corresponds to almost all frequencies (a Sinc shaped response), not just one (even if mostly concentrated near one), so any single bin modification will leak throughout the entire frequency response (except certain exact periodic frequencies).


I think you misinterpret the term inverse here.

Given a window sequence $w[n]$ of length $L$ with a region of support such that $ w[n] \not = 0 ,~~$ for $ 0 \leq n \leq L-1 $ , and a signal $x[n]$ you get the windowed signal $v[n] = w[n]x[n]$ , as each sample of the signal $x[n]$ is multipied by the corresponding sample of the window $w[n]$.

Inverting $v[n]$ means, given the windowed sequence $v[n]$ you obtain the original (unweighted) segment of $x[n]$ that resides in the region of suport range of the window $w[n]$ and that's given by: $$x[n] = \frac{v[n]}{w[n]} ~~,~~ w[n] \not = 0$$

The operation of recovering the original signal segment is as follows:
Let the original signal be the Matlab vector:

x = [1 2 3 0 5];   % original signal segment with 5 samples (L=5)

let the window be a Hamming type length of 5 samples:

w = hamming(5)';   % 5 sample Hamming window w = [0.08 0.54 1.0 0.54 0.08]

Hence the windowed signal is:

v = x.*w;          % v = [0.08 1.08 3.0 0.0 0.4] (windowed signal)

Now to recover x[n] from v[n] perform the following:

xrec = v./w;      % xrec = [1.0 2.0 3.0 0.0 5.0]; recovered

The operation $x_{rec}[n] = v[n]/w[n]$ is implemented per sample basis, such that $x_{rec}[0] = v[0]/w[0]$ or $x_{rec}[1] = v[1]/w[1]$ or in general
$$x_{rec}[n] = v[n]/w[n]$$

Now however the Hann window has the property that its initial and final samples are zero (in effect reducing the window length by 2) therefore an overlapping implementation is necessary to be able to recover $x[n]$ back from multiple $v[n]$ segments, as follows:

Let the original signal be the Matlab vector:

x = [1 2 3 4 5 4 3 2 1];   % original signal segment with 9 samples (L=9)

let's use a Hann window of length 5:

w = hann(5)';   % 5 sample Hann window w = [0.0 0.5 1.0 0.5 0.0]

First embed a zero to the first and last samples of x[n]:

xe = [0 x 0];   % xe = [0 1 2 3 4 5 4 3 2 1 0] (xe has 11 samples)

Now apply windowing into overlapping segments as such:

v1 = x[1:5].*w;          % v1 = [0.0 0.5 2.0 1.5 0.0] (first segment)
v2 = x[4:8].*w;          % v2 = [0.0 2.0 5.0 2.0 0.0] (second segment)
v3 = x[7:11].*w;         % v3 = [0.0 1.5 2.0 0.5 0.0] (third segment)

(notice the MATLAB indice operation x[i:j] to select portions of the input vector)

Now to recover $x[n]$ from $v_k[n]$ perform the following:
First note that each segment begins and ends in zero samples that you should discard, and you should consider only the remainig $L-2$ vetween those.

Then consider also the amount of overlap so that you don't recompute a previously computed input sample, then:

x1_rec = v1[2:4]./w[2:4];      % x1_rec = [1.0 2.0 3.0]; recovered first 3 samples of x[n]

x2_rec = v2[2:4]./w[2:4];      % x2_rec = [4.0 5.0 4.0]; recovered middle 3 samples of x[n]

x3_rec = v3[2:4]./w[2:4];      % x3_rec = [3.0 2.0 1.0]; recovered last 3 samples of x[n]


It is really quite simple for values between 1 and -1:

valueDBFS = 20*log10(abs(value))

If you have values between another range, for example 16bit, then it's:

valueDBFS = 20*log10(abs(value)/32768)

(because 16bit signed has values between -32768 and +32767, e.g. 2^15=32768)

And so the general formula is:

valueDBFS = 20*log10(abs(value)/maxValue)

The code examples are given in MATLAB code.


multiplication in frequency domain equals convolution in time domain mismatch

You missed the part where under the DFT, multiplication in one domain is equivalent to circular convolution in the other. Standard problem!

Just pad the signals you want to convolve with enough zeros on both sides, and things will look better.

Also, I hope you're just doing this as to exercise your matlab coding, because getting the phase out of two complex numbers just to multiply them is way, way, way more computationally intense than just directly multiplying the two numbers.


How does noise reduction for speech recognition differ from noise reduction that is supposed to make speech more "intelligible" for humans?

I don't really find papers that discuss this difference.

There are whole books on the subject:

Robust Automatic Speech Recognition 1st Edition

  Do speech intelligibility and speech quality correlate with the accuracy of speech recognition systems?

Usually no, usually noise reduction corrupts features in unpredictable way and reduces speech recognition accuracy.

  Are there objective measures that can evaluate how "good" a denoised speech signal will be for a speech recognition system, for example if also given the original clean speech? Or is the only way to find out how good your noise reduction technique is, to train the speech recognition system on the denoised data and look at the accuracy?

Second. Moreover feature-based noise reduction actually removes important information from the spectrum altogether so you can not repair an accuracy of the clean system. For that reason modern approach is to perform multi-style training on noisy data instead of using noise reduction algorithm beforehand. It ends in more accurate recognition.


en.wikipedia.org/wiki/Wiener_filter

what you are looking for is a Weiner  filter  that takes into account the standard deviations of your data. See http://www.astroml.org/book_figures/chapter10/fig_wiener_filter.html
 and the description in Wikipedia  for implementation.  Basically you will do what you were suggesting,  but with an additional term in the denominator of the ratio of fft's  that is the fft  of the noise.


Definition of sampling using delta or indicator function?

that equation following "Conceptual scheme" is misleading.  the $\implies$ arrow is not an equal sign.  this is what is true:

$$\begin{align}
x(t) \times \sum\limits_{n=-\infty}^{\infty}\delta(t-nT_s) &= \sum\limits_{n=-\infty}^{\infty}x(t)\times\delta(t-nT_s) \\
&= \sum\limits_{n=-\infty}^{\infty}x(nT_s)\times\delta(t-nT_s) \\
&= \sum\limits_{n=-\infty}^{\infty}x[n]\times\delta(t-nT_s) \\
\end{align}$$

because we simply define $x[n] \triangleq x(nT_s)$.

the difference between the dirac delta and the indicator function is that one integrates to an area of 1 and the other integrates to an area of 0.


Derivative of mutual information in EMMA algorithm for image alignment

I am not an expert on mutual information but I understand a fair bit of the mathematics here.
So, I read the paper and indeed there seems some typo in the equation $(17)$. I am yet to figure that out.
But I have a way for you to compute the derivate as per equation $(17)$. So, if we go back to the equation,

then it seems that equation $17$ is essentially subtraction of the two terms above. Now, the differentiation of $\frac{d}{dT}h(v(T(x)))$ and $\frac{d}{dT}h(u(x),v(T(x)))$ can be computed using equation $16$.

So, you could develop a function (in Python or any language of your choice) as per equation $16$ and just use that function twice, once for $h(v(T(x)))$ and $h(u(x),v(T(x)))$. And then, subtract the two values. This way you wouldn't have the dimensionality issue.
Hope this was sufficient. On a side note, this is essentially an optimization problem where your cost function to optimize is the mutual information $I$. In my work on optimization problems we often do this function subtraction trick.
P.S. You can even reverse engineer and try different variations for equation $17$ and then, figure out your issue of dimensionality.