Solving Scale and Permutation indeterminacies in Blind Source Separation

What i know in implementations of BSS (e.g jade) is that they first normalise the input data (to variance = 1) and then order the data according to scale. So this is unique and will produce consistent output ordering for the same input data (will add references shortly).

Have in mind however that this is just a consistent heuristic (any consistent heuristic will do) as BSS methods by definition are clue-less about scaling and permutation of the inputs.

References of permutation problem for jade:

The following heuristic is described in the comments section of the matlab version of jade:

% Usage:
%   * If X is an nxT data matrix (n sensors, T samples) then
%     B=jadeR(X) is a nxn separating matrix such that S=B*X is an nxT
%     matrix of estimated source signals.
%   * If B=jadeR(X,m), then B has size mxn so that only m sources are
%     extracted.  This is done by restricting the operation of jadeR
%     to the m first principal components.
%   * Also, the rows of B are ordered such that the columns of pinv(B)
%     are in order of decreasing norm; this has the effect that the
%     `most energetically significant' components appear first in the
%     rows of S=B*X.


but... "shaping the spectrum" is exactly the definition of a filter. I think you can justifiably call everything that shapes a spectrum a filter, so no, nothing but filtering shapes a spectrum.
Note that just because that's the most commonly analyzed case, filters are more than just linear filters.


How to detect corners in a binary images with OpenGL?

I just happened to be implementing something like this on OpenGL ES 2.0 using Harris corner detection, and while I'm not completely finished, I thought I'd share the shader-based implementation I have so far. I've done this as part of an iOS-based open source framework, so you can check out the code if you're curious as to how some particular step works.

To do this, I use the following steps:

Reduce the image to its luminance values using a dot product of the RGB values with the vector (0.2125, 0.7154, 0.0721).
Calculate the X and Y derivatives by subtracting the red channel values from the pixels left and right and above and below the current pixel. I then store the x derivative squared in the red channel, the Y derivative squared in the green channel, and the product of X and Y derivatives in the blue channel. The fragment shader for this looks like the following:

precision highp float;

varying vec2 textureCoordinate;
varying vec2 leftTextureCoordinate;
varying vec2 rightTextureCoordinate;

varying vec2 topTextureCoordinate;
varying vec2 bottomTextureCoordinate;

uniform sampler2D inputImageTexture;

void main()
{
 float topIntensity = texture2D(inputImageTexture, topTextureCoordinate).r;
 float bottomIntensity = texture2D(inputImageTexture, bottomTextureCoordinate).r;
 float leftIntensity = texture2D(inputImageTexture, leftTextureCoordinate).r;
 float rightIntensity = texture2D(inputImageTexture, rightTextureCoordinate).r;

 float verticalDerivative = abs(-topIntensity + bottomIntensity);
 float horizontalDerivative = abs(-leftIntensity + rightIntensity);

 gl_FragColor = vec4(horizontalDerivative * horizontalDerivative, verticalDerivative * verticalDerivative, verticalDerivative * horizontalDerivative, 1.0);
}

where the varyings are just the offset texture coordinates in each direction. I precalculate these in the vertex shader to eliminate dependent texture reads, which are notoriously slow on these mobile GPUs.
Apply a Gaussian blur to this derivative image. I used a separated horizontal and vertical blur, and take advantage of hardware texture filtering to do a nine-hit blur with only five texture reads on each pass. I describe this shader in this Stack Overflow answer.
Run the actual Harris corner detection calculation using the blurred input derivative values. In this case, I'm actually using the calculation described by Alison Noble in her Ph.D. dissertation "Descriptions of Image Surfaces". The shader that handles this looks like the following:

varying highp vec2 textureCoordinate;

uniform sampler2D inputImageTexture;

const mediump float harrisConstant = 0.04;

void main()
{
 mediump vec3 derivativeElements = texture2D(inputImageTexture, textureCoordinate).rgb;

 mediump float derivativeSum = derivativeElements.x + derivativeElements.y;

 // This is the Noble variant on the Harris detector, from
 // Alison Noble, "Descriptions of Image Surfaces", PhD thesis, Department of Engineering Science, Oxford University 1989, p45.
 mediump float harrisIntensity = (derivativeElements.x * derivativeElements.y - (derivativeElements.z * derivativeElements.z)) / (derivativeSum);

 // Original Harris detector
 //     highp float harrisIntensity = derivativeElements.x * derivativeElements.y - (derivativeElements.z * derivativeElements.z) - harrisConstant * derivativeSum * derivativeSum;

 gl_FragColor = vec4(vec3(harrisIntensity * 10.0), 1.0);
}

Perform local non-maximum suppression and apply a threshold to highlight the pixels that pass. I use the following fragment shader to sample the eight pixels in the neighborhood of a central pixel and identify whether or not it is the maximum in that grouping:

uniform sampler2D inputImageTexture;

varying highp vec2 textureCoordinate;
varying highp vec2 leftTextureCoordinate;
varying highp vec2 rightTextureCoordinate;

varying highp vec2 topTextureCoordinate;
varying highp vec2 topLeftTextureCoordinate;
varying highp vec2 topRightTextureCoordinate;

varying highp vec2 bottomTextureCoordinate;
varying highp vec2 bottomLeftTextureCoordinate;
varying highp vec2 bottomRightTextureCoordinate;

void main()
{
    lowp float bottomColor = texture2D(inputImageTexture, bottomTextureCoordinate).r;
    lowp float bottomLeftColor = texture2D(inputImageTexture, bottomLeftTextureCoordinate).r;
    lowp float bottomRightColor = texture2D(inputImageTexture, bottomRightTextureCoordinate).r;
    lowp vec4 centerColor = texture2D(inputImageTexture, textureCoordinate);
    lowp float leftColor = texture2D(inputImageTexture, leftTextureCoordinate).r;
    lowp float rightColor = texture2D(inputImageTexture, rightTextureCoordinate).r;
    lowp float topColor = texture2D(inputImageTexture, topTextureCoordinate).r;
    lowp float topRightColor = texture2D(inputImageTexture, topRightTextureCoordinate).r;
    lowp float topLeftColor = texture2D(inputImageTexture, topLeftTextureCoordinate).r;

    // Use a tiebreaker for pixels to the left and immediately above this one
    lowp float multiplier = 1.0 - step(centerColor.r, topColor);
    multiplier = multiplier * 1.0 - step(centerColor.r, topLeftColor);
    multiplier = multiplier * 1.0 - step(centerColor.r, leftColor);
    multiplier = multiplier * 1.0 - step(centerColor.r, bottomLeftColor);

    lowp float maxValue = max(centerColor.r, bottomColor);
    maxValue = max(maxValue, bottomRightColor);
    maxValue = max(maxValue, rightColor);
    maxValue = max(maxValue, topRightColor);

    gl_FragColor = vec4((centerColor.rgb * step(maxValue, centerColor.r) * multiplier), 1.0);
}

This process generates a cornerness map from your objects that looks like this:

The following points are identified as corners based on the non-maximum suppression and thresholding:

With proper thresholds set for this filter, it can identify all 16 corners in this image, although it does tend to place the corners a pixel or so inside the actual edges of the object.

On an iPhone 4, this corner detection can be run at 20 FPS on 640x480 frames of video coming from the camera, and an iPhone 4S can easily process video of that size at 60+ FPS. This should be a good deal faster than CPU-bound processing for a task like this, although right now the process of reading back the points is CPU-bound and a little slower than it should be.

If you want to see this in action you can grab the code for my framework and run the FilterShowcase example that comes with it. The Harris corner detection example there runs on live video from the device camera, although as I mentioned the reading back of corner points currently occurs on the CPU, which is really slowing this down. I'm moving to a GPU-based process for this, as well.


Hints

Define the frequency of the car as $f_0$, the heard frequency when it's approaching as $f1$ and the heard frequency when it's departing as $f_2$
Write the Doppler formula for both $f_1$ and $f_2$ using a receiver velocity of 0 (since the policeman is stationary with respect to the air) and a car velocity of $v_s$. Pay attention to the sign of the source velocity.
Write the ratio of both formulas and equal it to a major third. $f_0$ should drop out and you are left with a single equation for $v_s$
Solve for $v_s$


First of all, what is the order of your IIR filter? The highest order I have ever used was an order-10 IIR filter for a control loop application. I feel like it is unlikely that you need more that this.

Second, it is a good idea to split your filter in second-order-sections (SOS) and cascade them , this usually fix most issues.

https://www.dsprelated.com/showarticle/1137.php

You can also split your filter in second-order-sections and put them in parallel and add their ouputs too.


Here are a few things that aren't right:

1) Regardless of the Carrier Synchronization:
 You are using 0 and 90 degrees for your signaling. You would need an I and Q channel in order to decode the signal if that were the case. For BPSK (Binary) you should be using 0 and 180 degrees. Only if you are using QPSK should you use 90 and 270 degrees, then again its just mimicking using an I and Q channels.

There is a system block diagram below that will make this more apparent.

2) DPSK does not require any phase carrier synchronization. You pay for it by halving your SNR.

Encoding DQPSK

(for BPSK, eliminate Q and sqrt(2)

I-inphase channel Bit values
Q-quadrature bit values
C = (I + J*Q)/sqrt(2); // recombined I and Q channels
for n = 0 to length( C )
    x_diff = x_diff(n-1) *C(n);  //where X(n) <0 == 0
//X_diff is the differentially encoded channels

Decoding DQPSK

I - retrieved I channel
Q - retrieved Q channel
C =I + j*Q;
for n = 0 to length( C )
    decoded(n) =C(n)/decoded(n-1);  //where X(n) <0 == 0

4) You should never see a phase value because of thresholding at the sampler block

Transmitter Block Diagram

Receiver Block Diagram

(Remember you can ignore anything with the Q channel path if you are using BPSK, ignore the Convolutional encoding block and the input string block this whole thing will be easier to understand) .

The Sampler block acts like a Analog to Digital Converter. It samples the output of the I and Q channels at the bitrate your transmitting at (100bits/s, 200bits/s etc) it will give you the AMPLITUDE of the signal recovered.

Because, hopefully, you've used NRZ encoding, (-1, +1) to signal 0 and 1, you must threshold your signal. All values below Zero are 0's and all values above zero are 1's.  If you've used RZ encoding that is (0, 1) instead of (-1, +1) for your signaling you should threshold at 0.5.

Therefore you should NEVER work with a phase value.

5) Additional Block Diagram Information:

This is purely to be thorough and give readers a better idea about what is going on.

 1) The Users inputs a string ("Hello World" for example) this gets converted into its
    Ascii binary value.
3)This then gets put through a convolutional encoder, this is for feedfoward error
  correction. This is to reduce the transmission errors due to noise in the
  communication channel. Like power-lines or someone talking.
4)The binary values are converted into Non-return to Zero Encoding (NRZ) this means
  that 0's are converted to -1's and 1's are converted to +1's.
5) differential encoding is applied to eliminate the need for carrier (discussed above)
    synchronization.
6)The Signal is split into the I and Q channels (all even bits go to Q. All odd bits
   go to the I channel)
7) They are then shaped into pulses by the wave generator.
8) An Inter-Symbol Interference filter is applied (root raised cosine filter).
9)The I and Q channels are multipled by cosine and sine (same frequency, just shifted
    90 degrees) then recombined and transmitted.
10)The Receiver process is the same thing as the transmitter in reverse except:
   Convolutional Decoder: Use the Viterbi algorithm.
   Detector Block - this is where the NRZ thresholding takes place, values greater
   than Zero  are 1's values less than Zero are 0's

6) Working with M-ary (like Bin-ary, or trin-ary)

If you want to work with more data points in QPSK, you need to generate your constellation diagram.

Remember that all points on a QPSK system should lie on a circle.

Think of each point as a combination of the I and Q channels and not as purely a Phase value used for signaling.


Why is my MATLAB's bode plot wildly off?

I suspect that the reason for the difference is the number of samples used in the plot, since the extreme points are likely going to infinity (or very large numbers) such that if less points are used, those large values are simply missed in the computation. If the OP uses the same number of points in each, then the same result should be achieved- and for all the points shown if they were super-imposed on each plot it would be clear that for the points given, the results are the same.
My first suspicion before the OP added the plots was that it was due to the units of the axis. That is not the case but leaving this part of my answer below as it is another common reason for Bode plot differences:
The units of the frequency axis should be confirmed as they could be in units of:
Continuous domain units of frequency:
Hertz, or cycles/sec
radians/sec
Discrete domain units of (normalized) frequency:
cycles per sample (unique domain is $[-1,1)$)
radians/sample (unique domain is $[-\pi, \pi)$)


Proper generation of Complex white Gaussian noise using Matlab

Refer to the code below that generates some noise at a given SNR:

N = 100000;

% Generate some random signal
signal = randn(N, 1) + 1j*randn(N,1);

% Here, the signal power is 2 (1 for real and imaginary component)
signalPower_lin = 1/N*signal'*signal

% This corresponds to 3dB (assuming power=1 is 0dB)
signalPower_dB = 10*log10(signalPower_lin)

noisePower_dB = signalPower_dB - 30  % do 30dB SNR
noisePower_lin = 10^(noisePower_dB/10)

noise = sqrt(noisePower_lin/2) * (randn(N,1) + 1j*randn(N,1));

received = signal + noise;

The program outputs:

signalPower_lin =

   2.0042 + 0.0000i

signalPower_dB =

   3.0193 + 0.0000i

noisePower_dB =

 -26.9807 + 0.0000i

noisePower_lin =

   0.0020 + 0.0000i

Note that generating a complex noise of variance 1, you need to do

noise = sqrt(1/2) * (randn(N,1) + 1j*randn(N,1))

Since each component (real and imaginary) needs to have variance 1/2, such that their sum becomes 1.

To answer your points:
1) As a rule of thumb when to use 20 and when to use 10:
If you describe Powers or Energies, the factor is 10. If you describe amplitudes, the factor is 20. Then, in dB-scale both would have the same value, since:

$$
\begin{align}
P&=A^2 \text{ (Power = Amplitude squared)}\\
10\log(P) &= 10\log(A^2)\\
&=20\log(A)
\end{align}
$$

So, anytime you use a values which creates something meaningful when squared (like amplitude), use 20. Everytime, you have something that is meaningful when taken the square root from (like Power), use 10.

Note that in my code, I use noise = sqrt(noise_power/2) * randn(...) to generate the noise, i.e. I transform the noise power into the noise amplitude prior to multiplying the normal random variable (which corresponds to changing the factor 10 to 20 in the exponent). Also remember the rule
$$Var(aX) = a^2Var(X)$$
which is exactly the reason why you need to take the square root of the power, such that your resulting noise has the required power.

2) Here, I just transform noisePower_db into noisePower_linear, which is by definition linear=10^(db/10). Note that $10^{-SNR/10}$ is similar, since SNR=Signal/Noise, i.e. the noise is in the denominator (consider $10^{-1/x}=10^x$.


How to plot poles, zeros and impulse response of a transfer function

Two hints:

The standard way of writing a transfer function in the $\mathcal Z$ domain is a rational function in $z^{-1}$, not in $z$. Try this
Use tf2zp() instead of residuez()


I'll use the shorthand window for "window function".

With audio, any processing that creates something akin to pre-ringing or pre-echo will sound sloshy like a low bit-rate mp3. This happens when localized energy of a transient or an impulse is spread backwards in time, for example by modification of the spectral data in lapped transforms such as lapped modified discrete cosine transform (MDCT). In such processing, audio is windowed by overlapping analysis windows, transformed, processed in frequency domain (like data-compressed to a smaller bitrate), windowed again with a synthesis window and summed back together. The product of the analysis and synthesis window must be such that the overlapping windows sum to unity.

Traditionally the window functions used have been symmetric, and their width has been a compromise between frequency selectivity (long window) and time-domain artifact avoidance (short window). The wider the window, the more back in time the processing can spread the signal. A more recent solution is to use an asymmetrical window. The two windows used can be mirror images of each other. The analysis window drops from peak to zero fast so that impulses are not "detected" much in advance, and the syntheis window rises from zero to peak fast, so that the effects of any processing are not spread much backwards in time. Another advantage of this is low latency. The asymmetrical windows can have good frequency-selectivity, and can replace variable-sized symmetrical windows in audio compression, like a kind of cure-all. See M. Schnell, M. Schmidt, M. Jander, T. Albert,
R. Geiger, V. Ruoppila, P. Ekstrand,
M. Lutzky, B. Grill, “MPEG-4 Enhanced Low
Delay AAC - a new standard for high quality
communication”, 125th AES Convention, San
Francisco, CA, USA, preprint 7503, Oct. 2008 and another conference paper where they show also the magnitude of the Fourier transform of their window: Schnell, M., et al. 2007. Enhanced MPEG-4 Low Delay AAC – Low Bitrate High Quality Communication. In 122th AES Convention.

Figure 1. Illustration of use of asymmetrical windows in lapped analysis-processing-synthesis. The product (black dashed) of the analysis-window (blue) and the synthesis window (yellowish orange) sums to unity with the window from the previous frame (grey dashed). Further constraints are needed to guarantee perfect reconstruction when using MDCT.

Discrete Fourier transform (DFT, FFT) could be used instead of MDCT, but will in such contexts give redundant spectral data. Compared to DFT, MDCT gives only half of the spectral data while still enabling perfect reconstruction if suitable windows are chosen.

Here is my own asymmetrical window design (Fig. 2) suitable for lapped analysis-processing-synthesis using DFT but not MDCT with which it does not give perfect reconstruction. The window tries to minimize the product of mean-square time and frequency bandwidths (similarly to the confined Gaussian window) while retaining some potentially useful time-domain properties: nonnegative, unimodal with the peak at "time zero" around which the analysis and synthesis windows are mirror images of each other, function and first derivative continuity, zero-mean when the square of the window function is interpreted as an unnormalized probability density function. The window was optimized using differential evolution.

Figure 2. Left: An asymmetrical analysis window suitable for overlapped analysis-processing-resynthesis together with its time-reversed counterpart synthesis window. Right: Cosine window, with the same latency as the asymmetrical window

Figure 3. Magnitude of the Fourier transforms of the cosine window (blue) and the asymmetrical window (orange) of Fig. 2. The asymmetrical window shows better frequency selectivity.

Here is Octave source code for the plots and for the asymmetrical window. The plotting code comes from Wikimedia Commons. On Linux I recommend installing gnuplot, epstool, pstoedit, transfig first and librsvg2-bin for viewing using display.

pkg load signal

graphics_toolkit gnuplot
set (0, "defaultaxesfontname", "sans-serif")
set (0, "defaultaxesfontsize", 12)
set (0, "defaultaxeslinewidth", 1)

function plotWindow (w, wname, wfilename = "", wspecifier = "", wfilespecifier = "")

  M = 32; % Fourier transform size as multiple of window length
  Q = 512; % Number of samples in time domain plot
  P = 40; % Maximum bin index drawn
  dr = 130; % Maximum attenuation (dB) drawn in frequency domain plot

  N = length(w);
  B = N*sum(w.^2)/sum(w)^2 % noise bandwidth (bins)

  k = [0 : 1/Q : 1];
  w2 = interp1 ([0 : 1/(N-1) : 1], w, k);

  if (M/N < Q)
    Q = M/N;
  endif

  figure('position', [1 1 1200 600])
  subplot(1,2,1)
  area(k,w2,'FaceColor', [0 0.4 0.6], 'edgecolor', [0 0 0], 'linewidth', 1)
  if (min(w) >= -0.01)
    ylim([0 1.05])
    set(gca,'YTick', [0 : 0.1 : 1])
  else
    ylim([-1 5])
    set(gca,'YTick', [-1 : 1 : 5])
  endif
  ylabel('amplitude')
  set(gca,'XTick', [0 : 1/8 : 1])
  set(gca,'XTickLabel',[' 0'; ' '; ' '; ' '; ' '; ' '; ' '; ' '; 'N-1'])
  grid('on')
  set(gca,'gridlinestyle','-')
  xlabel('samples')
  if (strcmp (wspecifier, ""))
    title(cstrcat(wname,' window'), 'interpreter', 'none')
  else
    title(cstrcat(wname,' window (', wspecifier, ')'), 'interpreter', 'none')
  endif
  set(gca,'Position',[0.094 0.17 0.38 0.71])

  H = abs(fft([w zeros(1,(M-1)*N)]));
  H = fftshift(H);
  H = H/max(H);
  H = 20*log10(H);
  H = max(-dr,H);
  k = ([1:M*N]-1-M*N/2)/M;
  k2 = [-P : 1/M : P];
  H2 = interp1 (k, H, k2);

  subplot(1,2,2)
  set(gca,'FontSize',28)
  h = stem(k2,H2,'-');
  set(h,'BaseValue',-dr)
  xlim([-P P])
  ylim([-dr 6])
  set(gca,'YTick', [0 : -10 : -dr])
  set(findobj('Type','line'),'Marker','none','Color',[0.8710 0.49 0])
  grid('on')
  set(findobj('Type','gridline'),'Color',[.871 .49 0])
  set(gca,'gridlinestyle','-')
  ylabel('decibels')
  xlabel('bins')
  title('Fourier transform')
  set(gca,'Position',[0.595 0.17 0.385 0.71])

  if (strcmp (wfilename, ""))
    wfilename = wname;
  endif
  if (strcmp (wfilespecifier, ""))
    wfilespecifier = wspecifier;
  endif
  if (strcmp (wfilespecifier, ""))
    savetoname = cstrcat('Window function and frequency response - ', wfilename, '.svg');
  else
    savetoname = cstrcat('Window function and frequency response - ', wfilename, ' (', wfilespecifier, ').svg');
  endif
  print(savetoname, '-dsvg', '-S1200,600')
  close

endfunction

N=2^17; % Window length, B is equal for Triangular and Bartlett from 2^17
k=0:N-1;

w = -cos(2*pi*k/(N-1));
w .*= w > 0;
plotWindow(w, "Cosine")

freqData = [0.66697133904805994131, -0.20556692772918355727, 0.49267389481655493588, -0.25062332863369246594, -0.42388422228212319087, 0.42317609537724842905, -0.03930334287740060856, -0.11936153294075849129, 0.30201210285940127687, -0.15541616804857899536, -0.16208119255594669039, 0.12843871362286504723, -0.04470810646117385351, -0.00521885027256757845, 0.07185811583185619522, -0.02835116723496184862, -0.01393644785822748498, 0.00780746224568363342, -0.00748496824751256583, 0.00119325723511989282, 0.00194602547595042175];
freqData(1) /= 2;
scale = freqData(1) + sum(freqData.*not(mod(1:length(freqData), 2)));
freqData /= scale;
w = freqData(1)*ones(1, N);
for bin = 1:(length(freqData)/2)
  w += freqData(bin*2)*cos(2*pi*bin*((1:N)-1)/N);
  w += freqData(bin*2+1)*sin(2*pi*bin*((1:N)-1)/N);
endfor
w(N/4+1:N/2+1) = 0;
w(N/8+2:N/4) = (1 - w(N/8:-1:2).*w(7*N/8+2:N))./w(7*N/8:-1:6*N/8+2);
w = shift(w, -N/2);
plotWindow(w, "Asymmetrical");

You may want to use only every second sample of the window because it starts and ends at zero. The following C++ code does that for you so you don't get any zero samples except for in one quarter of the window that is zero everywhere. For the analysis window this is the first quarter and for the synthesis window this is the last quarter. The second half of the analysis window should be aligned with the first half of the synthesis window for calculation of their product. The code also tests the mean of the window (as a probability density function), and showcases the flatness of overlapped reconstruction.

#include <stdio.h>
#include <math.h>

int main() {
  const int windowSize = 400;
  double *analysisWindow = new double[windowSize];
  double *synthesisWindow = new double[windowSize];
  for (int k = 0; k < windowSize/4; k++) {
    analysisWindow[k] = 0;
  }
  for (int k = windowSize/4; k < windowSize*7/8; k++) {
    double x = 2 * M_PI * ((k+0.5)/windowSize - 1.75);
    analysisWindow[k] = 2.57392230162633461887-1.58661480271141974718*cos(x)+3.80257516644523141380*sin(x)
      -1.93437090055110760822*cos(2*x)-3.27163999159752183488*sin(2*x)+3.26617449847621266201*cos(3*x)
      -0.30335261753524439543*sin(3*x)-0.92126091064427817479*cos(4*x)+2.33100177294084742741*sin(4*x)
      -1.19953922321306438725*cos(5*x)-1.25098147932225423062*sin(5*x)+0.99132076607048635886*cos(6*x)
      -0.34506787787355830410*sin(6*x)-0.04028033685700077582*cos(7*x)+0.55461815542612269425*sin(7*x)
      -0.21882110175036428856*cos(8*x)-0.10756484378756643594*sin(8*x)+0.06025986430527170007*cos(9*x)
      -0.05777077835678736534*sin(9*x)+0.00920984524892982936*cos(10*x)+0.01501989089735343216*sin(10*x);
  }
  for (int k = 0; k < windowSize/8; k++) {
    analysisWindow[windowSize-1-k] = (1 - analysisWindow[windowSize*3/4-1-k]*analysisWindow[windowSize*3/4+k])/analysisWindow[windowSize/2+k];
  }
  printf("Analysis window:\n");
  for (int k = 0; k < windowSize; k++) {
    printf("%d\t%.10f\n", k, analysisWindow[k]);
  }
  double accu, accu2;
  for (int k = 0; k < windowSize; k++) {
    accu += k*analysisWindow[k]*analysisWindow[k];
    accu2 += analysisWindow[k]*analysisWindow[k];
  }
  for (int k = 0; k < windowSize; k++) {
    synthesisWindow[k] = analysisWindow[windowSize-1-k];
  }
  printf("\nSynthesis window:\n");
  for (int k = 0; k < windowSize; k++) {
    printf("%d\t%.10f\n", k, synthesisWindow[k]);
  }
  printf("Mean of square of analysis window as probability density function:\n%f", accu/accu2);
  printf("\nProduct of analysis and synthesis windows:\n");
  for (int k = 0; k < windowSize/2; k++) {
    printf("%d\t%.10f\n", k, analysisWindow[windowSize/2+k]*synthesisWindow[k]);
  }
  printf("\nSum of overlapping products of windows:\n");
  for (int k = 0; k < windowSize/4; k++) {
    printf("%d\t%.10f\n", k, analysisWindow[windowSize/2+k]*synthesisWindow[k]+analysisWindow[windowSize/2+k+windowSize/4]*synthesisWindow[k+windowSize/4]);
  }
  delete[] analysisWindow;
  delete[] synthesisWindow;
}

And the source code for the optimization cost function to be used with Kiss FFT and an optimization library:

class WinProblem : public Opti::Problem {
private:
  int numParams;
  double *min;
  double *max;
  kiss_fft_scalar *timeData;
  kiss_fft_cpx *freqData;
  int smallSize;
  int bigSize;
  kiss_fftr_cfg smallFFTR;
  kiss_fftr_cfg smallIFFTR;
  kiss_fftr_cfg bigFFTR;
  kiss_fftr_cfg bigIFFTR;

public:
  // numParams must be odd
  WinProblem(int numParams, int smallSize, int bigSize, double* candidate = NULL) : numParams(numParams), smallSize(smallSize), bigSize(bigSize) {
    min = new double[numParams];
    max = new double[numParams];
    if (candidate != NULL) {
      for (int i = 0; i < numParams; i++) {
        min[i] = candidate[i]-fabs(candidate[i])*(1.0/65536);
        max[i] = candidate[i]+fabs(candidate[i])*(1.0/65536);
      }
    } else {
      for (int i = 0; i < numParams; i++) {
        min[i] = -1;
        max[i] = 1;
      }
    }
    timeData = new kiss_fft_scalar[bigSize];
    freqData = new kiss_fft_cpx[bigSize/2+1];
    smallFFTR = kiss_fftr_alloc(smallSize, 0, NULL, NULL);
    smallIFFTR = kiss_fftr_alloc(smallSize, 1, NULL, NULL);
    bigFFTR = kiss_fftr_alloc(bigSize, 0, NULL, NULL);
    bigIFFTR = kiss_fftr_alloc(bigSize, 1, NULL, NULL);
  }

  double *getMin() {
    return min;
  }

  double *getMax() {
    return max;
  }

// ___                                                            __ 1
// |  \    |       |       |       |       |       |       |     / |
// |   \   |       |       |       |       |       |       |    /  |
// |    \_ |       |       |       |       |       |       |   /   |
// |      \|__     |       |       |       |       |       |  /|   |
// |       |  -----|_______|___    |       |       |       | / |   |
// |       |       |       |   ----|       |       |       |/  |   |
// --------------------------------x-----------------------x---|---- 0
// 0      1/8     2/8     3/8     4/8     5/8     6/8     7/8 15/16
// |-------------------------------|                       |-------|
//            zeroStarts                                   winStarts
//
// f(x) = 0 if 4/8 < x < 7/8
// f(-x)f(x) + f(-x+1/8)f(x-1/8) = 1 if 0 < x < 1/8

  double costFunction(double *params, double compare, int print) {
    double penalty = 0;
    double accu = params[0]/2;
    for (int i = 1; i < numParams; i += 2) {
      accu += params[i];
    }
    if (print) {
      printf("%.20f", params[0]/2/accu);
      for (int i = 1; i < numParams; i += 2) {
        printf("+%.20fcos(%d pi x)", params[i]/accu, (i+1)/2);
        printf("+%.20fsin(%d pi x)", params[i+1]/accu, (i+1)/2);
      }
      printf("\n");
    }
    if (accu != 0) {
      for (int i = 0; i < numParams; i++) {
        params[i] /= accu;
      }
    }
    const int zeroStarts = 4; // Normally 4
    const int winStarts = 2; // Normally 1
    int i = 0;
    int j = 0;
    freqData[j].r = params[i++];
    freqData[j++].i = 0;
    for (; i < numParams;) {
      freqData[j].r = params[i++];
      freqData[j++].i = params[i++];
    }
    for (; j <= smallSize/2;) {
      freqData[j].r = 0;
      freqData[j++].i = 0;
    }
    kiss_fftri(smallIFFTR, freqData, timeData);
    double scale = 1.0/timeData[0];
    double tilt = 0;
    double tilt2 = 0;
    for (int i = 2; i < numParams; i += 2) {
      if ((i/2)%2) {
        tilt2 += (i/2)*params[i]*scale;
      } else {
        tilt2 -= (i/2)*params[i]*scale;
      }
      tilt += (i/2)*params[i]*scale;
    }
    penalty += fabs(tilt);
    penalty += fabs(tilt2);
    double accu2 = 0;
    for (int i = 0; i < smallSize; i++) {
      timeData[i] *= scale;
    }
    penalty += fabs(timeData[zeroStarts*smallSize/8]);
    penalty += fabs(timeData[winStarts*smallSize/16]*timeData[smallSize-winStarts*smallSize/16]-0.5);
    for (int i = 1; i < winStarts*smallSize/16; i++) {
      // Last 16th
      timeData[bigSize-winStarts*smallSize/16+i] = timeData[smallSize-winStarts*smallSize/16+i];
      accu2 += timeData[bigSize-winStarts*smallSize/16+i]*timeData[bigSize-winStarts*smallSize/16+i];
    }
    // f(-1/8+i)*f(1/8-i) + f(i)*f(-i) = 1
    // => f(-1/8+i) = (1 - f(i)*f(-i))/f(1/8-i)
    // => f(-1/16) = (1 - f(1/16)*f(-1/16))/f(1/16)
    //             = 1/(2 f(1/16))
    for (int i = 1; i < winStarts*smallSize/16; i++) {
      // 2nd last 16th
      timeData[bigSize-winStarts*smallSize/8+i] = (1 - timeData[i]*timeData[bigSize-i])/timeData[winStarts*smallSize/8-i];
      accu2 += timeData[bigSize-winStarts*smallSize/8+i]*timeData[bigSize-winStarts*smallSize/8+i];
    }
    // Between 2nd last and last 16th
    timeData[bigSize-winStarts*smallSize/16] = 1/(2*timeData[winStarts*smallSize/16]);
    accu2 += timeData[bigSize-winStarts*smallSize/16]*timeData[bigSize-winStarts*smallSize/16];
    for (int i = zeroStarts*smallSize/8; i <= bigSize-winStarts*smallSize/8; i++) {
      timeData[i] = 0;
    }
    for (int i = 0; i < zeroStarts*smallSize/8; i++) {
      accu2 += timeData[i]*timeData[i];
    }
    if (print > 1) {
      printf("\n");
      for (int x = 0; x < bigSize; x++) {
        printf("%d,%f\n", x, timeData[x]);
      }
    }
    scale = 1/sqrt(accu2);
    if (print) {
      printf("sqrt(accu2) = %f\n", sqrt(accu2));
    }
    double tSpread = 0;
    timeData[0] *= scale;
    double tMean = 0;
    for (int i = 1; i <= zeroStarts*smallSize/8; i++) {
      timeData[i] *= scale;
      //      tSpread += ((double)i)*((double)i)*(timeData[i]*timeData[i]);
      double x_0 = timeData[i-1]*timeData[i-1];
      double x_1 = timeData[i]*timeData[i];
      tSpread += ((double)i)*((double)i)*(x_0 + x_1)*0.5 - ((double)i)*(2.0/3*x_0 + 1.0/3*x_1) + 0.25*x_0 + 1.0/12*x_1;
      double slope = timeData[i]-timeData[i-1];
      if (slope > 0) {
        penalty += slope+1;
      }
      tMean += x_1*i;
      if (timeData[i] < 0) {
        penalty -= timeData[i];
      }
    }
    double x_0 = timeData[0]*timeData[0];
    for (int i = 1; i <= winStarts*smallSize/8; i++) {
      timeData[bigSize-i] *= scale;
      double x_1 = timeData[bigSize-i]*timeData[bigSize-i];
      tSpread += ((double)i)*((double)i)*(x_0 + x_1)*0.5 - ((double)i)*(2.0/3*x_0 + 1.0/3*x_1) + 0.25*x_0 + 1.0/12*x_1;
      x_0 = x_1;
      tMean += x_1*(-i);
    }
    tMean /= smallSize;
    penalty += fabs(tMean);
    if (tMean > 0) {
      penalty += 1;
    }
    tSpread /= ((double)smallSize)*((double)smallSize);
    if (print) {
      printf("tSpread = %f\n", tSpread);
    }
    kiss_fftr(bigFFTR, timeData, freqData);
    double fSpread = 0;
    x_0 = freqData[0].r*freqData[0].r;
    for (int i = 1; i <= bigSize/2; i++) {
      double x_1 = freqData[i].r*freqData[i].r+freqData[i].i*freqData[i].i;
      fSpread += ((double)i)*((double)i)*(x_0 + x_1)*0.5 - ((double)i)*(2.0/3*x_0 + 1.0/3*x_1) + 0.25*x_0 + 1.0/12*x_1;
      x_0 = x_1;
    }
    if (print > 1) {
      for (int i = 0; i <= bigSize/2; i++) {
        printf("%d,%f,%f\n", i, freqData[i].r, freqData[i].i);
      }
    }
    fSpread /= bigSize; // Includes kiss_fft scaling
    if (print) {
      printf("fSpread = %f\n", fSpread);
      printf("%f,%f,%f\n", tSpread, fSpread, tSpread*fSpread);
    }
    return tSpread*fSpread + penalty;
  }

  double costFunction(double *params, double compare) {
    return costFunction(params, compare, false);
  }

  int getNumDimensions() {
    return numParams;
  }

  ~WinProblem() {
    delete[] min;
    delete[] max;
    delete[] timeData;
    delete[] freqData;
    KISS_FFT_FREE(smallFFTR);
    KISS_FFT_FREE(smallIFFTR);
    KISS_FFT_FREE(bigFFTR);
    KISS_FFT_FREE(bigIFFTR);
  }
};


Is the system $y[n]=x[n]+2=T\{x[n]\}$ an LTI-System?

A system that adds a constant is rarely linear. If you have two inputs $x_1$ and $x_2$, each one gets the constant, so their combination get the constant two times. While the input $x_1+x_2$ gets the constant only once.


What is the minimum sampling rate for QPSK or what is the bandwidth of QPSK signal?

As best as I can tell from the terminology, QPSK a.k.a. 4-PSK or 4-QAM can convey twice the number of bits of information as BPSK in the same bandwidth.

Given a one-sided bandwidth of $B$ and assuming sufficient S/N ratio so that no bits are corrupted, BPSK can transfer bits at a rate of $2B$.  That means, given the same one-sided bandwidth, QPSK can transfer bits at a rate of $4B$.

But the sample rate must be at least $2B$.  If the receiver is perfectly synchronized to the transmitter, there would be one sample per bit with no bit error, if the transmitter sends a bipolar $\operatorname{sinc}(\cdot)$ for each bit.

But for QPSK, each sample is a complex value with a real part and imaginary part, each part containing one bit sample if this is 4-PSK.

So if you use quadrature demodulation to bump the IF or RF down to baseband $i(t)$ and $q(t)$, and sample those signals, you will need one complex-valued sample per pair of bits.


How to get the volume level from PCM audio data?

There are two separate issues behind your question: the duration of the measurement and the algorithm used to do the actual metering. In the audio jargon, the term used for the perceived level measurement is loudness metering.

Time duration

The first issue you need to tackle is the length of the window used to do the metering. Usually, in audio they can be divided in three classes:

Momentary, which gives a loudness level that varies quickly over time (usually a few hundreds of milliseconds)
Short term, which winows of a few seconds
Integrated, which gives a loudness level on the whole audio file from beginning to end

The actual values for the lengths of the windows varies from the specification, but this EBU tech note specifies a duration of 400 ms for momentary, 3 s for short-term and the whole duration of the file for integrated, of course. Moreover, if you do live tracking of the metering, there is the issue of the refresh rate of the measurement as well, but you may want to set this aside in your case.
Finally, it may be useful to note that a rectangular window is used, contrarily to other cases in audio (short-term fourier transform for instance).

Loudness metering

To accurately represent the perceived level, you will need something more complex than a simple RMS measurement. Various methods have been developed over time, from the simple A-weighting mentioned in the comments (which is just filtering before computing the RMS value), to the LUFS (Loudness Units relative to Full Scale) measurement specified in the ITU-R BS.1770-4 recommendation. The latter is a complex algorithm, but it provides methods for measuring perceived loudness in a multi-channel audio file, which is your case.

Regarding your last questions:

Translating the measurement to dB is a good idea to produce a more understandable chart, and most audio products that implement loudness metering give the value in dB.
The windowing method already has inherent smoothing, because giving one value for a window of, say, 400ms time is already a kind of low-pass filtering. Regarding noise removal, most modern algorithms implement gating, i.e. an algorithm that removes all noise below a given level, which means that long periods of silence in the audio file will not lower the value measured when doing integrated measurements. It is then possible to compare the loudness from different sources (between a 2-hours movie and a 3-minutes rock song for instance).

Finally, please note (as it has been highlighted in the comments) that the value will be given in dB relative to full-scale, so it will not reflect the level of the actual audio source, except if you know the calibration of the digital-to-analog converter used for the initial recording.

There are a lot of audio plugins that implement metering, and one of the most recent ones is the Fabfilter Pro-L. You may not need such complexity for your measurement, but it can bee a good idea to take a look at what is done in the industry, the way the value is given and what a good loudness curve looks like.

You can find more comprehensive and extensive sources on loudness metering on the web, and especially on web sites of pro audio manufacturers, for instance TC Electronic or Waves


Finding the deterministic autocorrelation function (ACF) from its power spectrum

The power spectral density is the fourier transformed ACF. So what you have to to is to do the inverse fourier transform of the spectrum. And since the spectrum is a rectangle in frequency, it will be a $\frac{\sin(x)}{x}$ in the time domain.

Additionally:

  I know that the power spectrum is the DTFT of the ACF at 0

The signals power is given by $\phi_{xx}(\tau = 0)$:
$$
\phi_{example}(\tau) = \int_{-\infty}^{\infty} x_{1}(t)* x_{2}(t-\tau)\mathrm{d}x
$$
If you autocorrelate (correlate a signal with itself), your subscripts vanish ($ x_1, x_2$ become $x$) and if you calculate the ACF at $\tau = 0$ you will get the squared power of the signal:
$$
\phi_{example}(\tau = 0) = \int_{-\infty}^{\infty} x(t)* x(t)\mathrm{d}x = \int_{-\infty}^{\infty} x(t)^2\mathrm{d}x
$$

If your signal does not have a finite amount of energy, you will have to use the following definition of the ACF:
$$
\phi_{example}(\tau) = \lim_{T \to \infty} \frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} x(t)* x(t-\tau)\mathrm{d}x
$$


Hessian Matrix. Second partial derivative test

The function is fully approximated if one uses all the derivatives (see Taylor expansion). With using the Hessian only, we can only make a second degree approximation (because it is second derivative matrix), which is geometrically the same as using the second order polynomial.


Computing SSD/SSE using pattern with omitted pixel

By excluding the bottom right pixel, the authors only need to process 8 pixels which can be stored in eight SSE registers (xmm0 through xmm7).
https://en.wikipedia.org/wiki/Streaming_SIMD_Extensions#Registers


Discrete approximation to Gaussian filter

The difference between the two should be slight but the MATLAB one has the property that it preserves the signal energy content. That is a property that you expect to have when using a Gaussian filter.


Discrete wavelet transform; how to interpret approximation and detail coefficients?

Wavelet transforms can be more difficult to interpret than FFT at face value due to the various representations, nomenclature and output formats. I had to study more than 15 resources to get a good sense of the variety and which one is used by Pywavelets (which does not provide much theory or explanation in its documentation).

In order to grasp the meaning of cD and cA coefficients, it is helpful to run through a basic example wavelet transform calculation.
Here's a simple step-by-step calculation of what happens in a multi-level DWT (your example is basically the first level). In this representation, they concatenate cA and cD coefficients side by side.

When running DWT using a specific basis function, the signal is fed through (inner product with) a high pass filter (difference filter) and a low pass filter (smoothing / averaging filter), each of which is unique to the wavelet basis function.

For 'db2', the high/ low pass filtering each has two terms, and occurs with a step size (stride) of 2, therefore, after the filtering is completed, you also get a downsampling by 2 of the original signal. Actual length will depend on the filter length and the signal extension mode.
The high pass filtered result gives you the cD coefficients.
The low pass filtered result gives you the cA coefficients.

In your example, each of the 48 cD and 48 cA coefficients corresponds to 2 of the original 96 data points.
Therefore, to plot the cA and cD coefficients in time, just reduce time resolution by 2.

Hope this answers your question, if not, here're some more resources:

This has some simple equations and intuitive comparisons to FFT
This has some nice graphical representations and color coding to explain multi-level DWT
Another guide with detailed explanation and illustrations


Removing empty/low amplitude part from signal

If you know the frequency of the sine signal  already, you can correlate your [.......p-set...] [SINE-WAVE] [....q-set...] signal with the sine signal. The correlation will give you a peak where sine wave exactly matches with the signal. Since you know the sine signal length you can figure out the beginning of the sine signal in your original signal. I hope this helps you a bit.


Algorithm for Discrete Signal Decomposition

I'd just take a KIS (keep it simple) approach as a first step.

Define your (unknown) signal as
$$ s(t) = k A(t) + m B(t) + n(t)$$
and then just define the error:
$$ e(\tilde{k}, \tilde{m}) = \sum_{\forall t} \left| s(t) - \tilde{k} A(t) - \tilde{m} B(t) \right|^2 $$
and then your estimates $\hat{k}$ and $\hat{m}$ are just chosen as:
$$
(\hat{k}, \hat{m}) = \arg \min e(\tilde{k}, \tilde{m})
$$

The R code below does this, and the output seems to do the right thing. Here $n(t)$ is a Gaussian noise with zero mean and standard deviation of $0.1$.

The actual $k$ and $m$ values are at the green circle, the estimates are the blue diagonal bar beneath it (in the rightmost plot).

This appears to work for any values of $k$ and $m$, including zero.  You can probably extend it for more signals $C(t)$, $D(t)$, etc. too.

R Code Below

# 30067

prototype_length <- 40

A <- rep(0, prototype_length)
A[(prototype_length/4+1):(3*prototype_length/4)] <- sin(2*pi*seq(0,prototype_length/2-1)/prototype_length)

B <- rep(0, prototype_length)
B[(prototype_length/2+1):(prototype_length)] <- `^`(sin(2*pi*seq(0,prototype_length/2-1)/prototype_length),12)*3/8

par(mfrow = c(1,3), pty="s")
plot(A, type="l", col="blue", lwd=5)
lines(B, col="red", lwd=5)
title('Original Signals')

k <- sample(0:10, 1)
m <- sample(0:10, 1)

data <- k*A + m*B + 0.1 * rnorm(prototype_length, 0, 0.1)
plot(data, type="l")
title('Scaled & Added Signals')

err <- array(0,c(11,11))
kvals <- array(0,c(11,11))
mvals <- array(0,c(11,11))
for (khat in 0:10)
{
  for (mhat in 0:10)
  {
    err[khat+1,mhat+1] <- sum('^'(abs(data - khat*A - mhat*B), 2))
    kvals[khat+1,mhat+1] <- khat
    mvals[khat+1,mhat+1] <- mhat
  }
}

image(0:10,0:10,err, xlab="k hat", ylab = "m hat")
ix <- which.min(err)
lines(kvals[ix] + c(-0.2,0.2), mvals[ix]+ c(-0.2,0.2), col="blue", lwd = 10)
points(k, m, col="green", lwd = 10)
title('Error, actual values, and minimum')


What about this: if your three signals are balanced (I mean with equal amplitudes and $2\pi/3$ shifts), you can find the amplitude by summing the squares (check that it doesn't pulse). And from there, the arguments of the three sinusoids ($\omega t+\phi+k2\pi/3$). Actually you solve this geometric problem, where the three projections are given:

Then, repeating the computation at two close instants, you get the frequency. It may be useful to maintain a linear fit on the last few $n$ phase values to improve accuracy and ease handling of phase wraparound.


It really depends on the purpose of your project. If you want to do it from scratches, it will not be easy. However, if you want to do it for just one file, and the goal is to suppress the noise of the drone on that file, I recommend you to download free processing tool: audacity: http://www.audacityteam.org/download/, and use their feature where you can pick a period with a typical noise, analyze it and remove so analyzed noise from the whole file. I think that you may find their source code as well for free, which may include the implementation of the feature as well.
However, if you want to do it from scratches, there is a whole nice theory how to do it, but it is not easy.


Sequential processing of uncorrelated measurements in Kalman Filter

I think it is evident in both cases. It follows directly from
$$ P(A|B,C,D) = P(((A|B)|C)|D)$$
.i.e.
$$P(x[k]\!=\!\hat{x}[k]|\!\hat{z}_1[k]\!=\!z_1[k],\!\hat{z}_2\![k]\!=\!z_2[k])\!=\! P\Bigg(
  \Big(
      x[k]\!=\!\hat{x}[k]|\!\hat{z}_1[k]\!=\!z_1[k]
  \Big)|\!\hat{z}_2[k]\!=\!z_2[k]
\Bigg)$$
As long as you don't do a prediction update it should be fine.
For nonlinear system, you won't get the same exact estimate, because you will end up with a new dynamics, input, measurement and feedback matrices on every "single-measurement" update


This answer does not require any elaborate mathematics. Rather it requires a careful interperation of periodicity of DFT sequences considering their lengths.

Lets display the signal chain below:

$$ x[n] \longrightarrow \boxed{ 32DFT} \longrightarrow  X[k] \longrightarrow \boxed{ \uparrow 2 } \longrightarrow Y[k] \longrightarrow \boxed{ 64-IDFT } \longrightarrow  y[n] $$

Now the signals are:

Time sequence $x[n]$ has a length of $32$ samples defined for $n=0,1,...,31$. For any other values of $n$, we interpret $x[n]$ as the periodic extension by the help of modulus operator as $x[n] = x[(n)_{32}]$. For example $x[35] = x[ (35)_{32} ] = x[3]$.

The frequency sequence $X[k]$ has a length of $32$ samples defined for $k=0,1,...,31$. For any other values of $k$, we interpret $X[k]$ as the $32$ point periodic extension of $X[k]$ similarly.

After the expansion block we have:

DFT frequency sequence $Y[k]$ has a length of $64$ samples defined for $k=0,1,...,63$. For any other values of $k$, we interpret $Y[k]$ as the periodic extension by $64$ samples similarly.

The IDFT signal $y[n]$ has a length of $64$ samples defined for $n=0,1,...,63$. For any other values of $n$, we interpret $y[n]$ as the periodic extension by $64$ samples similarly.

Now, simply write the $64$ point inverse DFT equation to define $y[n]$ :

$$
\begin{align}
y_{64}[n] &= \frac{1}{64} \sum_{k=0}^{63} Y[k] e^{-j \frac{2\pi}{64} k n } ~~~&,~~~ \text{ for } k,n = 0,1,...,63 \\
&= \frac{1}{64} \sum_{k=0}^{63} X[k/2] e^{-j \frac{2\pi}{64} k n } ~~~&,~~~ \text{ for } k,n = 0,1,...,63 \\
&= \frac{1}{64} \sum_{m=0, k=2m}^{31} X[m] e^{-j \frac{2\pi}{64} 2m n } ~&,~ \text{ for } m=0,1,..,31 ~~ n = 0,1,...,63 \\
&= \frac{1}{2} \left( \frac{1}{32} \sum_{m=0}^{31} X[m] e^{-j \frac{2\pi}{32} m n } \right) ~~&,~\text{ for } m=0,1,..,31 ~~~ n = 0,1,...,63 \\
y_{64}[n] &= \frac{1}{2} x[(n)_{32}] ~~~&,~~~ \text{ for } n = 0,1,...,63 \\
\end{align}
$$

The expression in the parenthesis above is simply the $32$ point inverse DFT of the sequence $x[n]$. Note that this sequence has a length of $32$ points and we can express it with the modulus notation $x[n] = x[(n)_{32}]$ as in the last line.

And the last line describes the $64$ point signal $y_{64}[n]$ as a two times repetition of $32$ point $x[n]$, divided by $2$. In other words, $y[0] = x[0], y[15]=x[15], y[40] = x[8], y[60]= x[28]$ etc...


How is the sound converted to matrix in Matlab?

The function audioread doesn't generate any values, it just reads audio samples stored in a file. If you want to generate the sound of a guitar, you need to look into sound synthesis, such as the Karplus-Strong method.


Recommended approach for extracting specific edges from image?

You might find my, or the other answers, from here helpful.

After advice about detecting focus quality of objects in a photo detected using YoloV3

I would have put this in a comment, but I can't include a pic there.

Perhaps you could come up with a better verbal description of how to define the edges of interest.


What is the Reference tap of an equalizer (FIR Filter)

I haven't checked the code, but that parameter refers to the delay of the reference signal when computing the error used for updating the filter coefficients. If $y[n]$ is the output of the linear equalizer, and $r[n]$ is the reference signal (e.g., the training signal when in training mode), the error is computed as

$$e[n]=y[n]-r[n-d]$$

where $d$ is the delay defined by d=RefTap-1. This means that if the channel introduces a pure delay and the parameter RefTap is set to match that delay, the equalizer doesn't need to do anything. Otherwise, the equalizer would need to try to compensate for that delay, which is impossible for a causal filter (because it would try to implement an advance in time). So that parameter is indeed very important in order to achieve the best possible equalization for a given equalizer filter length. Since in practice the channel delay is often unknown, a reasonable choice is to set the parameter RefTap to half the equalizer length. Then, on average, the best use is made of the equalizer's weights. If the performance is not satisfactory, the number of weights needs to be increased (and the parameter RefTap must be changed accordingly).


Stopping Viola Jones Face Detection early

Viola-Jones works in stages.  Only areas of the image that have a high enough "maybe" score get passed to the next stages.  If you implement your own Viola-Jones detector, then you should be able to stop between each stage and apply color/lighting corrections before going on to the next stage.

As an alternative to implementing your own, you might take a look at the openCV implementation.  You might be able to modify it for your needs.

The source is here for the cascade detector, which is the openCV Viola-Jones implementation.  You will need to examine the area after the comment //---------------------------------------- Classifier Cascade

That is where the actual detector is implemented.  I think the interesting places are operator and ocl_detectMultiScaleNoGrouping, and I really think operator is where you would need to start.


OFDM symbol acquisition block in GNU Radio gives rotating constellation even if there is no carrier offset. How is this possible?

Got me at that one! The "OFDM symbol acquisition" block is in fact not from gr-digital (where your other OFDM blocks come frome), but from gr-dtv, where it is used to capture DVB-T signals, if I remember correctly. It might be very DVB-specific!
Let us have a look at the dvbt_rx_8k.grc example from gr-dtv (or, at least, the top half):

So your understanding is absolutely right, OFDM symbol acquisition just does some syncing.

This time, I'll leave a bit more up to your own research (as finishing this up as a nice answer will probably take quite some time), and will concentrate on the steps I took to figure out how the block works. I hope this is most helpful to you, and of course helpful for future answer-seekers trying get the general grasp on how to find out how GNU Radio blocks "tick". Looking at your flow graph, I'm pretty sure you've done most of this already, but as I mentioned, I'd like to illustrate the "investigative principle" for others.

I figured out how the underlying block is called. Typically, the default name of a new block of that type gives that away, or the documentation type. So dtv_dvbt_ofdm_sym_acquisition was what I found.
Then I went to the GNU Radio source repo clone and used git to figure out where to look (I'll annotate the result with a bit of hints):

$ git ls-files|grep ofdm_sym_acquisition
gr-dtv/grc/dtv_dvbt_ofdm_sym_acquisition.xml ##.xml: the file describing how this block looks in GRC
gr-dtv/include/gnuradio/dtv/dvbt_ofdm_sym_acquisition.h ##gr-something/include/: public interface, the one accessible and documented in Doxygen
gr-dtv/lib/dvbt/dvbt_ofdm_sym_acquisition_impl.cc ##_impl.cc,h: the actual implementation, which has all the "ugly" details, private variables and stuff
gr-dtv/lib/dvbt/dvbt_ofdm_sym_acquisition_impl.h

For us, the _impl is of most relevance. So I opened it and looked for a work or general_work method right away (for difference between these, refer to the Guided Tutorials part 4)
The (general_)work method is the central workhorse of a GNU Radio block: GNU Radio's scheduler will call it whenever there are input items to process (given that there's enough output buffer space). So I skimmed it. Pasting it here, trying to delete all lines that I didn't care about at first. Find my "mental" annotations in between.

    dvbt_ofdm_sym_acquisition_impl::general_work (int noutput_items,
…
    {

The usual stuff: input items are in in, output goes to out:

        const gr_complex *in = (const gr_complex *) input_items[0];
        gr_complex *out = (gr_complex *) output_items[0];
…

Ah-Ha! We loop over something, let's keep an eye out for the i.

        for (int i = 0; i < noutput_items; i++) {
…

ok, judging by the variables name this is the case of us not having "locked" on to sync sometime. Let's concentrate on this and ignore the other case for now.

          if (!d_initial_acquisition) {

This looks a lot as if the interesting work is done in ml_sync. d_cp_length definitely points to us having to deal with a system that uses cyclic prefixes.

            d_initial_acquisition = ml_sync(&in[d_consumed], 2 * d_fft_length + d_cp_length - 1, d_fft_length + d_cp_length - 1, \
                &d_cp_start, &d_derot[0], &d_to_consume, &d_to_out);
            d_cp_found = d_initial_acquisition;
          }
          else {
…
            }
          }

I decide to break here and take a look at ml_sync.

The name is pretty revealing already; it's probably going to be a Maximum Likelihood Synchronizer.

So, going up in the source, we find ml_sync. It's pretty well-commented and uses (VOLK = Vector optimized library of kernels) functions that have names like magnitude_squared; I think I'll leave the lecture up to you and give you my quick understanding of it:

volk_32fc_magnitude_squared_32f(&d_norm[low], &in[low], size); calculates the mag² == norm of the input, starting at the low position, saves it in d_norm
The "moving sum" loop calculates the "amount" of correlation and the signal energy in windows of size d_fft_length; that makes sense. A high power signal can have a very high correlation coefficient with something that is not really "similar", so you need to "normalize" later on. Results are stored in d_gamma and d_phi, respectively.
The magnitudes of the windowed-correlation vector elements in d_gamma are stored in d_lambda, and something linearly dependend of d_phi is subtracted from those.
Peak search in d_lambda.
From the peak position, a frequency correction is calculated. From the correlation coefficients, phase-derotation is done.

So what we have here is an autocorrelation-based synchronizer that locks onto the cyclic prefix. Nice!

Comparing this with Schmidl&Cox that was used in ofdm_rx.grc, calculating the autocorrelation is of course pretty CPU-intense, and we need as much Cyclic Prefix as we can get. However, this might be (is) a suitable synchronizer for the typical broadcast channel, where FFTs are long, and so are prefixes, and the trade-off between being fast at decoding packets and not missing a single packets is biased in favor of a RX-maximizing receiver operation characteristic.

To the question why this doesn't work in your setup, I must admit that it's a bit over my head. My intuition tells me that there might be ambiguities due to your 1/4 payload/prefix ratio, and you consistently lock onto the middle of two consecutive frames or so.

The point is: try!


Low frequency and high frequency noise of ecg

The frequency content of baseline wander is usually in the range below 0.5 Hz, while power line noise is somewhere around 50/60 Hz accompanied by some harmonics.


[Python ]How can I improve my 1D Kalman Filter estimate?

I believe ZHUANG's answer is correct, so I went searching for some source code in Python that implemented it.  The closest I could find was a 2D example that uses velocity as well.  The reference is at the top of the listing.  I have reformatted and restructured the code to make it more readable to me and likely more efficient.  I also adapted it to the OPs function.

I am a newbie to Kalman filters.  I am using this as my learning tool.  It can probably be improved upon greatly.  I am posting it hoping that it will be helpful to the OP as well as others.

Ced

# Derived from:
# https://stackoverflow.com/questions/13901997/kalman-2d-filter-in-python
#   and
# https://dsp.stackexchange.com/questions/48343

import numpy as np
import matplotlib.pyplot as plt

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

#---- Set Parameters

        N = 50

        x = np.matrix( '0. 0. 0. 0.' ).T
        P = np.matrix( np.eye( 4 ) ) * .0001000 # initial uncertainty

        R = 0.1**2

#---- Set up Helper Matrices

        F = np.matrix( '''
            1. 0. 1. 0.;
            0. 1. 0. 1.;
            0. 0. 1. 0.;
            0. 0. 0. 1.
            ''' )

        H = np.matrix( '''
            1. 0. 0. 0.;
            0. 1. 0. 0.''' )

        motion = np.matrix( '0. 0. 0. 0.' ).T

        Q = np.matrix( np.eye( 4 ) )

        m = np.matrix( '0. 0.' ).T

#---- Create the Test Case

#        true_x = np.linspace( 0.0, 10.0, N )
#        true_y = true_x**2

#        observed_x = true_x + 0.05 * np.random.random( N ) * true_x
#        observed_y = true_y + 0.05 * np.random.random( N ) * true_y

        sz = (N) # size of array

        true_x = np.linspace( 0.0, 10.0, N )
        true_y = true_x + 2 * true_x**2

        observed_x = true_x
        observed_y = true_y + np.random.normal( 0, 7.5, size=sz )

#---- Run the Sequence

        kalman_x = np.zeros( N )
        kalman_y = np.zeros( N )

        for n in range( N  ):
            m[0] = observed_x[n]
            m[1] = observed_y[n]
            x, P = kalman( x, P, m, R, motion, Q, F, H )
            kalman_x[n] = x[0]
            kalman_y[n] = x[1]

#---- Plot the Results

#        plt.plot( observed_x, observed_y, 'ro' )
#        plt.plot( kalman_x, kalman_y, 'g-' )
#        plt.show()

        plt.figure()
        plt.plot( true_x, observed_y,'k+',label='noisy measurements')
        plt.plot( true_x, kalman_y,'b-',label='a posteri estimate')
        plt.plot( true_x, true_y,color='g',label='truth value')
        plt.legend()
        plt.xlabel('t')
        plt.ylabel('f(t)')

        plt.show()

#=============================================================================
#def kalman( x, P, measurement, R, motion, Q, F, H ):
def kalman( x, P, m, R, motion, Q, F, H ):
        '''
        Parameters:
        x: initial state
        P: initial uncertainty convariance matrix
        measurement: observed position ( same shape as H * x )
        R: measurement noise ( same shape as H )
        motion: external motion added to state vector x
        Q: motion noise ( same shape as P )
        F: next state function: x_prime = F * x
        H: measurement function: position = H * x

        Return: the updated and predicted new values for ( x, P )

        See also http://en.wikipedia.org/wiki/Kalman_filter

        This version of kalman can be applied to many different situations by
        appropriately defining F and H
        '''

#---- UPDATE x, P based on measurement m
#            distance between measured and current position-belief

        y = m - H * x

        S = H * P * H.T + R  # residual convariance
        K = P * H.T * S.I    # Kalman gain

        x = x + K * y
        I = np.matrix( np.eye( F.shape[0] ) ) # identity matrix
        P = ( I - K * H ) * P

#---- PREDICT x, P based on motion

        x = F * x + motion
        P = F * P * F.T + Q

#---- Exit

        return x, P

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


Why are there so many windowing functions?

Aside from reduction of spectral leakage, there is a one major trade-off to be made when choosing a window function. Below you can see a figure with various parameters. Two of them are most important:

Main Lobe Width
Peak Side Lobe Level

Width of the main lobe affects resolution of your analysis. I am sure you know that multiplication in time domain is convolution in frequency domain. By default you always use rectangular window, which has frequency response in a form of $\mathrm{sinc}$ function.

For example it might be the case that two pure sinusoids (Fourier Transform being Dirac deltas) with frequencies being very close to each other, peaks can be very hard to distinguish and thus become one. On figure below I compared 4 window functions for two sinusoids separated exactly by 3 frequency bins. You can see how wide is the main lobe of Nuttall window and how it affects the separation, whereas Rectangular window is the window with narrowest possible width of main lobe (i.e. lowest Normalized Equivalent Noise Bandwidth - NENBW)

On the other hand level of first side lobe affects the ability to recognize adjacent frequencies when spectral leakage is present and side lobes of frequency with high amplitude are almost "covering" other frequency next to it with much lower amplitude:

And here is the trade-off... You can't get very low level of side lobes and very narrow main lobe, i.e. it's not possible to have window function with main peak like in rectangular window, and very low level of side lobes as in Nuttal window. Here is the comparison of some windows:

Apart from parameters mentioned above, there is more features such as:

Side Lobe Roll-off rate (how fast side lobes decrease - for example see Bartlett and Nuttal above)
Existence of analytical frequency response equation (useful when performing theoretical calculations)
Summing of window samples to constant value in time domain when doing overlap (for example Hann window with $33\%$ of overlap)

This is extremely useful when you want to perform STFT and then reconstruct time domain signal.

So now to answer your question. Why do we have so many window functions?

Depending on task whether it is simple FFT analysis of: closely spaced sinusoids, white noise, impulsive sounds, or even FIR filter design, different type of window is recommended. Obviously in 90% of cases people are using Hamming/Hanning window or no window at all. Nevertheless for some specific signals you might find that extracted features depending on window choice tend to be more discriminative.

For more info about Window Functions please refer to:

  F. J. Harris - On the Use of windows for Harmonic Analysis with the
  Discrete Fourier
  Transform

  G. Heinzel, et al - Spectrum and spectral density estimation by the
  DFT, including a comprehensive list of window functions and some new
  flat-top windows

  M. Cerna, A. F. Harvey - The Fundamentals of FFT-Based Signal
  Analysis and
  Measurement

  H. A. Gaberson - A Comprehensive Windows
  Tutorial


Unclear inverse Z-transform of $G(z) = \frac{1-p}{z-p}$

As we've already seen in your previous question about that paper, the authors of that paper are not very much into details when it comes to math and signal processing. What Eq. $(10)$ shows is definitely not the inverse $\mathcal{Z}$ transform of

$$G(z)=\frac{1-p}{z-p}\tag{1}$$

which is

$$g[n]=(1-p)p^{(n-1)}u[n-1]\tag{2}$$

as you've figured out by yourself. What Eq. $(10)$ in the paper does show is the step response, given by

$$s[n]=\sum_{k=-\infty}^{n}g[k]=(1-p)\sum_{k=1}^np^{(k-1)}=1-p^n,\quad n\ge 0\tag{3}$$

So if you apply a step $\tilde{\mu}\cdot u[n]$ at the input, the system's response is

$$\mu[n]=\tilde{\mu}\left(1-p^n\right),\quad n\ge 0\tag{4}$$

as given by Eq. $(10)$ in the paper.


What are some of the disadvantages of parametric PSD estimation methods?

Can we trust a spectrogram obtained this way even though the signal can't be modelized or doesn't fit the correct framework of the used method ?

Short answer: anything based on a model will work only when applied to the model. Not quite sure which other answer you were expecting here!

Of course, some "model mismatches" are more benign than others: For example, some classes of estimators are very robust against noise not fitting the original model, but break down when the assumed number of signal components isn't what was specified (example: ESPRIT). Others fail the other way around (Bartlett's method will not like you having negative SNR).


Why these cross-connections in the allpass filter?

This is called a lattice structure implementation of all-pass filters.

$k_1$ and $-k_1$ are scalar gains. The input to the delay element is $$a[n]=x[n]-k_1a[n-1]$$ taking the $z$-transform, the transfer function is
$$A(z)(1+k_1z^{-1})=X(z)\Rightarrow A(z)=\frac{X(z)}{1+k_1z^{-1}}$$ and the output is
$$\begin{align}
y[n]&=k_1a[n]+a[n-1]\\
\end{align}$$
Thus
$$\begin{align}
Y(z)&=A(z)(k_1+z^{-1})\\
&=X(z)\frac{k_1+z^{-1}}{1+k_1z^{-1}}
\end{align}$$
$$H(z)=\frac{Y(z)}{X(z)}=\frac{k_1+z^{-1}}{1+k_1z^{-1}}$$
So it is an all-pass filter ($|H(e^{j\omega})|=1,\ \forall \omega$) with mirrored coefficients in numerator and denomerator.


Can spectral density be a complex quantity?

and the same goes with the spectral density (Fourier transform of the auto-correlation function).

No, that's not the case.

Since the autocorrelation is a hermitian symmetric function for any $S$, its Fourier transform is always real.

  If the above analysis valid? as the process is not wide sense stationary.

If the process is not WSS, then you can't just proclaim $E[S^*(t)S(t-\tau)]$ to be dependent on only one variable (usually, $\tau$), and hence, a (1D) Fourier transform doesn't make much sense.

  If the analysis is not valid, is there some way to handle this kind of situation?

Depends! You might want to define/find coherency times and do Short-Time Fourier Transforms within that.

Your system, in fact, is just a phase shifted impulse response – as such, a phase-delay representation, which might be derived from a Frequency Shift-Delay plane, might be more helpful in analyzing things. You'll find such a Frequency Shift-Delay plane in what is called scatter function in wireless communications, representing the Doppler and path coefficients of a wireless channels in motion.

But in your case: Is trying to understand PSD or PSD-equivalents really useful? Don't you just want to build a parametric estimator for $\mu_G$, $\sigma_G^2$ and $\theta$ instead?


Variance of the Discrete Gabor Transform for practical implementation

What I look for in such parameter choice problems is empirical experiments. Whatever works best for my application is the correct value to use.

You can develop a parameter search algorithm looking for the optimal values in some training data set. Keep in mind that what works for one application may not perform well in others.

I'm not aware of mathematical arguments estimating  gabor's $\theta$. These estimates are typically conservative anyway.


Negative exponential signal's energy and power

An ideal exponential signal $x(t)=e^{at}$ , which extends from $-\infty$ to $\infty$ has infinite energy and infinite power, as for real $a >0$ ( and similarly for real $a < 0$) you have

$$
\mathcal{E}_x = \int_{-\infty}^{\infty} e^{2at} dt = \lim_{t \to \infty} \frac{1}{2a} e^{2a t} \to \infty
$$

and

$$
\mathcal{P}_x = \lim_{T \to \infty} \frac{1}{T} \int_{-T/2}^{T/2} e^{2at} dt = \lim_{T \to \infty} \frac{e^{a T}}{2aT}  \to \infty
$$

Therefore It's neither energy nor power signal. Note that it's not a practical signal and exists only for mathematical purposes.


Ktuncer, there are a number of methods you can use here. One method that I would recommend is to use a Discrete Wavelet Transform, (DWT), and in particular, look at the Daubechies Wavelet. I would pick, say, Daub-14 / Daub-Tetra.

Basically what you really need to do is 'trend' your signal, and then from there, do a min or max pick. That will get rid of your outliers. A daub-14/daub-tetra wavelet transform can help you do this, and this helps especially since you do not know the nature of your signal. (Using daub-14, you can accurately represent polynomial signals of degree 14/2 = 7, and it looks like you wont need more than that).

The computation of this wavelet transform essentially 'compresses' your energy into a few indicies. Those indicies represent weights on basis vectors. The rest of the weights will (ideally) be near zero. When you have noise in your signal (as you do), those weights that were normally zero have some weights now, but you can simply zero them out and 'denoise' your signal. Once that is complete, you can then do a simple max/min detect.

There are more details involved, you may send me an email if you like to discuss how to implement it. I have done similar work on this before.

EDIT: Here are some images illustrating Daub-Tetra Denoiser:


Will binarization / thresholding of image yield better results when performed in a linear color space?

Disclaimer: What I am about to say is a general rule of the thumb, there are exceptions.

In short: Linear space is more useful for computer vision while non-linear space is more useful for showing the final image to humans, without outputing a quantitative analysis.

This highly depends on your problem.

Linear space is more related to physical measurements. Thus, thresholding in it is more natural when you are dealing with measurements such as: counting instances, location of objects, etc. Noise behaves in a predictable way - variance of the signal is linearly dependent on the mean value of the signal. There are many other properties that make your algorithms more straightforward and easier to implement. Among other things, when you are selecting a threshold, you can use the linearity of the signal.

On the other hand, the non-linear space is more related to human perception of light intensity. Thus, it is more suited for image processing problems in which there is no measurement, but rather showing a good image to the observer, or measuring something in a way that humans do.

Examples:

Linear space -
Counting cells in microscopic images, finding defects in a wafer, tracking cars on a road. (Some image processing operations must be done in linear space due to the physics of the process involved, but the final image is almost always in non-linear domain)

Non-Linear space - Performing image processing pixel crunching operations to show a nice image on your tablet, attempting to find defects on a carpet texture, performing image segmentation of medical images in order to assist a human interpreting a scene.


How to Combine / Fuse 2 Least Squares Estimates?

I'm unsure if this is the answer you are looking for, but why not save and share $|H_i|^2$ in addition to the least squares estimates?

If you have $w_1^*=\frac 1 {|H_1|^2} H_1^t d_1$, $w_2^*=\frac 1 {|H_2|^2} H_2^t d_2$, and also know $|H_1|^2$ and $|H_2|^2$, you get the total least squares estimate with:

$w^*=\frac 1 {|H_1|^2+|H_2|^2} (H_1^t d_1 + H_2^t d_2)=\frac{|H_1|^2}{|H_1|^2+|H_2|^2}w_1^*+\frac{|H_2|^2}{|H_1|^2+|H_2|^2}w_2^*$.


Harmonic Product Spectrum algorithm. I dont understand one step

"8000 / buffer_size"  is frequency resolution, decide by  sample length or "buffer_size",

you know the sampling interval is "1/8000", the whole time length 'T' of the sequence is "buffer_size/8000", then '1/T' is frequency resolution,

in your code,if 'k' is 0, represent 0hz,is the signal's dc part,if 'k' is 1,represent '1/T' or "8000 / buffer_size" hz, then you can see 'k' correspond with frequency bin "k*8000 / buffer_size" .


Why is this DFT of a real symmetric signal resulting in complex valued coefficients?

You defined the signal vector as x = [1 2 3 2 1]. Since the DFT is defined by

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

the command fft(x) computes the DFT of the signal

$$x[n]=\delta[n]+2\delta[n-1]+3\delta[n-2]+2\delta[n-3]+\delta[n-4]$$

This signal is not symmetrical with respect to $n=0$, and it is not equal to the signal you computed the DTFT of.

If you want to compute the DFT of the signal

$$x[n]=\delta[n+3]+2\delta[n+1]+3\delta[n]+2\delta[n-1]+\delta[n-3]$$

you have to periodically continue it in the interval $n\in [0,N-1]$ (with $N=5$):

x = [3 2 1 1 2]

This results in a real-valued and symmetrical DFT, as expected:

X = fft(x);
X =

   9.00000   2.61803   0.38197   0.38197   2.61803


Are there any recursive online max/min filters for time-series

I think you are looking for a "peak" detectors. They are not filters in the conventional sense. Typically the output is increased when the incoming sample is either above a certain threshold or above the current output and decreased if it's smaller. The increase and decrease are often just simple exponential with time constants determined by your application. The increase time constant is typically called "attack" time and the decrease time constant is "decay" time.


You could form a statistical test, based on the autocorrelation of the potentially-white sequence. The Digital Signal Processing Handbook suggests the following.

This may be implemented in scilab as below.

Running this function over two noise sequences: a white noise one, and a lightly filtered white noise one, then the following plot results.  Script for generation of each realization of the noise sequences is at the end.

The mean of the statistic for the white noise is 9.79; the mean of the statistic for the filtered noise is 343.3.

Looking at a chi-squared table for 10 degrees of freedom, we get:

and we see that there is no significance level at which 9.79 (in the table) that the white noise isn't white. We also see that the value of 343.3 is very likely to be non-white (comparing it to the 23.2093 value in the $p=0.01$ significance column).

function R = whiteness_test(x,m)
    N = length(x);
    XC = xcorr(x);
    len = length(XC);
    lags = len/2+1 + [1:m];
    R = N*sum(XC(lags).^2)/XC(len/2+1).^2;
 endfunction

X = rand(1,1000,'normal');
Y = filter(1,[1 -0.5],X)
R = [R; whiteness_test(X,10)];
R2 = [R2; whiteness_test(Y,10)];


DCT of 4x4 checkerboard matrix - what's the correct result?

As far as I can see, there are two little mistakes. First of all, you should use dct2() instead of dct(). And second, you shouldn't use the 'transpose' function because you want the DCT matrix not only transposed but also conjugated. So you should do the following:
$$\tt{T=dctmtx(4); DCTB = T*B*T'}$$
And this should equal $\tt dct2(B)$.

EDIT: As correctly pointed out by Peter K., the DCT matrix is real-valued, so indeed the 'transpose' operation and the ' operation give the same result. So I only noticed ONE little mistake ...


closed-form analytical expression for DFT coefficients

The DFT of a vector of size L is L samples of the Fourier transform of the signal:

$$
X[k] = X(e^{j\theta})|_{\theta = 2\pi k/L}
$$

So now we need to calculate the FT of $x[n]$:

$$
X(e^{j\theta}) = \sum_{n=0}^{n=M-1} e^{-j n \theta}
$$

This is a geometric series which results in:

$$
X(e^{j\theta}) = \frac{1 - e^{-jM\theta}}{1 - e^{-j\theta}}
$$

Evaluate on the DFT samples to give:

$$
X[k] = \frac{1 - e^{-j M 2 \pi k/L}}{1 - e^{-j 2 \pi k / L}}
$$

You can extract half phase from top and bottom to give:

$$
X[k] = e^{-j(M-1)\pi k/L} \frac{\sin(M \pi k/L)}{\sin(\pi k/L)}
$$


Finding the convolution $a^n u[n]$ with $b^n u[n]$

Just do it the hard way by writing out the individual formulas for the general convolution sum
\begin{align}
y[n] &= \sum_{k=-\infty}^\infty a[k]b[n-k]\\
&= \sum_{k=0}^\infty a[k]b[n-k] &\scriptstyle{\text{because }a[k]=0~\text{whenever } k < 0,}\\
&= \sum_{k=0}^n a[k]b[n-k] &\scriptstyle{\text{because }b[n-k]=0~\text{whenever } k > n,}
\end{align}
which gives
\begin{array}{rclcl}
y[0] &= &a[0]b[0] &= &1\\
&&&= &\displaystyle\frac{b^1-a^1}{b-a}\\
y[1] &= &a[0]b[1] + a[1]b[0] &= &b+a\\
&&&= &\displaystyle\frac{b^2-a^2}{b-a}\\
y[2] &= &a[0]b[2] + a[1]b[1] + a[2]b[0] &= &b^2 + ba +a^2\
&&&= &\displaystyle\frac{b^3-a^3}{b-a}\\
&&&= &\displaystyle\frac{b^3-a^3}{b-a}\\
\vdots &= &\ddots~~~\scriptstyle{\text{Hey, Ma! I think I see a pattern developing here}}\\
\\
y[n] &= &a[0]b[n] + a[1]b[n-1] + \cdots + a[n]b[0]
&= &b^n + b^{n-1}a + \cdots + ba^{n-1} + a^n\\
&&&= &\displaystyle\frac{b^{n+1}-a^{n+1}}{b-a}
\end{array}
if you can recall the formula for the sum of a geometric series.


Is Savitzky-Golay filter a special case of an FIR lowpass filter?

Yes. The Savitzky-Golay filter is a weighted sum of a finite number of data samples. So it is indeed an FIR filter (albeit, not causal, because it uses samples from the past and the future). And it is indeed a low-pass filter (smoother) as seen from the frequency response of the filter on the same page.


Reconstruct FIR frequency response from the coefficients

You can try to window the filter with some window function, i.e. Gaussian, to get rid of small coefficients (taper them). Although it won't really work very well and you might really want to think about redesigning your filter - it's the best solution in my opinion.

Here is an example for a bandpass filter, showing that it is not really working well and is affecting the frequency response. Everything is done in Python so presumably  it's meeting your requirements of not spending a single penny.

Code is as follows:

import numpy as np
import matplotlib.pyplot as plt

# Some sampling frequency
fs = 48000.0
# Size of FFT analysis
N = 1024

def fir_freqz(b):
    # Get the frequency response
    X = np.fft.fft(b, N)
    # Take the magnitude
    Xm = np.abs(X)
    # Convert the magnitude to decibel scale
    Xdb = 20*np.log10(Xm/Xm.max())
    # Frequency vector
    f = np.arange(N)*fs/N

    return Xdb, f

if __name__ == "__main__":

    # FIR filter coefficients
    b = np.array([-0.0159147603287189,  0.000745724267485,
                   0.0404251831063147,  0.019042013872129,
                  -0.0644535904607569, -0.054523709591490,
                   0.0787623967281351,  0.105430811100048,
                  -0.0645610841865355, -0.148306808873938,
                   0.0257418551415616,  0.166568575042643,
                   0.0257418551415616, -0.148306808873938,
                  -0.0645610841865355,  0.105430811100048,
                   0.0787623967281351, -0.054523709591490,
                  -0.0644535904607569,  0.019042013872129,
                   0.0404251831063147,  0.000745724267485,
                  -0.0159147603287189])

    # Window to be used
    win = np.kaiser(len(b), 15)
    # Windowed filter coefficients
    b_win = win*b

    # Get frequency response of filter
    Xdb, f = fir_freqz(b)
    # ... and it mirrored version
    Xdb_win, f = fir_freqz(b_win)

    # Plot the impulse response
    plt.subplot(211)
    plt.stem(b, linefmt='b-', markerfmt='bo', basefmt='k-', label='Orig. coeff.')
    plt.grid(True)
    plt.hold(True)
    plt.stem(b_win, linefmt='r-', markerfmt='ro', basefmt='k-', label='Windowed coeff.')
    plt.plot(win*b.max(), '--g', label='Window function')

    plt.title('Impulse reponse')
    plt.xlabel('Sample')
    plt.ylabel('Amplitude')
    plt.legend()

    # Plot the frequency response
    plt.subplot(212)
    plt.plot(f, Xdb, 'b', label='Orig. coeff.')
    plt.grid(True)
    plt.hold(True)
    plt.plot(f, Xdb_win, 'r', label='Windowed coeff.')

    plt.title('Frequency reponse')
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude [dB]')
    plt.xlim((0, fs/2)) # Set the frequency limit - being lazy
    plt.legend()

    plt.show()


Definition of Ideal Low pass filter (Time Continous)

An ideal low-pass filter has frequency response $H(f)$ equal to the brickwall function: $$H(f)=\mathbb I[|f|<f_c]$$ where $\mathbb I[s]$ is the indicator function, which is equal to 1 if the statement $s$ is true and 0 otherwise, and $f_c$ is the cutoff frequency. Since this frequency response is real, its phase is zero.

This ideal filter is not stable and non-causal. It also has zero delay; it corresponds to $t_0=0$ in the first definition you list above. Now consider the case where we want an ideal LPF, but with some delay $t_0>0$. If the filter input is a sinusoid of frequency $f_0<f_c$, then the filter output is  delayed a time $$-\frac{-2\pi f_0t_0}{2\pi f_0}=t_0.$$

In other words, the first filter in your question is ideal with a delay $t_0$ that you can measure or specify; the second filter is ideal with delay equal to 0. As pointed out in another answer, this helps to make the filter "more causal" -- it allows you to implement a filter that approximates the ideal.

See here for another explanation of phase and delay.


My understanding of your question is that you have a small misunderstanding here. In the definition it is not said time and space. A signal can vary in time or space. Some signals vary with time, as your example human voice that varies over time in air pressure (or equivalently voltage), some vary with space, like image and some vary with both time and space such as video. I feel,  your misunderstanding in human voice example is that you believe there should be an space dimension too in the game.


How to pick coefficients for Fractional Delay Filters?

I usually do this by creating a low pass filter that entirely passes through the signal that I want to delay.  I create the LPF "manually" by creating a windowed sinc function.  Something along the lines of-

filt = sinc(-80:.8:80);
filt = filt .* hamming(length(filt)).';

This gets you a filter that passes about 80% of the nyquist region (the 80% is set by the .8 increment in the sinc indexing).  To get a fractional delay you simply include the fractional delay in the calculation of your filter.  Adding a constant to the sinc indexing delays the filter by $\frac{constant}{increment}$ samples.  For instance, in the case of the filter above, adding 0.1 would create a fractional delay of $\frac{0.1}{0.8}$ samples, or $0.125$ samples.

filt = sinc((-80:.8:80) + .1);
filt = filt .* hamming(length(filt)).';

This works because filtering is a linear operation and delaying the filter is equivalent to convolving your original LPF with a fractionally delayed delta function.  Since linearity implies that the associative property holds, convolving your signal with the delayed LPF is equivalent to convolving your signal with a fractionally delayed delta function and then filtering it with the "undelayed" LPF.

Note: When creating the "undelayed" LPF you must make sure that the indices are symmetric.  The easiest way to do that is to make sure the start and end points are integer multiples of your indexing increment.  That is why I used +/- 80 for an increment of 0.8.


Matlab: Estimating power spectral density of an experimental data?

MATLAB has built-in functions taking care of the steps you mention in 2. You can check out the pwelch function (here) which uses Welch's method for PSD estimation.(here)
The choice of segment length, number, overlap and windowing function presents a trade-off between bias and variance.


The area under an impulse function is one. So if you multiply a continuous time waveform with an impulse, the result will be a weighted impulse with the area equal to whatever the waveform was at that moment in time.


My answer is related to this question
Why root raised cosine filter can eliminate intersymbol interference (ISI) ?

A classical system with Nyquist pulse $p(t)$ is:

The equivalent baseband signal has bandwidth limited in $[-1/2T,1/2T]$ and a sampling rate $F_s > 1/T$ is enough to avoid aliasing.

We normally sample at rate $t=kT$. I wonder what would happen if suddenly I sample faster than the normal sampling rate, for example at $t = kT/2$ instead of $t = kT$ ?

In the frequency domain, the composite filter $g(t) = p(t) \star p^*(-t)$ has Fourier transform $G(f)$ satisfies $\sum_k G(f-k/T) = T, \forall f$. Double the sampling rate $\sum_k G(f-2k/T)$ is not constant any more. This is true at least for typical pulses such as sinc or raised cosine.
In the time domain, thank Dan Boschen for the visualization below, the sample at impair $k$ of $kT/2$ is contributions of two adjacent symbols.

Thus, is it true that for free-ISI-Nyquist-pulse-shaped signal, the largest sampling rate is the value of nominal (Nyquist) bandwidth $1/T$ ?

One can say that we still have correct data by taking only one per two samples (downsampling). What if I sample at rate $F_s > 1/T$ that satisfies Shannon-Nyquist sampling theorem to avoid aliasing ? We still have ISI.

OK another one can say that the criterion $F_s > 1/T$ is in the sense that the entire information after being sampled is intact, just take a baseband filter to get back the continous time spectrum and then sample it at the correct sampling rate at $t=kT$. Have I answered my questions by myself ?


Situations where complexity is too big to exceed linearity and gaussianity

Well, let's look at the two issues: 1) linearity and 2) Gaussianity.

Linearity

If you're imaging moving 3D objects (people) with a single camera, then you're working with a 2D projection of those 3D objects.  That dimensionality reduction can cause non-linearities to appear.

Take a 2D to 1D example:  an object moving in a circle in 2D.  The object is moving with constant velocity and constant acceleration.  All nice and linear.  Project that object onto a 1D surface (a line) and the velocity and acceleration are now (co)sinusoidal.

One of the underlying assumptions in Kalman filtering (well, the sort usually presented) is that the system you're modeling is linear.  When dimensionality reduction (3D $\rightarrow$ 2D) is involved, this may not be true.

Gaussianity

When a system is linear, any noises tend to be additive.  When Gaussian noise is added to Gaussian noise (or, for enough samples, noise of just about any distribution is added to noise of just about any distribution), the result is also Gaussian.

If the system is non-linear, then the operations involved cause Gaussian noise to be transformed into potentially non-Gaussian noise.

For example, if you divide one Gaussian variable with another, the distribution of the result is Cauchy.

What's so bad about that?  The mean (average) of a Cauchy random variable does not exist.  The variance of a Cauchy random variable does not exist.

The Kalman filter is all about tracking means and variances of variables in the system.  If the system is non-linear, it's possible that the relevant quantities' means and variances do not exist, so application of a Kalman filter is meaningless.


Texture-like measures for quantifying density of data in binary images

One of the most obvious that comes to mind is MPEG-7's Edge Histogram Descriptor.

In this, you divide the image into blocks and consider the edge pixels across different angles forming the complete local histogram. This is organized the same way as a Gabor transform features at angular orientation at different angles. This is fairly a very good indicator to classify and search pictures with similarity in their structures and hence edges.

See this reference 2 which is a great introduction to MPEG 7 color and texture. See Reference 3 for a more detailed coverage for the Local Edge descriptors.

Thomas Sikora, The MPEG-7 Visual Standard for Content
Description—An Overview IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY, VOL. 11, NO. 6, JUNE 2001
B. S. Manjunath, Jens-Rainer Ohm, Vinod V. Vasudevan, and
Akio Yamada Color and Texture Descriptors IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY, VOL. 11, NO. 6, JUNE 2001
Dong Kwon Park, Yoon Seok Jeon, Chee Sun Won Efficient Use of Local Edge Histogram Descriptor Proceeding MULTIMEDIA '00 Proceedings of the 2000 ACM workshops on Multimedia Pages 51-54


Calculate Collision force using Audio Signal

Answer:
Whenever the impact happens, that portion of audio has to be detected( if you want to automate it) or else you can identify those sections manually.

To automate it you need to identify characteristics of the audio during the impact, for ex: the frequency contents, or time spread. A good approach to combine these two characteristics is to either do a wavelet decomposition (morlet wavelet) or do a short time Fourier transform.

Once you have a method to identify these regions then you would need a reference (may be using an already available audio meter) to map these signals to audio (in DB), this would depend on how you quantize the analog audio, Sampling rate, FFT size etc. From this identified reference you can then map other audio impacts to suitable audio decibels.


Why is the Fourier transform valid only for absolutely integrable signals?

Fourier on exponential functions can indeed make sense. In theory, under specific conditions. And  in practical DSP, not likely.

There are so many conditions under which Fourier transformations exist that whole books are devoted to them, like   D. C. Champeney, A Handbook of Fourier Theorems,  Cambridge University Press, 1987. And the domain remains an open topic, in other words not all conditions are known, under which Fourier series are unique for instance.

Simply put: if $x(t)$ is absolutely integrable ($x\in L_1$ space), so is $x(t)e^{-jwt}$, by the virtue of the dominated convergence theorem, since $|x(t)e^{-jwt}|\le | x(t)|$. As such, this is a very direct condition to satisfy.
By the way,
note that $e^{-jwt}$ IS NOT absolutely integrable, an interesting paradox.

In signal processing applications and in statistics, as orthogonality or energy preservation is important, it is common to further restrict the domain to $x\in L_1 \cap L_2$. And this is not necessary for discrete signals, as  $\ell_2 \subset \ell_1$.

As exponential functions $\exp(-at)$, $a\neq 0$, are not integrable on the whole $\mathbb{R}$-line, they should be truncated or windowed to be used in practice.

And if you REALLY want to define it, you can go to some distribution theory, see for instance Fourier transform of (real) exponential, where you can restrict the space of tests functions in the spirit of windowing.


You can do a lot for "free" by working with python/numpy.  You can work with .wav files.  Most computers have some kind of sound proceseeing hardware for audio I/O, but you may want to get a sound card if you're computer doesn't already have this.  You can record audio for processing using Windows Sound Recorder.  Numpy has tools for importing and exporting .wav files.

If you really want to learn about DSP coding, get a development Kit from TI or Analog Devices.  These kits include a board with DSP, computer I/O (often via USB) and development tools for generating and debugging code.  Analog Devices seems to have the least expensive development kits.  Try something in the EZ-KIT family from Analog Devices.


Extracting one single hash per video and expecting it to work for subsegment matching is unrealistic.

The approach - which is the one used in audio identification systems but could be extended to video - is to extract many hashes, one for each small segment of content (say a few seconds of video). The matching process consists in matching your "query" sequence of hashes to the sequence of hashes extracted from your video. This is quite similar to a text indexing/retrieval problem, and traditional data structure for such problems (like inverted indices) would work. You can probably build your implementation on top of an existing text indexing system like Lucene/Solr. See for example, the echoprint server code, built on Solr.

This addresses the sequence aspect. Remains the problem of generating robust hashes for images. What you are doing is probably going to still be sensitive to minor differences in compression/resolution - You'll need a very rough quantization of a small number of DCT / wavelet coefficients for it to be robust.


Defending Continuous values using Error Correction Codes

In communications engineering this problem is called "Unequal error protection" and you can find many references on google. It can be a very involved topic.

However, a simple method is to group your bits of each sample into N groups (e.g. if you hve a 16-bit sample, you have N=4 groups, of 4 bit each). Then, you can encode each group with a different code rate (e.g. by puncturing the bits) and form a frame from these bits. Since the code performance increases with its length, you can increase performance when you sacrifice latency. Though, this really depends on your particular application.


Well, the Angle modulation basically means adding a phase function to the carrier that depends on time.

Let's call that instantaneous phase deviation function as:

$$
\Phi(t)
$$

and our modulated signal is

$$
m(t) = A(t) * cos(\omega t + \Phi(t)).
$$

As you can see, if we eliminate the amplitude modulation (constant A(t)), we have a phase shift that varies with time in the modulated signal formula. And hey, why the heck we don't use this phase shift to transmit our information? It can be noise resistant and easy to demodulate at the same time! It's a terrific idea!

So, if we come to your question, the answer is obvious. We are modulating our information to carrier's instantaneous phase function, therefore it's angle.

If you are looking for further information and modulation techniques, you can take a look at frequency shift keying or phase shift keying modulations. Angle modulation is not limited to these modulation techniques, however they are easy to understand and give you a deeper understanding about angle modulation.

As always, have a nice day.


Is the path loss exponent $\alpha$ is always 2 in friis equation model?

Path loss exponent is not always equal to 2, only in the free space scenario is this true. In general the received power $P_r$ is proportional to $d^{-\alpha}$, where $\alpha$ is somewhere around 2-6 depending on the situation. The choice of the path loss exponent should match the scenario you're trying to simulate.

For example, for communication at a high altitude with little obstructions then $\alpha$ is probably close to 2 since the situation is approximately free space condition.

For communication in a building, within the same floor, then $\alpha$ may be around 2-3. This is because walls other obstructions attenuate the signal power further than the free space condition.

Things can get even worse if we again consider communication in a building but now if the transmitter is on the first floor and receiver is on the third floor, then $\alpha$ could be as high as 4-6 depending on construction of the building, etc.


The plot of instantaneous power of the Dirac function

The problem becomes as I try to calculate the power in time domain, $P(t) = \left\lvert\delta(t)\right\rvert^2$

As you figured out yourself, that's impossible: the power of the Dirac impulse isn't defined, it diverges (goes to infinity).
Hence, you especially can't plot the power "correctly"; it's zero anywhere but at $t=0$, where it diverges.
But quite honestly, this should come as no surprise: you can't correctly plot $\delta(t)$, either (the value at $t=0$ isn't defined just as much as its square isn't).
So, whatever question can be answered by that plot: I'm afraid you'll have to find a "better" question.


What is the best codec for streaming of audio in this case?

Opus would be a good choice. It claims to be the superior solution for Internet (i.e. lossy transmission with inconsistent latency) music and speech streaming applications. Specifically, it provides tight control over encoding latency when that's important. It's based on technology from two previous codecs: SILK, developed by Skype for speech transmission, and CELT, a codec developed for low-latency streaming audio.

One nice aspect of Opus is that it is an open standard, so a reference implementation is available that has source code for both floating-point and fixed-point computations. It also claims to have support for multiple architectures, including x86, ARM, and MIPS, which might be useful.

As far as quality goes, there are a number of detailed comparisons available between Opus and other audio codecs like AAC or MP3. Here is a plot that seems to summarize it well:

Basically, the makers of Opus claim superior audio quality to its competitors over a wide range of bitrates. Especially in the pretty-low-bitrate range of 16-32 kbps, it seems to have no peer at the present time. This page has a number of live tests that you can run in your browser, showing its performance with both speech and audio, its robustness to packet loss, and its ability to change bitrates (and therefore audio quality) on the fly.


Sounds like a job for dedicatedly identifying these surges and then actually subtracting them.

I'd start with a low-pass filter to find the slightly time-varying mean of the signal.

Use that to define lower and upper thresholds above or below you count something as surge.

Identify the samples lying outside the thresholds. Find a signal model for surges, e.g. a (truncated) sin x/x function or a quadratic function. Fit that model to the surge samples and then subtract from the signal.