What type of low-pass filter to use for gyroscope filtering in control application?

The lower the cut-off frequency of the anti-aliasing filter, the higher the delay thus degrading your phase margin. You need to be really careful when implementing an anti-aliasing filter in control loops applications.

First question, what is the amount of noise above 500 Hz? Based on the picture, you have noise at about 3 kHz, the main purpose of your anti-aliasing filter is to remove this spike, I would not worry too much about content between 500 Hz and 3 kHz.

Second question, what is the bandwidth of your control loop? 100 Hz? This is not the same as your sampling frequency. Typically it should be BW <= Fs/10. In your case BW <= 100 Hz.

Third question, what amplitude of aliasing can your application live with? If the bandwidth of your control loop is 100 Hz, aliasing occuring at f > 100 Hz will not really impact your performance, you only need to focus on aliasing that occuring at f <= 100 Hz. Even if there's aliasing occuring, if it less than say 1% of the signal of interest, you can probably neglect it.

Fourth question, can you run your control loop at 32 kHz? You would not need an anti-aliasing filter then. (I'm assuming that there's already an analog anti-aliasing filter in the gyroscope).

Based on the information we have, here's what I would do :

If you cannot sample faster than 1 khz, the biggest problem is that the 3-kHz spike will alias to DC (or near DC) . So you need to filter it. I would design a notch filter that removes your 3 kHz spike. This notch filter will add a smaller delay compared to a low-pass filter with a 500 Hz cut-off.

However, if you can sample faster, say 4 kHz, then the 3 khz spike will alias at 1 kHz which should be outside the bandwidth of your controller. In that case, I would only design a simple low-pass filter with a high cut-off frequency.

Edit : You mention you have an FPGA. In that case keep your sampling at 32 kHz and design a controller with a Bandwidth of like 100 Hz. In that case you don't need to worry about aliasing and noise since the noise will be outside the bandwidth of the control loop.


The Alamouti scheme guards against Rayleigh fading channels.  If you were to send two symbols within the same channel within the same coherence time interval, they will experience the same fading.  In order to get diversity gain, you need to send them in different channels (different antennas or different frequencies or distant time slots)


Why $H(A)=H(C)$ where $C$ is $A$ with an additional parity bit?

Because the entropy represents information quantity, or if being measured in bit, the smallest number of bits per symbol we need to represent a source.

The source $A$ contains $4$ equiprobable symbols hence it is obvious that we need $\log_2(4) = 2$ bits per symbol to represent the source.

The source $C$ is simply $A$ with its parity bits hence no new information is added. We have more bits per symbols, but the number of information bits is unchanged.

In communications, these additional bits are added to scope with error of transmissions.


Is convolution the only way to apply filters?

For this discussion it's important to restrict the class of filters to linear time-invariant (LTI) filters. Their input-output relation is described by the standard convolution sum (or, in continuous-time, convolution integral) that you've probably come across. So the operation of any LTI system can be described by a convolution. In discrete time you have

$$y[n]=\sum_{k=-\infty}^{\infty}h[k]x[n-k]\tag{1}$$

where $x[n]$ is the input signal, $h[n]$ is the system's impulse response, and $y[n]$ is the output signal. An example where $(1)$ is implemented directly is the transversal filter, which has a finite impulse response (FIR).

However, this doesn't mean that an LTI system can only be implemented by directly implementing the convolution sum ($1$). Any infinite impulse response (IIR) must be implemented recursively, because the convolution sum $(1)$ can't be computed exactly in finite time (since it's an infinite sum). But mathematically, any recursive implementation of an IIR filter computes the same thing as the convolution sum.

The convolution sum can also be implemented in the frequency domain, where blocks of data are transformed using the DFT, these frequency domain data are modified according to the desired impulse response, and after applying an IDFT, the data blocks are combined to give the desired output signal. Two such methods are the overlap add and the overlap save methods.

In sum, all LTI systems are described by a convolution sum (or integral), but there are many ways to implement that sum, and more often than not, the convolution sum is not evaluated directly.


Separate Complicated Signal into Exponential Components

Hilbert transform method from your earlier question and Jason R's answer should still work. What you are interested in is extracting the instantaneous phase of your function. For a signal $x(t)$ this is given by $\phi(t) = \angle x_a(t)$ where $x_a(t)$ is the analytic signal of $x(t)$.

Digitally you can do this using a Hilbert transform too. Here's a code snippet in Python:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
plt.style.use('presentation')

if __name__=='__main__':

  t = np.linspace(0,3,100)
  x = np.cos(10*t*t)

  y = np.imag(hilbert(x))

  x_a = x + 1j*y
  inst_phase = np.angle(x_a)

  true_phase = 10*t*t
  true_phase_wrapped = (true_phase-np.pi)%(2*np.pi)-np.pi

  plt.figure(1)
  plt.plot(t, x)
  plt.plot(t, inst_phase)
  plt.plot(t, true_phase_wrapped)
  plt.legend((r'$x(t)=\cos(10t^2)$', r'$\angle x_a(t)$',
      r'$10 t^2$ wrapped to $[-\pi,\pi]$'))

  plt.figure(2)
  plt.plot(t, np.unwrap(inst_phase))
  plt.plot(t, true_phase)
  plt.legend(('unwrapped inst. phase', 'true phase=$10t^2$'))
  plt.show(block=False)


Why does my stop-band filter attenuate other frequencies?

Filters always have an inherent "roll-off" in their frequency response, because you can't practically realize a pass-band that's a perfect rectangular function. For a low-pass filter, the point where the magnitude of your frequency response drops to -3dB is called the pass-band and everything beyond is called the stop-band (technically, everything beyond the corner frequency, but we'll take the corner frequency to be the -3dB level). How fast your frequency response attenuates beyond the pass-band is dependent on the filter length.

If you look at the frequency response of your filter, ht, you'll see that it drops to -6dB at 1000 Hz:

So it makes sense that the power drops by 6dB after filtering, which you see in your figure as a halving in amplitude.

If you had looked at the documentation for the fir1 function that you used, you would've realized this too (emphasis mine):

  B = fir1(N,Wn) designs an N'th order lowpass FIR digital filter
      and returns the filter coefficients in length N+1 vector B.
      The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0
      corresponding to half the sample rate.  The filter B is real and
      has linear phase.  The normalized gain of the filter at Wn is
      -6 dB.

Now to create sharper filters, with responses that get pretty close to a rectangle, you'll have to use IIR filters, which come with their own host of problems with stability, etc., but definitely an option. You can see my answer here for some ideas on implementing a discrete-form 2, second-order section IIR filter that gives very sharp corners. The example there is for a band-pass filter, but you can read up the docs for the functions used there and implement a low-pass version yourself.


Decomposition of a signal to slow and fast components

Your problem is pervasive in signal and image processing: being able to separate a trend (with some assumptions) from a signal (with some other assumptions) and potentially additional noise. This is an instance of (semi-)informed source separation: a sine from a drift,  a texture from a cartoon image, a biological information from instrumental artifacts, etc, that requires a more involved modeling and optimization framework.

Your choices are related to fundamental questions related to the model you can build based on physical and mathematical assumptions, aside from the algorithmic tools you can program.

1) Low-pass filtering:
Standard low-pass filtering is linear and time-invariant. Linearity is very consistent with "signal = tonic + phasic". But phasic could be (on theory) negative, as "$2+(-1) = 1$". As your graphic does not include x/y scales, "above" is clear, but i am not sure of what you can model with "no points in the recorded signal should be below the tonic component", since additional noise can be fluctuating too. But positivity is a property which can be taken into account, with some non-linearity caveats.

2) Duration:
"that are at least as far as the duration of phasic response" : your assumption seems to me to related to a time-delay, yet i cannot interpret it precisely so far.

3) High-order derivatives:
"A third approach would be to work on first or second derivative space and characterize the system in this domain, as it is supposed to be less sensitive to low-variations" Indeed, slow variations are expected to be with limited derivatives. This could result in a "sparsity" prior assumption.

4) Causality: "rely more on the history" there is some causality assumption behind that. This too is a tough stuff.

All your four approaches seem consistent to me, and may help you ponder and state your assumptions. They are otherwise termed as “morphological component analysis (MCA)”, “geometric separation” or “clustered sparsity”, adding positivity. Such topics are addressed for instance in "Blind source separation of positive and partially correlated data" http://dx.doi.org/10.1016/j.sigpro.2005.03.006.

You may find additional history for instance in the introduction of (mong other papers "Chromatogram baseline estimation and denoising using sparsity (BEADS)"
http://dx.doi.org/10.1016/j.chemolab.2014.09.014
which also deals with positivity, low-pass-filtering and derivative sparisty, yet the causal approach is not considered. So you can specify your model in more precise terms (you did do talk about the noise yet), you may find some tracks to solve your questions


So far what you are proposing sounds like a reasonable approach. However, I don't think you will know how well it works until you try it, just like you have tried SIFT.

I have a question though.  Why are you restricting yourself to DCT?  There are lots of representations that have been used for texture classification: co-occurrence matrices, local binary patterns, etc.  The fact that you have only found one paper on using DCT for texture classification would suggest that this is not the most commonly used feature for this problem.  I would recommend that you widen your literature search to see what other features people have used, and how well they have worked.


The MAP criterion minimizes the average probability of error, which in your case reads
\begin{equation}
\begin{aligned}
P_e &= P \left( -1 \right) P \left( Y>\theta \mid -1 \right) +
P \left( +1 \right) P \left( Y<\theta \mid +1 \right) \\
&= 0.2 \int_{\theta}^{+1} f_Y \left( y \mid -1 \right) \mathit{dy} +
0.8 \int_{-1}^{\theta} f_Y \left( y \mid +1 \right) \mathit{dy} \,.
\end{aligned}
\end{equation}
Therefore, the optimum threshold, call it $\hat{\theta}$, is optimum in terms of minimizing $P_e$.
Exploiting the previous equation and the definitions of the random variables $X$ and $N$, you get
\begin{equation}
\begin{split}
P_e &= 0.2 \times 0.25 \times \left(1-\theta\right) +
0.8 \times 0.25 \times \left(1+\theta\right) \\
&= \frac{1}{20} \times \left(1-\theta\right) +
\frac{4}{20} \times\left(1+\theta\right) .
\end{split}
\end{equation}
Minimizing $P_e$ w.r.t. $\theta$, you obtain $\hat{\theta} = -1$.
An easier way to solve your problem would be: look at the plot, and select $\hat{\theta}$ such that it minimizes the sum between the right tail of the left function (i.e., the grey area) and the left tail of the right function. In this way you easily get the solution, but the analytical way is, in my opinion, preferable to understand the concepts.


Gaussian pyramid: why needs the image to be downsampled

In the context of a gaussian pyramid, why is the image downsampled separately although the numbers of pixels are decreased through smoothing already?

The Gaussian filter adjusts the bandwidth of the content of the image. The downsampling adjusts the spatial resolution of the image.

The Gaussian filter is a low pass filter. After its application to an image, only the frequencies within its pass band will be present in the resulting image.

An image of dimensions $M \times N$ pixels can depict $\frac{M}{2} \times \frac{N}{2}$ lines per pixel of spatial frequencies.

If you apply a low pass filter at $\frac{M}{4}, \frac{N}{4}$, you have reduced the spatial frequency requirements for that signal by half. Where before you needed $\frac{M}{2} \times \frac{N}{2}$ lines per pixel, now you know that spatial frequencies in your image do not exceed $\frac{M}{4}, \frac{N}{4}$ lines per pixel. The rest of the bandwidth that an image of that resolution can support goes unused.

At that point, you can reduce the image resolution to accommodate the "new" limits of the bandwidth without any loss.

From another point of view, you use the low pass filter to ensure that the resulting image will not contain any spatial frequency that could be aliased if it was to be plainly subsampled (because ultimately, resizing an image means that it will have to be resampled and "fit" into a smaller space. This smaller space means less bandwidth). Once you ensure that no frequency will be aliased, then it is safe to resize.

  An image with 8 x 5 pixels smoothed with a gaussian filter of 3 x 2 pixels will create a 6 x 4 pixel image if the filter can only be slided within the boundaries of the original image without crossing the edges.

This is to do with choosing only the "valid" values in the output of the filter. Values for which all the weights of the two dimensional impulse response are applied to image pixels.

Hope this helps.


How to reverse frequency spectrum of a 2D image?

I think that you should see it as doing it in one dimension and then the other, meaning you'd invert the values of every other pixel in every row, then once it's done do it again for every column, which I suppose means you can do it directly in 2D by inverting pixels in a checkered pattern like this:

1 -1
-1 1

As for your FFT/IFFT question yes it should be lossless, save for small computation errors, and I'm not sure but it might only work with even sizes so the Nyquist and DC components can switch places. However in both cases you'll have to deal with negative pixel values, so that kind of limits you to your own algorithms, unless you do $f(v) = 0.5+0.5v$ on all your values and then undo it.


Exponential average with time constant of slow, fast and impulse

Exponential averaging is performed using a first order IIR filter according to the difference equation:
\begin{equation}
y[n] = (1-\alpha)\;x[n]+\alpha \; y[n-1]
\end{equation}
where $x[n]$ and $y[n]$ are the input and output signals at time $n$ and the value of the single coefficient $\alpha$ is given by the equation:
\begin{equation}
\alpha = \exp( \frac{-1}{f_s \tau} )
\end{equation}
where $f_s$ is the sampling frequency and $\tau$ is the time constant.

In this light, it can be seen that exponential averaging filters are simply low-pass filters, the roll-off of which is -6dB/oct.  The only free parameter is the cut-off frequency, which is specified by the time constant.  Longer time constants will yield a lower cut-off frequency, as can be see in the Matlab code and plot below.

figure; hold all
fs = 8000;
tau = 10.^(-1:-1:-4);
alpha = exp(-1/fs./tau);
color = 'brgc';
for index = 1 : length( alpha )
    [Hz fVec] = freqz( [1-alpha(index)], [1 -alpha(index)], 8192, fs );
    plot( fVec, 20*log10(abs(Hz)), color(index), 'LineWidth', 2 );
    legText{index} = sprintf( '$\\tau = 10^{%d}$', log10(tau(index)) );
end
hLeg = legend(legText,3);
set( hLeg, 'interpreter', 'latex' );
axis( [8000/8192 4000 -70 10] )
ylabel('Amplitude [dB]')
xlabel('Frequency [Hz]')
set( gca, 'xscale', 'log' )

Notice the parallel -6dB/oct slopes (characteristic of single pole filters) and the falling cut-off frequency with lengthening time constants.

The correct choice of time constant depends on how fast the variation of the desired average is.  In order to capture a quickly varying mean, only signal content above the fastest frequency of variation should be attenuated -- the time constant should be chosen accordingly.  This means that if the signal mean varies very slowly, long time constants are needed, and vice versa.  I find that trial and error works best to determine the optimal value, but formulas to calculate the cut-off frequency do exist, e.g. see this link.

Good luck!


Is there a mathematical method to determine if noise is Gaussian?

There are several statistical tests if a time series is Gaussian, although in statistics, the term "tests for normality" is usually how you search for them.

The Nist EDA site is a good place to look and the probability plot is better for shorter data sets than the sample histogram.

http://www.itl.nist.gov/div898/handbook/eda/section3/probplot.htm

Near the bottom of the page, there are references to q-q plots, KS, Chi squared, and other goodness of fit tests.  You can find ample information about them on the web and replicating here isn't going to add anything.

Matlab has qqplot and prob plot in the Statistics toolbox, and the qqplot with a single argument is specific to Gaussian distributions. SAS has all these tests. R has the tests.

I recommend this book, written by 2 Engineers, and they cover several tests including for things like independence, and stationarity.  The book is oriented towards the practical, minimum of mathematics.

  Bendat, Julius S., and Allan G. Piersol. Random data: analysis and
  measurement procedures. Vol. 729. John Wiley & Sons, 2011.

The wrinkle of these tests is that they don't conform to a Signal plus Noise scenario.  The tests generally assume that the time series is all Gaussian or not. A constant mean isn't a problem. Signals are not usually Gaussian and a simple test can't tell the difference.

Signal processing operations such as a DFT, tend to manifest central limit theorem effects on data, so you need to be aware that even linear transformations will not preserve a non Gaussian pdf.

It should be also noted that from a practical perspective, Gaussianity isn't black and white.  Algorithms that have Gaussian assumptions usually work well even if the Gaussianity assumption is not strictly valid.  Things like bi-modality and non-symmetry are more important to know about. Cauchy (heavy tails) like noise and multiplicative noise are also important to know about.


Explaining / Communicating spectra: Generating Visualizations

Personally, I value aesthetics and consistency in my documents, so I use Latex whenever I can. I use Inkscape if I'm in a hurry, though.

I have found that the pgfplots package (built on top of tikz) reduces the time needed to code a plot substantially, especially once you get the hang of it. As a simple example, the following filter response:

was generated with the following Latex code:

\documentclass{standalone}
\usepackage{pgfplots}
\begin{document}
\begin{tikzpicture}
\begin{axis} [
        height=5cm, width=7cm, xlabel={$f$}, ylabel={Filter A},
        xmin=0, xmax=300, ymin=0, ymax=1.6, no markers,
        xtick={100,175,250}, ytick={0.25,0.5,0.75,1},
        axis lines=middle
    ]
    \addplot coordinates {(0,1) (100,1) (250,0)};
    \addplot[black,dotted] coordinates {(0,0.5) (175,0.5) (175,0)};
\end{axis}
\end{tikzpicture}
\end{document}

Excluding boilerplate (which can be reused between plots), this is just two lines of "real" code (the \addplot lines). The benefits are the usual latex quality, plus font/style/size consistency with the rest of the document.

As a further example I tried a (simplified) version of the last plot in the question:

which requires just a bit more code:

\documentclass{standalone}
\usepackage{pgfplots}
\usepgfplotslibrary{fillbetween}
\begin{document}
\begin{tikzpicture}
\begin{axis} [
        height=5cm,width=9cm,xlabel={$f$},axis lines=middle,
        xmin=-120,xmax=200,ymin=0,ymax=1.8,no markers,
        xtick={80,160},xticklabels={$\frac{f_s}{2}$,$f_s=\frac{1}{T}$},
        ytick={0},yticklabels={},ylabel={$S_s(f)$}
 ]
    \addplot[black] coordinates {(-100,0) (-20,1) (0,0)};
    \addplot[black] coordinates {(0,0) (20,1) (100,0)};
    \addplot[black,name path=A] coordinates {(60,0) (140,1) (160,0)};
    \addplot[black,dashed,name path=B] coordinates {(80,0) (80,1.5)};
    \addplot[blue] fill between[of=A and B,soft clip={domain=60:80}];
\end{axis}
\end{tikzpicture}
\end{document}

It is true that generating these plots takes time. When working on anything but the most trivial plot, I constantly need to refer to the manuals and to tex.stackexchange.com. I think the examples above show that the pgfplots package removes some of the pain; once you get practice with the most common commands, it shouldn't take too long to create nice plots.

A final point: while I don't (yet) have a program to automatically generate spectrum plots, it shouldn't be too hard to do. I do have a script to generate pgfplots scripts from Matlab/Octave, which you can see here. This script takes two or more Matlab arrays and a plot description, writes a latex file, and runs pdflatex on it. While laborious, coding it wasn't particularly difficult. It should be possible to modify it to take spectrum descriptions (maybe in the edges, amplitudes format of the filter design commands) and generate the plots automatically.


If $x(t)$ and $X(f)$ are Fourier transform pairs then the value

$$\int_{-\infty}^{\infty}|X(f)|df$$

can serve as an upper bound on the magnitude of $x(t)$:

$$|x(t)|=\left|\int_{-\infty}^{\infty}X(f)e^{j2\pi f t}df\right|\le \int_{-\infty}^{\infty}|X(f)|df$$


How to calculate the frequency response of LTI system based on measured impedance

this is more about electronics, but i think it can be asked here.

when you are looking for a frequency response of an impedance, i believe this means you have a two-terminal (a.k.a. one-port) device.  so, instead of an input voltage and output voltage to relate by a transfer function (and from the transfer function you get a frequency response), what you have is a voltage applied to the two terminals and a current measured, or a current shoved into (and out of) the terminals and the resulting voltage measured.

the complex impedance of the two-terminal device, as a function of frequency, is the frequency response the device assuming current is the input and voltage is the output.  you can plot it on a dB scale, but without a reference resistance, you don't have a reference for 0 dB on that frequency response plot.  but you have relative dB for the different impedances at different frequencies.

maybe ask the EE guys this question.  i dunno.


VOFDM (Vector OFDM) advantages and disadvantages

Referring to VOFDM as described in this reference, some thoughts on that (I wrote this summary some years ago, due to an issue in my own research). We can understand VOFDM as quite similar to interleaved SC-FDMA, (this means that a single carrier occupies several frequency bins, which are not grouped together but spread around the full bandwidth of the system). So, assuming the system has $N=KM$ OFDM subcarriers (FFT window length N, N frequency bins), the frequency indices occupied by an interleaved SC-FDMA carrier would be $[l, l+K, l+2K, l+(M-1)K]$, where $M$ is the number of frequency bins per SCFDMA carrier and $K$ is the number of different carriers possible in the full bandwidth. $l=0,\dots,K-1$ is the index of the SC-FDMA carrier.  Then, after FFT at the receiver, the received signal for the $l$ carrier is given by

$$Y_l=H_ld_l+n_l$$ where $Y_l$ is the received signal at the frequency bins for the $l$ carrier, $d_l$ is the transmit data for the $l$ carrier and $n_l$ is AWGN. $H_l$ is the equivalent channel matrix, which can be considered a circulant matrix.

For the considerations below, VOFDM is actually doing exactly the same thing (compare eq (4) of this reference). So, how can we compare OFDM and VOFDM?

The big advantage of OFDM is low complexity at the receiver. A single DFT plus one-tap equalization is sufficient for the detection. In contrast, VOFDM requires additional M-point DFTs on each carrier to perform one-tap equalization.
In OFDM, each frequency bin can be treated independent from the others (it's orthogonal). This allows great implementation speedups in terms of parallel processing. VOFDM can be only parallelized to $K$ "threads".
The complexity aspect comes even more into play, when we consider MIMO transmission (e.g. spatial multiplexing). With today's hardware we actually have the chance to use a sphere-decoder in real-time on each carrier to get the optimal receiver performance. For VOFDM, the dimension of the sphere-decoding would be M-times larger, which quickly becomes infeasible for today's hardware.
VOFDM can exploit frequency diversity, since one data symbol is spread around multiple, far apart, frequency bins. However, in order to harvest this diversity, you'd at least need an LMMSE equalizer. However, even better performance can be obtained by employing a sphere-decoder on the equivalent channel $Y_l=H_ld_l+n_l$. But, this requires very high implementation complexity.  Also, it needs to be investigated, if additional (smart) soft-output decoding of the (certainly) encoded input data stream can gain the same/similar diversity in the OFDM case.
In principle VOFDM can achieve better FER performance when considering iterative receivers due to its frequency diversity and concatenation of a channel code and the constellation constraint of the equivalent channel. However, again, to achieve this performance high calculation complexity needs to be accepted.
regarding robustness against synchronization errors and channel estimation errors, I'd assume both systems are similar. However, Channel estimation itself can be more challenging for VOFDM, since no orthogonal frequency bins exist which can be nicely used for dedicated pilot frequency bins.

So, to summarize: VOFDM is a nice idea and has potential to perform better in terms of FER than OFDM. But, it needs a lot more complexity.


Why is the difference of gaussians scale space scale invariant?

The term "scale-invariant" means the following here.
Let's say you have image I, and you have detected a feature (aka an interest point) f at some location (x,y) and at some scale level s.  Now let's say you have an image I', which is a scaled version of I (downsampled, for instance).  Then, if your feature detector is scale-invariant, you should be able to detect the corresponding feature f' in I' at the corresponding location (x',y') and corresponding scale s', where (x, y, s) and (x', y', s') are related by the appropriate scaling transformation.

In other words, if your scale-invariant detector has detected a feature point corresponding to someone's face, and then you zoom in or out with your camera on the same scene, you should still detect a feature point on that face.

Of course, you would also want a "feature descriptor" which would allow you to match the two features, which is exactly what SIFT gives you.

So, at the risk of confusing you further, there are two things that are scale-invariant here.  One is the DoG interest point detector, which is scale-invariant, because it detects a particular type of image features (blobs) irrespective of their scale.  In other words, the DoG detector detects blobs of any size.  The other scale-invariant thing is the feature descriptor, which is a histogram of gradient orientation, which stays more or less similar for the same image feature despite a change in scale.

By the way, the difference of Gaussians is used here as an approximation to the Laplacian-of-Gaussians filter.


The Meaning of the Terms Isotropic and Anisotropic in the Total Variation Framework

In the Total Variation framework we define 2 flavors:

$$ \text{Isotropic TV} \; {TV}_{ {L}_{2} } \left( X \right) = \sum_{ij} \sqrt{ { \left( {D}_{h} X \right) }_{ij}^{2} + { \left( {D}_{v} X \right) }_{ij}^{2} } $$

$$ \text{Anisotropic TV} \; {TV}_{ {L}_{1} } \left( X \right) = \sum_{ij} \sqrt{ { \left( {D}_{h} X \right) }_{ij}^{2} } + \sqrt{{ \left( {D}_{v} X \right) }_{ij}^{2} } $$

Where $ {D}_{h} $ and $ {D}_{v} $ are the horizontal / vertical finite differences operators.

Quoting Isotropic and Anisotropic Total Variation Regularization in Electrical Impedance Tomography:

  It is worth noticing that referring to functional as ”anisotropic” is slightly misleading, because the weights corresponding to all coordinate directions are equal.

The term Anisotropic is used because this flavor separate Vertical and Horizontal components as done in Anisotropic Diffusion. Hence it doesn't have rotational invariant term as the Isotropic Flavor. As you know, sum of quadratic terms is connected to Ellipse and with equal weights it is a Sphere / Circle which is Isotropic.

Great reading about the subject is Laurent Condat - Discrete Total Variation - New Definition and Minimization (DOI Address). It analyze the properties of different TV definitions regarding their spatial smoothing.


Is there any difference in Implementing sigma delta modulator using filters and state space model in FPGA?

Probably not much advantage either way, possibly unless you've got a really high-order modulator (in which case someone else will have to answer!)

If the state-space representation is a single-input single-output system, and if it's linear, and if it's running in steady-state, and if data path widths aren't an issue, then there's no difference between it and a "more traditional" IIR filter implementation such as a set of biquads.

If you've got some super-high-order modulator, or if someone's gone and invented a nifty nonlinear way to get more performance out of a sigma-delta modulator when I wasn't looking, then maybe -- but I doubt it.


Frequency response and sampling theorem for triangular function

You're right that a triangular function is not band-limited and that it can't be sampled without introducing aliasing. However, this doesn't mean that linear interpolation can't be useful. The usefulness of linear interpolation depends on the data to be interpolated and on the desired interpolation factor. If the data have a low pass character, i.e., if they are already sufficiently oversampled, then linear interpolation can be good enough if the interpolation factor is not too large.
The obvious advantage of (bi-)linear interpolation is its simplicity and computational efficiency. In image processing it performs better than the common nearest-neighbor interpolation. However, if there are no real-time requirements and if computational efficiency is no big issue, then the standard interpolation method in image processing is the more complex bicubic interpolation.


Finding an Image's Orientation With DFT Frequency Amplitudes

Have you heard of the Hough Transform? It is a possible alternative for tackling your problem.


What is sound output delay of a band pass filter at different freqencies?

The lower the frequency, the longer the period of oscillation.  In any filter using memory (IIR or causal FIR), the filter will have to evaluate a signal for a longer period of time in order to gather information about the same number of periods of a lower frequency signal component than for a higher frequency signal component.  And a filter can't filter a particular spectral component unless it has gathered enough information about that spectral component to "know" whether that particular component is in the stop-band or pass-band, which takes longer (other factors being roughly equivalent) for a lower frequency signal component.

Or you could design a linear phase FIR with the same delay for a lower frequency bandpass, but since it has gathered less periods of information about that lower frequency, it's quality of filtering (width of the transition band, etc.) will be less.


How to choose the "best" measurment (from a given set) as input for a kalman filter?

Question: Which parameter is suitable to indicate how "good" the
  measurement fits to the Kalman filter?

To estimate a quality of association you can use likelihood function. The likelihood considers not only residual but also uncertainty and represented as scalar value:

$$\mathcal{L} = \frac{1}{\sqrt{2\pi S}}\exp [-\frac{1}{2}\mathbf y^\mathsf T\mathbf S^{-1}\mathbf y]$$

  Problem: For an object tracking scenario with multiple objects and
  multiple tracks, I want to choose the "best" assignment
  object<->track.

The question which you are interested is called data association and as was said there are a lot of solution and this question is too broad. Here are some of methods (you can google it and find more information):
 - Global nearest neighbour (in Euclidean space);
 - Global strongest neighbour;
 - Joint probability data association;
 - Multiple hypothesis tracking (the best choice).


Expressing DTFTs in terms of one another with similar time domain signals

The problem is asking you to find an analytical expression given $x_1[n]$ and $x_2[n]$. There is no need to compute anything, let alone code it up!
The problem wants you to identify that given the discrete-time domain functions, you use properties of the DTFT to write $X_2(e^{j\omega})$ in terms of $X_1(e^{j\omega})$. You've already identified that $x_2[n]$ contains a shifted version of $x_1[n]$.
You can use the time-shift property which says that if the DTFT of $x[n]$ is $X(e^{j\omega})$ then
$$\mathcal{F}(x[n - m]) = e^{-j\omega m}X(e^{j\omega})$$
Keeping in mind how the domain of $x_2[n]$ is defined, you can use the DTFT properties to express $X_2(e^{j\omega})$ in terms of $X_1(e^{j\omega})$.


Does Wavelet Packets (WPT) decompose signal on dyadic scales like DWT?

I suppose  you are talking about $2$-channel DWT. WPT may have  slightly different  interpretations. With DWT, one only allows to decompose frequencies "downward ", ie after lowpass and subsampling. WPT allows to further decompose, at each level, after (lowpass or highpass) and subsampling. In that sense, WPT are a generic scheme, of which the DWT denotes a specific branch. The decision to decompose or stop at each level is governed by rules "external" to wavelets: frequency target, entropy, etc.

Another vision is the uniform WPT, or full wavelet packet, where at each level, both outputs are further decomposed, generating a uniform frequency splitting.

So basically, you are correct at question 1. For question 2, when a DWT subband is approximately located (in reduced frequencies) between two successive dyadic bounds $[2^{-j-1},2^{-j}]$, a WPT can be located between any two dyadic rational bounds $[k_1 2^{-j_1},k_2 2^{-j_2}]$, with $j, j_1,j_2$ nonnegative integers, $k_1 \le 2^{j_1}$ and $k_2 \le 2^{j_2}$.

For instance, the upper DWT subband only (approximately) localizes stuff in $[1/2,1]$, while you can reach  $[5/8,3/4]$ with a WPT. However, the aliasing comes into play, and this $[5/8,3/4]$ is disturbed by aliased components from other frequencies.  Aliasing is considered more severe for high frequency wavelet packet coefficients.

With $M$-band DWT and WPT, both structures remains $M$-adic at each level, but the WPT tends to split the subband axis onto $M$-adic fractions, with form $k/M^j$.


What is difference between OFDM and GFDM?

There is a great interactive demonstration of GFDM by Xiaojie Wang and Stephan ten Brink explaining the differences between OFDM and GFDM.

To sum it up, GFDM uses generalized pulse shaping filters that can span

a subset of $N_\mathrm S$ subcarriers out of a total number $N_\mathrm{S,tot}$ of subcarriers and
a set of $N_\mathrm Y$ symbols

OFDM is a special case of GFDM where $N_\mathrm Y = 1$, $N_S = N_\mathrm{S,tot}$, and the pulse shaping filter is rectangular in frequency domain. This can be implemented using an IFFT at the transmitter and an FFT at the receiver.

The implementation of GFDM is generally not that simple. I could think of an implementation using $N_\mathrm{S,tot}/N_\mathrm S$ IFFTs of length $N_\mathrm{S,tot}$ where at every IFFT just a subset of $N_\mathrm S$ inputs is used and the others are zero padded to realize the frequency allocation of one subset of subcarriers. Each IFFT output would be followd by a specific pulse shaping filter spanning $N_\mathrm Y$ symbols. Finally, all filter outputs would be added up and converted to an analog signal using an digital-to-analog converter (DAC).

The purpose of GFDM is to achieve, amongst others:

lower out-of-band radiation (better spectral efficiency)
higher robustness against carrier frequency offset (CFO)
higher robustness against sampling time offset

when compared to OFDM.

GFDM is also motivated by the fact that in Long-term evolution (LTE) mobile communications, users are allocated so-called physical ressource blocks (PRB). PRBs consist of a set of adjacent subcarriers and symbols.


The answer is posted as a reply to Otsu's method problem? . So it is now a duplicate I guess.


Non-Linear, Non-Stationary spectral analysis methods! When and where?

Non-linearity and non-stationarity are non-properties. Without more details, they do not say much about the methods that may perform well, and moreover the  choice depends a lot on what you really do: analysis, feature extraction, enhancement, filtering, component separation, restoration? What follows are typical sets of tools you could use:

Your moving-window-Fourier definitely has some STFT flavors, and there are many non-linear extensions: bilinear or quadratic time–frequency distributions (Wigner-Ville, Choi–Williams, Rihaczek)
Non-linearity may come from changing the classical energy preservation into other norms, like robust DFTs (under an $\ell_1$ norm minimization), see for instance Robust $L$-estimation based forms of signal transforms and time-frequency representations, see also Using different algorithms/formulas to create spectrograms
Some methods also alleviate the hypothesis that signal formation is linear: variance-stabilizing transforms,  cepstral analysis or homomorphic signal processing (linear transforms applied in a $\log$ domain)
Generalizations of linear approaches in polynomial or rational fraction modeling, with a recent interest in Rational Optimization for Nonlinear Reconstruction with Approximate $\ell_0$ Penalization.

All of the above can be combined somehow, like EMD at different scales in Intrinsic multi-scale analysis: a multi-variate empirical mode decomposition framework. EMD has many uses and some known limitations at the same time.


how to understand frequncy response and lowpass filter?

$H(j\omega)$ comes from the more general response $H(s)$ which comes from the Laplace transform using the complex unit s which has real and imaginary components given as $$s= \sigma+j\omega$$

This may be too much to swallow if you haven't been exposed to the Laplace Transform yet, but let me make an attempt to simplify it: The Laplace transform decomposes analytic single valued functions into a continuous function of s by correlating the function with all values of $e^{st}$  If you are asking about j you are probably equally confused about $e^{st}$. Note that $Ae^{j\theta}$ is simply a function with magnitude A and angle $\theta$, so if s is only imaginary ($s=j\omega$), then $e^{st}=e^{j\omega}$ and the result is a function whos phase is changing with time. I like to view it as a phasor (vector) of constant magnitude rotating around the complex plane, thus we see that the real part of this is a cosine wave and the imaginary part is a sine wave resulting in Euler's Identity $e^{j\omega}= cos(\omega)+jsin(\omega)$.  When s has real components ($\sigma$) then the time domain function also decays or grows depending on the sign of $\sigma$. See the plot below demonstrating this:

Now to relate this to a low pass filter see this next plot below. This shows the surface plot of H(s) for a low pass filter with two poles in the left half plane:

The frequency response is H(s) when we restrict s to just be imaginary such that $s=e^{j\omega}$, therefore we are seeing the result when we only apply fixed magnitude sine and cosine waves to the input of the filter, omitting all the cases where the input signal may decay or grow- which makes sense when we want to measure the frequency response: apply a fixed frequency at the input and measure how much this one frequency has been attenuated at the output (magnitude response), and how much the signal has shifted in phase (phase response). And repeat this for every frequency (which just means sweep s along the $j\omega$ axis from $-\infty$ to $+\infty$.  This is in fact the Fourier Transform as indicated in the plot above showing an expected response for a low pass filter (in this case one that has peaking near the resonance points based on the pole locations.

Here is another example very close to yours but showing the same relationship of H(s) and $H(j\omega)$

And finally the Bode Plot showing the magnitude and phase versus frequency for this low pass filter, which is specifically the magnitude and phase of H(s) when $s=j\omega$ or otherwise stated more directly as $H(j\omega)$. Electrically what is happening is the capacitor becomes a smaller and smaller impedance as the frequency increases ($X_c = \frac{1}{2\pi fC}$). This capacitor is in parallel with what is typically a resistive load. At low frequencies, DC for example, the capacitor is an infinite impedance or open circuit, as if it is not even there and the voltage at the load will simply be $V_{in}\frac{R_L}{R+R_L}$ but as the frequency increases the output will be $V_{}in\frac{R_L||X_C}{R+R_L||X_C}$ where $||$ represents the parallel combination of the two impedances.  What we see from this is that the voltage at the load will decrease as the frequency increases (low pass filter). At high frequencies (RF and microwave) we deal more with reflection coefficients in that the high frequencies are actually reflected back to the source, but I won't get into that here.


You should parameterize the path as a parameter of time.
You can do that off line with accurate measurements of the path.
Then use Non Linear Least Squares to find the best match between the reads of the IMU to the model.

From there you'd be able to match the speed.


Decimating by last interpolating factor? (Multi-stage)

The idea is that you'll be saving computational steps in the filtering part of the decimation. Think about it, if you have 800 samples then downsample by a factor of 8, now you only have to low pass filter 100 samples for the first stage. If you instead tried to downsample by a factor of 2 first, then you'd have 400 samples to filter after the first stage. The order matters because of the number of computations taken in the filtering step.


How to reconstruct original signal from sampled signal?

To implement a digital-to-analog converter, all you need is an ideal low-pass filter to filter out the periodic frequency response that $\varOmega \geq \varOmega_s/2$.
$$
H(j\varOmega) =\left\{
\begin{aligned}
T,\ \ \ \  |\varOmega|<\varOmega_s/2 \\
0,\ \ \ \  |\varOmega|\geq\varOmega_s/2
\end{aligned}
\right.
$$
The impulse response of an ideal low-pass filter is a sinc function
$$
h(t) = \frac{\sin(\pi t/T)}{\pi t/T}
$$
Convolving the sampled signal with this impulse response gives the analog signal.
$$
\begin{aligned}
y_a(t) =& \sum_{n=-\infty}^{\infty} x_a(nt) h(t-nT)\\
=&\sum_{n=-\infty}^{\infty} x(n) \frac{\sin[\pi (t-nT)/T]}{\pi (t-nT)/T}
\end{aligned}
$$
f1=2;
f2=5;

fs=100;
Ts=1/fs;
t=0:Ts:1;
xt=cos(2*pi*f1*t)+cos(2*pi*f2*t);
figure
plot(t,xt)

fs=12;
Tn=1/fs;
tn=0:Tn:1;
xn=cos(2*pi*f1*tn)+cos(2*pi*f2*tn);
hold on
stem(tn, xn)

y = dac(xn, Tn, t);
plot(t,y)

function y = dac(x, T, t)
m = 0:length(x)-1;
y = zeros(1, length(t));
for ii = 1:length(t)
    h = sinc((t(ii)-m*T)/T);
    y(ii) = sum(x .* h);
end
end


A use of such window functions would be to window data that "goes straight to business" starting at the first sample and then gradually dies off, like an impulse response of some system. See What is the meaning of half window functions? A smooth rise in the window function would be detrimental as it would attenuate the important peak, unless the rise in the window is quick enough and the peak delayed enough. However, I think it would be a stretch to call something like that a half-window, if it is not a half of some window. So start directly at 1 at the peak of a symmetrical window I think.

Then again, beware that a sine window a.k.a. cosine window may also be called half-sine window or half-cosine window, as it represents half a cycle of a sinusoid (sine or cosine).


Virtual Microphone using Mic Impulse Response

Will this work? Well, it will do something, but it will not magically make one microphone sound like another. A microphone is much more than a signal receiver with an impulse response. The character of a microphone is determined by its directional characteristic and how it couples to diffuse and direct field components, the kind of coupling to the pressure field, the nonlinearity of the membrane and integrated preamplifier, the noise of the membrane and the electronic components.

You will only be able to capture a small part of all that with an impulse response, and very likely exactly not the part of a microphone that is responsible for the praise of its sound.

But even if you ignore all that, you're still left with the problem of deconvolution that is not trivial, especially in the presence of noise of a potentially unknown character. Your spectral division approach is rather naive and will work to a certain degree. If you want to get good results, you will have to take a harder route and you will probably not even in that aspect get what you imagine. Also, you will always have to start with a very good microphone. A bad microphone does not allow for such a simple improvement usually.

In short, even if I sound overly pessimistic: No, this will not work how you probably expect it to. If you want to do this to start a research project then you're good. But it's not going to be a plug and play product for microphone aficionados.


Mathematically or information-wise, you're doing the same thing.

However, filtering in time domain means convolving the signal with a filter kernel, so the algorithmic complexity would be something like $O(n^2)$ in terms of float additions with $n$ the length of the signal.

In frequency domain, filtering is an element-wise multiplication which algorithmic complexity is at most $O(n)$. You also have to consider the complexity of switching between time-domain and frequency-domain, meaning using an FFT algorithm. The FFT complexity is most of the time $O(n\log(n))$.

So the overall complexity of frequency domain filtering would then be $O(n\log(n))$, which is way more efficient than $O(n^2)$.

Considering for example a 1024-long signal, time-domain filtering would need more than 1 millions operations, while frequency-domain filtering would need approximately 3000 !


... Instrument Recognition using pyAudioAnalysis. I was wondering if results could be improved by also analyzing the ADSR shape of the samples.

In a pattern recognition setting, additional features could lead to better classification results, provided that they can better "differentiate" the phenomenon being "measured".

Because of psychoacoustics, human beings are capable of distinguishing timbre from the first few milliseconds. As an extreme example, if you listen to a guitar and a similar piano tone, half way through they have started playing, it is very difficult to tell them appart. For example, here is Yngwie Malmsteen sounding like a flute, a violin, a guitar (or two) just by modulating the Attack via the guitar's volume (notice his right hand after the "without echo" part).

Now, for some piano tones and some guitar tones it might be difficult to tell them appart just by their spectral content because they are in similar ranges. But if you hit some bass tones on the piano, really hard, they will Sustain differently because the strings are now less of a string and more of a rod, which vibrates differently (than a taut string). In that case of course, the timbre will start having differences too.

So, ADSR based features would help in instrument recognition when the spectral characteristics of two or more instruments are near enough to be mixed up by a classifier using just spectral features. (For more information, please see this link).

How you are going to do this is another issue because it's not as easy as just running some metrics on some recordings. The instruments would have to be recorded "clean" in a room with minimal reflections and under controlled conditions. Certain instruments will sound differently for the same tone if it was to be played softly or loudly and the range of each instrument would have to be taken into account to decide on the Sampling Frequency. All of theses parameters will have an impact to your classification performance.

So, assuming a set of clean recordings, you could obtain the envelope of the waveform  to then try and infer the ADSR phases from that curve and use these figures in your classifier. For a brief overview of automatic ways of extracting the envelope please see this link while an example (using the analytic method of deriving the envelope) in Python is available from this link.

Hope this helps.


What is the advantage of using Bessel function as the basis in a Hankel Transform?

The Hankel Transform includes the Bessel function as part of its definition. I wouldn’t say the Bessel has an advantage. I would say that the Hankel Transform has an advantage where a two dimensional signal has radial symmetry, or saying it another way, the signal in $x$ and $y$ can be expressed as a product of two functions in $r$ and  $\theta$ through a change of variables.
$x=r\cos \theta$ $y=r\sin \theta$

$$f(x,y)dx dy \Rightarrow g(r) h( \theta) r dr d \theta $$

The Hankel Transform is equivalent to a Fourier Transform on the radius $r$ variable.

Any advantage over STFT or Wigner would need to be explained in the literature you read.  Im not aware of any general use of Hankel grams, or why the Hankel would make any obvious relevance outside of two dimensional geometries with rotational symmetry.


Basicly DEMON is looking for peridicidies in the signal envelope, like the swish ,swish, swish, ....sound that a propeller of a big ship makes.  The movie Das Boot, has some nice examples.  It's a good way to extract balde rate.  Modern skewed multiblade propellers don't produce much for DEMON to see. I would have expected to see something in DTIC which I interpret more as this is so old, why bother.

You need to make some choices like the band you expect to see a modulated envelope. It needs to be in the right place and wide enough, but not too wide.  While not obvious in any papers, running more than one DEMON is permitted.

You can thank Australia and find a block diagram in:

http://digext6.defence.gov.au/dspace/bitstream/1947/9119/1/WRE-CPD-TM-169%20PR.pdf

but it will be up to you how to interpret in terms of FFTs.

Also something called a modulation processor in:

Nielsen, Richard O. Sonar signal processing. Artech House, Inc., 1991.
which is unfortunately out of print, but useful if you can find it.

If your theoretic interest needs to be scratched,

R. O. Nielsen, "Cramer-Rao lower bounds for sonar broad-band modulation parameters," in IEEE Journal of Oceanic Engineering, vol. 24, no. 3, pp. 285-290, Jul 1999.
doi: 10.1109/48.775290
keywords: {amplitude modulation;frequency estimation;phase estimation;signal classification;sonar detection;sonar signal processing;Cramer-Rao lower bounds;DEMON processing;amplitude modulated broadband signal;blade frequency;discrete frequency lines;low SNR;modulation frequency;modulation level;modulation phase;parameter estimation;passive sonar signal processing;propeller noise;propulsion shaft;rotational frequency;sonar broadband modulation parameters;source classification;Amplitude estimation;Amplitude modulation;Frequency estimation;Frequency modulation;Parameter estimation;Phase estimation;Phase modulation;Propellers;Signal processing;Sonar detection},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=775290&isnumber=16844

Don't expect an easy time.  A lot of tweaking on real data will be required.

Using "Modulation Processor" as a search term as opposed to DEMON would help Google.


Fluctuation of autocorrelation of a signal due to signal's noise

If you view $V(t)$ as desired signal $kt$  affected by noise $R(t)$, the clearly as you keep increasing $k$, this has the effect of increasing signal energy. Hence if the noise power is kept constant then $V(t)$ tends towards being more deterministic (in a loose sense) with improving SNR, and autocorrelation hence is more deterministic with reduced fluctuations.


Magnitude-squared Coherence calculation inconsistence

For any single chunk (window) of data the coherence will, as you observed, be 1. In order to properly estimate coherence you must average the spectra and cross-spectra for multiple windows, and THEN calculate coherence.

The auto-spectra XX and YY can be averaged the conventional way. For the cross-spectrum XY you must average the real and imaginary components first before calculating XY = sqrt( XY[imag]^2 + XY[real]^2 ).

Does that help? Averaging over 8 windows usually yields reliable estimates.


Nyqvist-shannon sampling theorem  only says that your sampling frequency should be greater than twice the bandwidth of the signal and not twice the carrier frequency of the signal. Hence you can modulate your 20 MHz signal at any carrier frequency you need and still get back your original signal by sampling  at greater than 40 MHz. Hence you need not sample at 2GHZ for recovery of signal, you can do away with a sampling frequency greater than 40 MHz itself  For a more clear understanding refer to bandpass sampling.

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


Questions about reproducing results of a paper on Communication Systems and Wireless Communication

Let's tackle this piece by piece:

Path loss

Your path loss model gives the path loss in dB. Convert to linear scale. You can directly read $\alpha$ and $G$ from the result.

Haven't read that paper (paywalled), but it looks to me like your shadowing should be part of the fading  $\beta$. Make sure your channel model is compatible with the one used in the paper at all! There's no guarantee that different sources use channel models that are.

Power Spectral Density to Power

Well, assuming white noise (check that assumption is correct for this paper!), the noise power in a channel is of course bandwidth of that channel times noise power spectral density.

  should I do that in the dB or Watts units ?

dB is not a unit!

Generally, the actual units don't matter, as long as you do all the calculations correctly, so it really doesn't matter whether you're working in W or (MeV per minute).

Your first question, however, points to you not being 100% familiar with what dB and dBm are. I strongly recommend reading the Wikipedia article on decibel if that is the case!


Why is pitch detection necessary in pitch correction?

I don't know what algorithm you are talking about exactly, as you didn't state that in your question, but my idea of pitch correction consists of detecting the actual pitch and shifting it to the right one.

You ask why pitch detection is necessary. Well, to be able to write a program that automatically corrects the pitch of a recording (i.e. in opposition to a scenario where you are manually checking what frequency is the right one for each instant of time) you should first have an array containing all possible frequencies corresponding to musical notes. That is, a list containing all the values that would be correct if detected, so that they wouldn't have to be changed. If the singer has accurately hit a $\mathrm{C_3}$, then why would you auto-tune it?

Then, you should perform pitch detection using STFT. For each instant of time, there will be a dominant frequency. One (primitive) way to do this would be to find the frequency that shows the maximum amplitude. When you find that maximum (pitch detection), you should compare the frequency where it's located to the array of musical notes you created before. Then, assuming the singer is not that bad and the tone he's hitting is approximate to the desired one, one should take the detected frequency to the closest musical note (pitch shifting). For example, if the maximum is detected at $\mathrm{124 \ Hz}$, one would go for a $\mathrm{B_2}$ ($\mathrm{123.47 \ Hz}$) instead of a $\mathrm{C_3}$ ($\mathrm{130.81 \ Hz}$). I believe there are some things about logarithmic perception that should be taken into account here, but that's not relevant in this question. Take into account that there are pitch detection algorithms, I just used this easy method as an example to make it illustrative.

That's why detection has to be made. You have no way to automatically "round" the pitch to the right one if you don't know what actual pitch the singer is in.


Estimates sub-pixel shift directly in phase region

I am unable to comment on your Matlab code, but eq. (14) seems straightforward to me.  You have a phase difference field $\theta$ which depends on the two spatial wavenumber components $k_1$ and $k_2$, and which formes a wrapped plane.  The shift you seek is the slope of this plane, expressed as two scalar components $a$ and $b$.  Since phase wraps, the slopes are decomposed into integer and fractional parts, $a = a' + a''$ etc.  In eq. (14) you obtain the residual phase field by subtracting the phase due to the integer part $(a', b')$.  You need the actual $k$ values when you compute $(a'k_1 + b'k_2)$, not just the $k$ indices.  Computing the $k$ values is similar to making a frequency axis for a 1D Fourier spectrum analysis.


To expand on Justme's reply:

From the looks of your code you seem to have confused what constitutes an LTI system.

In your code the LTI system that I see (y1=2*z) scales an undelayed input signal by 2. The additional 2*cos(pi/t) is NOT part of the system, it is in effect just another sinusoid that you have added to your signal - thus giving it an additional frequency component.

An LTI system can only ever do two things to an input complex signal - it can scale the magnitude and it can modify the phase of frequency components of the signal. It cannot modify the fundamental frequencies of the sinusoidal components of a signal. (NB: It can however create resonances that self-oscillate - this is however cause by recursion in the system)

Returning to your code once again, we see that you have checked the linearity property, but have also included a scaling of the additional sinusoid that you have misattributed as part of the system. If you remove this sinusoid in your system function, you will find that the linearity property holds for your system

Edit: Fixed some technical wording issues


Try to remember what happens when you multiply a signal with a cosine function with a frequency of Fc.

Namely, what's the shifting operator in Frequency domain?

If we're talking on visualization, so the operator in MATLAB is called fftshift()

Enjoy...


jpg is a compressed storage format and no matter how binary your image is, you will almost always end up with gray values in the resulting image (due to compression).
Please save it as png and try again.

(Also note that you still have 0s and 255s as the dominant bins, because jpg generally degrades around the edges, which constitute a small portion in the image.)


Properties of convolution same as LTI systems?

The author simply means that the response of an LTI system $h[n]$ can be found by convolving the input sequence $x[n]$ with the $h[n]$. And in doing that all properties of convolution will hold true for LTI systems. Like, frequency response of an LTI system will simply become $X(e^{j\omega}).H(e^{j\omega})$, because convolution in time means product in $\omega$ domain.

You can look through all the properties on convolution operation at : ConvolutionProperties, and be sure that all of them will be held by an LTI System as well.

For example : Associativity : $$(x*h_1)*h_2 = x*(h_1*h_2)$$
In context with LTI systems, it means that if $x[n]$ is the input to two cascaded LTI systems $h_1[n]$ and $h_2[n]$, then we can actually, find out the equivalent LTI system by convolving the 2 cascaded systems as $(h_1*h_2)[n]$ and then find the output of this equivalent system by applying input $x[n]$ to it.

Similarly, you can see the relation between all of the convolution properties in context of an LTI system.


As this page https://ccrma.stanford.edu/~jos/ReviewFourier/FFT_Convolution.html mentions, when you perform cyclic convolution (FFT convolution) you have to add zeros to your signal $A$ to the length $$N_Y = N_A + N_A -1,$$ where $N_Y$ is the length of the output convolution and $N_A$ is length of your signal. In your case you are performing the FFT Correlation, which is almost the same operation and you should extend your signal $A$ with zeros as well.

So for example:

a = np.concatenate((a,np.zeros(len(a)-1))) # added zeros to your signal
A = np.fft.fft(a)
S = np.conj(A)*A
c_fourier = np.fft.ifft(S)
c_fourier = c_fourier[:(c_fourier.size//2)+1]

After that you should get the same results as for numpy auto-correlation.


fvtool is interactive so it doesn't let you recover an data in a script. Fortunately you can use the freqz function to get the same frequency response. Something like this might work:

[h,w] = freqz(b,a);  % h is frequency response, w normalized frequency
mag = 20*log10(abs(h));
[minval,minidx] = min(mag);
minfreq = w(minidx);
[maxval,maxidx] = max(mag);
maxfreq = w(mcxidx);


Bandwidth of the Derivative of a Message

Hint:

Checking a table of Fourier transform properties you'll find that the Fourier transform of the derivative is given by

$$\mathcal{F}\left\{\frac{dm(t)}{dt}\right\}=j\omega M(\omega)\tag{1}$$

where $M(\omega)$ is the Fourier transform of $m(t)$. Now if $M(\omega)$ is band-limited, what does this tell you about the Fourier transform of $dm(t)/dt$?


Calculating the first derivative of an image using DFT

They are not the same. Using 1D notation,

the discrete-time (backward) first difference is $x[n] - x[n-1]$ whose frequency domain DTFT equivalent is

$$ x[n]-x[n-1] \leftrightarrow X(e^{j\omega}) - e^{-j \omega} X(e^{j\omega}) =X(e^{j\omega})(1- e^{-j \omega})  $$

which becomes
$$ x[n]-x[n-1] \longleftrightarrow X[k](1 - e^{-j \frac{2 \pi}{N} k})$$ using the DFT to implement it.

The FIR impulse response of the discrete-time system that implements the first difference is therefore,

$$ h[n] = \delta[n] - \delta[n-1]$$

The first derivative of a continuous-variable function $x(t)$ is $x'(
t)$ and in CTFT domain it becomes :
$$ x'(t) \longleftrightarrow j\Omega X(\Omega) $$

where the analog system frequency response is

$$H_d(\Omega) = j \Omega $$

which is not implementable in digital form, but a bandlimited approximation to it is attained under a sampling period of $T$ that yields an equivalent discrete-time frequency response of a (bandlimited) differentiator as

$$ H_d(e^{j \omega}) = j \frac{\omega}{T} ~~~, ~~~~ |\omega| < \pi$$

The associated discrete-time (IIR) impulse response is
$$ h_d[n] = \text{IDTFT} \{ H_d(e^{j \omega}) \} = \text{IDTFT} \{ j \frac{\omega}{T}  \} $$

Practically you would truncate and window $h_d[n]$ before using.

So they are not the same.


How do I quantify the amount of information being lost by filtering a digital signal?

what metric can I use (or even theory or branch of analysis if the answer isn't a simple one) to determine the amount of filtering necessary?

You can use what you already have, that is, the "strengths" of the components that result from the analysis. You can sort the coefficients in decreasing order and then obtain their cumulative sum. Finally, you can set a threshold and stop adding components once the cumulative strength of the chosen components goes (for example) above 90%. In this way you would be selecting the $N_{sc}$ strongest components and rejecting components that have higher signal-to-noise ratios (or in other words are closer to some background noise level). If this has to be a low pass filter, then after you have isolated the first $N_{sc}$ components that contribute to 90% (for example) of the total signal strength then you can set the low pass filter cut off frequency at the component with the highest frequency (or bin).

  For any given analysis window of say 500 samples, how do I determine programatically how to preprocess my signal before trying extract features?

Please see previous answer, given the frequency (or scale) limits imposed by the length of the window (here 500 samples), you might have to adjust thresholds and cut offs. Please note that if you try to apply the above algorithm per frame (of 500 samples each), then, inevitably, you will get some sort of modulation between frames which, depending on your application you will have to manage somehow. This is because the signal strength will be changing between windows of 500 samples, which means that the potential cut off of the low pass filter will be changing as well, which means that each frame will bear different amounts of filtering. The easiest way to deal with this would be to increase the window overlap (at the expense of more frames and therefore more calculations).

  Given a sample window with a particular spectral distribution, does an optimal cutoff for a bandpass filter exist? (or a notch frequency)

Only under specific assumptions about the spectral distribution of the noise. In the absence of this, you can study your signal and try to create a Matched Filter that is tailored to the spectral characteristics of the signal and therefore anything outside that would be filtered out.

  if so, does it always exist?

Please see previous answer

  if an optimal value exists, how do I obtain it?

Please see previous answer

  What is the relationship between this and the amount of information in the signal?

This "Amount of information" is a bit unclear. Information has a specific meaning, so are you refering to the amount of bits that are expected to be able to encode completely your "clean" waveform? Or entropy derived by the spectrum maybe? Or spectral flatness? Or, are you talking about the actual information that this pulse sequence would carry if it was to be seen as a communications signal (?). As a general comment, the low pass filter would lower the entropy of the signal in its output...but this is too general a statement to be useful in your context perhaps (?).

Hope this helps. Happy to ammend the response if more data becomes available.


Can linear transforms of vector space be seen as LTI systems?

The question of abstract nature so a rigorous answer could require more elaboration, however on the surface one can at least observe the following.

An LTI system is also a linear transform, which maps its input signal (vector) $x(t)$ or ($ x[n]$) into an output signal (vector) $y(t)$ (or $y[n]$), where both the input and output belong to the same signal space, one kind of vector space.

For such transforms and vector spaces the time invariance property can be checked easily such as $y(t-d) = \mathcal{T} \{x(t-d) \}$, but for other vector spaces such as the n-tuples or matrix spaces etc, the operation of time invariance is not trivial to check, at least intuitively. The shift of a vector (an element of the set) should be defined first.

May be the following links could help further:

link1

link2


Good Reference Problem to Test Filtering/Estimation Algorithms

This answer might give you one example problem and a KF implementation to solve it.

The usual problems are kinematics problems, where you want to track position, velocity, and acceleration in one or more dimensions.

There are some nonlinear problems, too, like frequency / phase tracking but that might be overkill. See how it works on linear problems first.


Detecting changes between two similar images

The first part of the problem that you're describing is called registration.
Assuming that the displacement is indeed small, and your problem has no time limits, you can use normalized cross-correlation on the images to find the displacement. It will also take care of amount of light normalization. (That means that you are assuming that the transformation between the two images is only displacement.)

After you found the displacement, you can normalize both of the images by their mean intensity, and do some simple pixel-by-pixel comparison.

Ideas for improvement and variations:

Use keypoints descriptors and matching to find transformation between two images, to save runtime.
Use different models of transformation for different scenarios - rotation for in-plane camera rotation that is orthogonal, scale for zooming scenarios, perspective for out-of-plane rotations.
Use smarter comparison between the two registered images, for example, only in the places where there are edges.
Create a noise model by checking your camera characteristics in order to know what is a significant change. For example, more than two standard deviations of noise.


Finding specific spots in an audio signal

What you are describing is exactly how to plot a spectrogram. Finding the spots can be done with a classic peak-finding function in any software.

If you are looking for coding the algorithm from scratch, here is a method, which is the most intuitive but not the fastest in term of calculations :

Consider your spectrogram as a $n\times m$ array of numbers, with each number position referenced as $(i,j)$. You have to find which points are greater or equal to all of their neighbors.

In pseudo-code it would give something like :

for i from 1 to n
  for j from 1 to m
     if :
       number(i,j) > number (i-1,j-1)
   and number(i,j) > number (i-1,j)
   and number(i,j) > number (i-1,j+1)
   and number(i,j) > number (i,j-1)
   and number(i,j) > number (i,j+1)
   and number(i,j) > number (i+1,j-1)
   and number(i,j) > number (i+1,j)
   and number(i,j) > number (i+1,j+1)
     then :
       save (i,j) as a peak value
   end if
end for
end for

Of course, this is not very sharp since it does not consider whether you are finding a "big" peak or just a smaller one. In that case, you can improve it by setting a minimum offset value like :

if number(i,j) > number(i-1,j-1) + offset

Note that this is just a simple example of a peak finding algorithm and I didn't took the boundary points specific case into account. More efficient methods exist but I can detail each one here.


Duration of downsampled signal in analog frequency axis

If you plot spectrum of y[n] along y-axis and frequency(in Hertz or in any other suitable  unit) along x-axis, there will be infinite copies of the frequency scaled(by a factor of M) version of the original spectrum(i.e, spectrum of original signal x(t)) each around kFs, where k is an integer and Fs = sampling frequency.


If you would write the frequency-response in terms of projections through inproducts. How would you do it?

the impulse-response and frequency-response of a signal

Only systems have impulse responses and frequency responses, signals don't. I assume that what you mean here.
A System  describes the relation ship between its input signal and its output signal. For an LTI system, that relationship can be captured through either the transfer function or the impulse response. The idea is similar

Project the signal onto the some orthogonal basis
Run the basis function through the system
Assemble output basis function into the output signal. Since it's an LTI system you can simply sum or integrate.

Of course that only makes sense if the bases functions are "easier" to run through the system than some arbitrary signal. For the transfer functions the basis functions are complex exponentials and for the impulse response they are time shifted dirac deltas.
Update
Elaborating on step 2: The basis functions have an index or parameter. For example for a complex exponential, the parameter is frequency. Now you calculate the output of the system with the basis function as a an input using the parameter as a variable. Specifically you calculate the output for a complex exponential with frequency $\omega$  as a variable. The result is a function of frequency: That's the transfer function.
So if you have a signal $x(t)$ that can be expressed as the sum of some basis functions b_k(t) like
$$x(t) = \sum a_k \cdot b_k(t) $$
You can calculate the output as a weighted sum of the outputs of the basis functions.
$$y(t) = T\begin{Bmatrix}
x(t)
\end{Bmatrix} = T\begin{Bmatrix}
\sum a_k \cdot b_k(t)
\end{Bmatrix} = \sum a_k \cdot T\begin{Bmatrix}
b_k(t)
\end{Bmatrix} $$


Exponential decaying step response in LTI System

You're forgetting the discontinuity at $t=0$. Your step response is actually

$$Y(t)=(1+e^{-0.1t})u(t)\tag{1}$$

where $u(t)$ is the unit step function. Differentiating $(1)$ gives

$$h(t)=-0.1e^{-0.1t}u(t)+(1+e^{-0.1t})\delta(t)=-0.1e^{-0.1t}u(t)+2\delta(t)\tag{2}$$

where $\delta(t)$ is the Dirac delta impulse. Convolving $(2)$ with a unit step gives the original step response:

$$\begin{align}\int_{-\infty}^{\infty}h(\tau)u(t-\tau)d\tau&=-0.1u(t)\int_{0}^{t}e^{-0.1\tau}d\tau+2\int_{-\infty}^{t}\delta(\tau)d\tau\\
&=u(t)(e^{-0.1t}-1)+2u(t)\\&=(1+e^{-0.1t})u(t)\end{align}$$


Why can the FFT always be mirrored in the middle of the x axis?

The DFT of a real signal is conjugate symmetric. For example, if your DFT result at, say, 2Hz was $1+j5$, then your DFT result at -2Hz would be $1-j5$. This is conjugate symmetry.

Of course, when you take the absolute magnitude, the result in both cases is the same, which is why you see this mirror image.

Going deeper, the reason for why a real signal can be decomposed into symmetric complex conjugate parts, is because a real signal oscillates on the real axis, and this can be thought of as the resultant of two phasors rotating around the center of the complex plane, in opposite directions.

Hope that helped.


Filtering a digital signal online in real-time using python

Do I have to filter the whole (or at least a huge bit) of the signal every time a few new samples came in or is there a way (like the sliding DFT) where it is possible to efficiently determine the new part of the filtered signal?

Digital filters don't work like that -- basically, classical FIR or IIR can work on every single new sample. You should really read up on what these filters are, and how people model them.

  I would like to use a butterworth filter

Well, there's plenty implementations of that out there,

  I am currently using scipy's butter and lfilter

of which you already know one!

Now, a butterworth filter is a recursive thing, so to calculate the next part of your sampled signal, you'll need the last state. That's exactly the "filter delay state zi" that lfilter returns and can take in on the next call as zi parameter.

  but I don't know how to use it to get a constant signal.

I think you mean "to achieve continuous filtering".

Now, that being said, the point is that you're setting yourself up for writing your own streaming architecture. I wouldn't do that. Use an existing framework. For example, there's GNU Radio, which lets you define signal processing flow graphs in Python, and also is inherently multithreaded, uses highly optimized algorithm implementations, has a lot of in- and output facilities, and comes with a huge library of signal processing blocks, which can be written in Python or C++, if you happen to need to do that.

For example, a flow graph that takes in samples from a sound card, butterworth-filters them and writes them to a file is:

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
##################################################
# GNU Radio Python Flow Graph
# Title: Butterworth Test
# Generated: Mon Feb  8 16:17:18 2016
##################################################

from gnuradio import audio
from gnuradio import blocks
from gnuradio import eng_notation
from gnuradio import filter
from gnuradio import gr
from gnuradio.eng_option import eng_option
from gnuradio.filter import firdes
from optparse import OptionParser

class butterworth_test(gr.top_block):

    def __init__(self):
        gr.top_block.__init__(self, "Butterworth Test")

        ##################################################
        # Variables
        ##################################################
        self.samp_rate = samp_rate = 48000

        ##################################################
        # Blocks
        ##################################################
        # taps from scipy.butter!
        self.iir_filter_xxx_0 = filter.iir_filter_ffd(([1.0952627450621233e-05, 0.00013143152940745496, 0.0007228734117410033, 0.0024095780391366808, 0.005421550588057537, 0.008674480940892064, 0.010120227764374086, 0.008674480940892081, 0.005421550588057554, 0.0024095780391366955, 0.0007228734117410089, 0.00013143152940745594, 1.0952627450621367e-05]), ([1.0, -4.4363862740719835, 10.215121830052535, -15.374408118154847, 16.57333784740102, -13.325056987818655, 8.133543488903097, -3.77641064765334, 1.3181452681671835, -0.3361758629961047, 0.05930166356243964, -0.0064815521348275, 0.00033130678123743994]), False)
        self.blocks_file_sink_0 = blocks.file_sink(gr.sizeof_float*1, "", False)
        self.blocks_file_sink_0.set_unbuffered(False)
        self.audio_source_0 = audio.source(samp_rate, "", True)

        ##################################################
        # Connections
        ##################################################
        self.connect((self.audio_source_0, 0), (self.iir_filter_xxx_0, 0))
        self.connect((self.iir_filter_xxx_0, 0), (self.blocks_file_sink_0, 0))

def main(top_block_cls=butterworth_test, options=None):

    tb = top_block_cls()
    tb.start()
    try:
        raw_input('Press Enter to quit: ')
    except EOFError:
        pass
    tb.stop()
    tb.wait()

if __name__ == '__main__':
    main()

Note that this code was auto-generated from a graphical flow graph that I just clicked together using the gnuradio-companion program:

If you want to know more about how to implement signal processing flow graphs in Python, go for  the GNU Radio Guided Tutorials.

EDIT: I liked @Fat32's answer pretty much! What he describes as a double buffering architecture is pretty close to what GNU Radio does:

An upstream block produces samples in sample chunks of arbitrary sizes, writes them to the output ring buffer (which is represented as an arrow in above picture), and notifies its downstream blocks that there's new data.

The downstream block gets notified, checks whether there's enough space in its output buffer to process the samples that are in his input ring buffer (which is the same as the output buffer of the upstream block), processes these. When finished, it informs the upstream block(s) that it has used up input ring buffer (which can then be reused by the upstream blocks as output), and the downstream blocks about new samples being available.

Now, GNU Radio being multi-threaded, the upstream block might already be producing samples again; in a normal GNU Radio application, nearly all blocks are "active" simultaneously and things scale pretty well on multi-CPU machines.

So, GNU Radio's main job is to give you this buffer infrastructure, the notifying and thread housekeeping, the clear signal processing block API and something to define how everything is connected, so you don't have to write what Fat32 describes in her/his post yourself! Notice that doing sample stream marshalling isn't that easy to do properly, and GNU Radio takes the hardness out of it and lets your concentrate on what you want to do: DSP.


Why do the DTFT and FFT give me completely different results for magnitude at a specific frequency?

HINT:
Try introducing the sampling frequency $F_s$. Then $f$ in the relation in Equation $(1)$ is relative to the sampling frequency
$$
\omega = 2\pi f\quad\text{where}\quad f = \frac FF_s\tag{1}
$$
Where the units are as follows:

$\omega$ in [radians/sample]
$f$ in [cycles/sample]
$F_s$ in [samples/second]
$F$ in [cycles/second] or [Hz]

Your discrete time vector for your signal should be as shown in equation $(2)$
$$
t = nT_s\quad\text{where}\quad T_s = \frac 1F_s\quad\text{and}\quad n = 0, \ldots, N - 1\tag{2}
$$
You seem to use the frequency $F$ in [Hz] instead of the relative frequency $f$ in [cycles/sample].


What are the properties of continuous-time improper systems?

Note that the imaginary axis is not inside the ROC because $|H(j\omega)|\to\infty$ for $\omega\to\infty$. So the system is definitely unstable, regardless of causality, which confirms your intuition.

Concerning causality, the ROC is neither a right half-plane nor a left half-plane (both because of the pole at $s\to\infty$), which is consistent with a (generalized) function concentrated at $t=0$.


What is the true meaning of a minimum phase system?

The relation of "minimum" to "phase" in a minimum phase system or filter can be seen if you plot the unwrapped phase against frequency.  You can use a pole zero diagram of the system response to help do a incremental graphical plot of the frequency response and phase angle.  This method helps in doing a phase plot without phase wrapping discontinuities.

Put all the zeros inside the unit circle (or in left half plane in the continuous-time case), where all the poles have to be as well for system stability.  Add up the angles from all the poles, and the negative of the angles from all the zeros, to calculate total phase to a point on the unit circle, as that frequency response reference point moves around the unit circle.  Plot phase vs. frequency.  Now compare this plot with a similar plot for a pole-zero diagram with any of the zeros swapped outside the unit circle (non-minimum phase).  The overall average slope of the line with all the zeros inside will be lower than the average slope of any other line representing the same LTI system response (e.g. with a zero reflected outside the unit circle).  This is because the "wind ups" in phase angle are all mostly cancelled by the "wind downs" in phase angle only when both the poles and zeros are on the same side of the unit circle line.  Otherwise, for each zero outside, there will be an extra "wind up" of increasing phase angle that will remain mostly uncancelled as the plot reference point "winds" around the unit circle from 0 to PI.  (...or up the vertical axis in the continuous-time case.)

This arrangement, all the zeros inside the unit circle, thus corresponds to the minimum total increase in phase, which corresponds to minimum average total phase delay, which corresponds to maximum compactness in time, for any given (stable) set of poles and zeros with the exact same frequency magnitude response.  Thus the relationship between "minimum" and "phase" for this particular arrangement of poles and zeros.

Also see my old word picture with strange crank handles in the ancient usenet comp.dsp archives: https://groups.google.com/d/msg/comp.dsp/ulAX0_Tn65c/Fgqph7gqd3kJ


SQNR when multiplying 2 16 bit signed integers and truncating result?

When you multiply the two 16 bit numbers you get the "true" 32 bit number and then you truncate to get a 16 bit number.  The noise that is introduced when you truncate is a random value from $0$ to $2^{16}-1$ ($65535$), which averages out to $32767.5$ and has an average power (average of the square) of $1.43E9$.  If you were rounding instead of truncating it would be a random value from $-(2^{15})$ to $2^{15}-1$, which averages out to approximately $0$ (no bias), and has an average power of $3.58E8$.  That's why it's better to round than truncate, but that's a story for another day.

The truncation noise, therefore, is $1.43E9$.  You get the signal power by multiplying your truncated data by $2^{16}$ (to make the 16 bit values into 32 bit values), squaring them, and then calculating the mean.  You then calculate the SQNR by dividing the signal power by the truncation noise, $1.43E9$.


Compressive Sensing: Reconstruct Gap in Antenna Array

In general: no. What you can and cannot reconstruct will depend a lot on your model.

If you assume that you have:

point sources in the far field of the antenna array
these point sources emit narrow-band, mutually uncorrelated signals
you can observe them over a long window in which they are stationary
then you can use co-array processing to virtually fill the missing elements between your two arrays.

Let's see how it works for a single source and let us emit noise for now for clarity (it's not hard to add it later). The observation vectors in your two subarrays can be expressed as $\mathbf x_i(t) = \mathbf a_i \cdot s(t) $, $i=1,2$, where $\mathbf a_i$ are the array steering vectors and $s(t)$ is the emitted signal. Let $\mathbf x(t) = [\mathbf x_1(t)^T, \mathbf x_2(t)^T]^T$ be the vector of all observations. Then, its covariance matrix for a single source is given by $$\mathbf R = \mathbb{E}\{\mathbf x(t) \mathbf x(t)^H\} = \mathbf a \cdot \mathbf a^H \cdot P_s, $$ where $P_s$ is the variance of $s(t)$ (remember it was assumed stationary) and $\mathbf a = [\mathbf a_1^T, \mathbf a_2^T]^T$. Now, the trick is to vectorize $\mathbf R$. Using the Kronecker-vectorization identity, this gives $$\mathbf r = {\rm vec}\{\mathbf R \}= (\mathbf a^* \otimes \mathbf a) \cdot P_s,$$ where $^*$ denotes complex conjugation and $\otimes$ is the Kronecker product. Now, this looks like a single observation with an array that has an array steering vector given by $\mathbf a^* \otimes \mathbf a$, which we refer to as the co-array of $\mathbf a$.

For a ULA with half-wavelength spacing, your array steering vectors take the form $\mathbf a_1 = [{\rm e}^{0}, {\rm e}^{\jmath \mu}, {\rm e}^{2\jmath \mu}, \ldots, {\rm e}^{(N-1) \jmath \mu}]^T$ and $\mathbf a_2 = [{\rm e}^{D \jmath \mu}, {\rm e}^{(D+1)\jmath \mu},  \ldots, {\rm e}^{(D+N-1) \jmath \mu}]^T$, where $\mu = \pi \sin(\theta)$ is the electric angle. Therefore, $\mathbf a^* \otimes \mathbf a$ contains terms ${\rm e}^{p \jmath \mu}$, where $p=-(N-1), \ldots, (N-1)$ (from $\mathbf a_1^* \otimes \mathbf a_1$ and $\mathbf a_2^* \otimes \mathbf a_2$) and $p=-D-N+1, ..., D-N+1$ and $p=N-D+1, ..., N+D-1$ (from $\mathbf a_1^* \otimes \mathbf a_2$) [give or take a $\pm$, I guess you get the idea].

Given that, for $D<N$, you've created a virtual ULA that is completely filled. I'll let you do the math how many elements it has in total.

Note that you will in general need to apply spatial smoothing if you want to estimate more than one source.

That said, compressive sensing is the wrong term here. It is co-array processing, which is somehow linked to sparsity-based processing (which became popular in the compressive sensing field, leading to many people confusing the two).


how to detect short pause in a speech with noise?

You can use Voice Activity Detector(VAD) for detecting the pauses in speech and there location(and duration). if your input signal is not very noisy and noise is not varying much you can use a fixed threshold on the energy(calculated for short frames, ex-20 ms), so if energy is above that threshold you declare that frame as speech else pause. if input signal is very noisy(varying noise) you can use adaptive threshold and do the same,i.e. you can have short term(lasting for 10-12 frames) and long term energies(lasting 40-50frames) and long term energy can be used as threshold. once this classification is done for entire input signal you can get duration's of pauses in speech


You should consider doing an STFT on the accelerometer signal. This will allow you to visualize the frequency of the signal. If you can determine the frequency of the brakes, you should be able to set a threshold for that specific braking frequency and count each time it is exceeded. You would also be able to measure how long the braking occurs.

How are you processing your signal? Matlab and python have the specgram function, which will allow you to see the actual frequency braking. From there you will be able to determine the brake frequency as well as a threshold.


I need to preprocess raw ecg data in R, here is a sample already standardized. I'm not an expert in signal processing nor experienced in working with medical data,...

Not being an expert on how the heart works and how its phases manifest themselves on the ElectroCardioGram (ECG) is not a problem. But it would help immensely if you mentioned what is the purpose of this, i.e what do you want to do, what are your specifications and constraints?

For example "We put the driver in the car, they go through the track a few times and then we want to know what was their heartbeat at a specific point in the track".

Excellent, then, do not even use ElectroCardioGraphy (ECG) because:

The electrical environment of the racing car is incredibly noisy already. For comparison, you are trying to measure a signal that has a max amplitude of a few mV through long leads that have to run from the chest of the driver to the ECG box. Just a short distance away from this long electrical loop, you have sparks going on and off, relays turning big loads on and off, electrical motors or stepper motors turning. All of these things create fields that the loop will pick up and they will contaminate your signal.
ECG measures the electrical activity of the heart via contact electrodes. Which means that the skin of the driver would have to be prepared, ideally shaved, cleaned with alcohol, then have a conductive gel applied to it and finally apply the electrodes and fix them with medical tape. Throughout the track, I would not expect the driver to pay attention to the electrodes or, to not sweat. Even in controlled conditions, if you tap the electrodes with your hand, you can observe changes in the signal.

Instead, use Pulse Oximetry. This is the way that most gadgets (e.g. smart watches, fitness wrist bands, etc) obtain Beats per Minute from the heart today. Pulse oximetry needs to have contact with the skin, but not because it measures electricity. Very briefly: It measures the attenuation of light that is proportional to the blood flow below it. Because the heart is pumping blood in a rythmic way, the flow is rythmic as well and this signal can be used to derive a relatively accurate BPM count.

It will also be much easier to cross reference this with the rest of the telemetry because it will simply be a time series of BPM measurements. Pulse oxymeters for people to manage their heart conditions are more accurate than gadget type instruments. But then again, they might get a bit of getting used to by the driver. (For an example, see this)

  Below is a plot of how the data looks like,

No way.

Here (and here ) is what typical ECG signals look like.

To extract BPM, you need to "detect" those spikes which correspond to heart contractions and then measure the timing between them (also known as R-R interval, because the main spike of the ECG is the R peak).

Very briefly: QRS Detectors work by Pre-filtering -> Rectification -> Integration -> Thresholding -> QRS Complex pulse sequence -> R interval extraction -> Beat Detection. For much more detailed information about this, have a look at Moody's work.

What is depicted in the plot is saturation and possibly lots of noise from the contact. There is no trace of "typical" ECG complexes to work with.

  Question 1) What strategy is needed in terms of subsampling, filtering to preprocess the data to clean get a clean ECG signal (google how it should look).

Why subsample? The signal is already in a bad state, why make it even worse? You could try a bandpass filter between 0.1 - 100Hz. This would reduce DC and cut out any high frequency noise but you have a lot of spikes which means that most of time times you will be "hearing" the ringing of the filter.

  Question 2) Would this strategy work for EMG data also?

EMG is even worse because it picks up signals at higher frequencies. What is the application?

If you are trying to see the timing of the driver as to WHEN do they start to accelerate (i.e. leg muscles start working) versus when the car actually accelerates or brakes (i.e. break signal from telemetry), then maybe you can get away with very simple processing using thresholds. But if you are trying to make more fine measurements, you might be back into the same problem of lots of noise from various sources. In the case of EMG you will have to work around those, which means, twisted pair wires, shielding, as short lengths as possible, very good contact with the driver's skin, very good preparation of the driver's skin, perhaps even slightly different electrodes that attach on to the skin rather than simply be in contact with it.

Hope this helps.

EDIT:

I am sorry, this is an unsalvageable situation. I have a couple of comments only,but I do not believe they will be very much of help.

So, here is the full signal you shared (thank you, by the way):

As you can see, there is nothing usable here. My perception is that the ECG amplifier is either saturating (because of induced voltage) or the contact with all three electrodes is not stable. If you notice a Holter monitor, which is used to "catch" abnormal heart functioning wherever it may happen, you will see that the leads are as short as possible with the device being mounted ON the chest. If we were to use the same device in an application similar to yours, we would have to use shielded leads and make sure the pads are firmly attached on the body.

Going a bit closer, it almost looks like "reversed bits"

What I mean by this is: Say for instance, you created your own ECG circuit. So, take an off the shelf instrumentation amplifier, set the gain and then feed that to an "Arduino"ish type of board with some 10bit Analog To Digital (ADC) converter that establishes a serial like channel to a pc or stores values internally. Because you have 10bits (or more) and you would have to send two bytes over the wire (or store two bytes) you would have to decide if they are sent/stored/received/retrieved Big Endian or Little Endian.

So, when you read a binary file from the disk and you get the endianess wrong, your integer values, the levels recorded by the ADC, would look a little like this second figure right here. There are glimpses of what looks like a signal inter dispersed with wild swings from the minimum to the maximum value. Also, your "zero" is not at 0, but maybe you used just two ECG leads.

I am just leaving this here as a comment in case it is what is going on, though I doubt it.

If you look more closely, the signal is just all over the place:

In fact, if you have telemetry from the rest of the car, try to plot together the RPM of the engine and the ECG signal that you have. What I would expect to see is the "density" of disturbances in the signal to correlate with the RPM of the engine, assuming that most of this interference is coming from the electrical system of the car (ignition, injection, throttle stepper motors, etc).

Hope this helps.


Normalisation of different length signals

In audio, "normalization" means making all values be between 0 and 1. So

recording = recording / max( abs(min(recording)), abs(max(recording)) );


Can I take piecewise FFT of a larger input sample?

You need the Cooley-Tukey algorithm, which can be used to express the DFT of size $N=N_1N_2$ by DFTs of sizes $N_1$ and $N_2$. It works a bit like what you sketched in your question, but you forgot the twiddle factors. This answer explains the details of the algorithm and it also includes a simple Matlab implementation.

EDIT (see comments): If you really want to compute the FFTs of blocks of consecutive data you need to split the computation of even and odd frequency indices (assuming the DFT length $N$ is even):

$$X[2k]=\sum_{n=0}^{N/2-1}\left(x[n]+x[n+N/2]\right)W_{N/2}^{nk}\\
X[2k+1]=\sum_{n=0}^{N/2-1}\left(x[n]-x[n+N/2]\right)W_N^nW_{N/2}^{nk}\tag{1}$$

where $W_N=e^{-j2\pi/N}$. Note that you can of course compute the DFTs of the length $N/2$ sub-blocks separately instead of adding (or subtracting) the blocks before the DFT (as in (1)), but this is computationally less efficient.


Time stretching allowing for pitch change

Resampling does appear to be the answer for the OP. The point is that you will not change the speed of the sample as you will not change the sampling rate by resampling just the number of samples, so therefore the time represented by each sample must change after resampling. An example will best clarify this :

We have 10 seconds of audio data sampled at 44.1 KHz.

This is a total of 441,000 samples over 10 seconds.

We resample to twice the number of samples by inserting zeros between each sample and then low pass filtering out the images created (standard interpolation techniques that will minimize distortion).

Now we have 882,000 samples.

The thing is we don't actually change the sampling rate! So therefore time must have stretched.

882,000 samples / 44,100 samples/sec =  20 seconds duration!

This will result in a time stretching changing all frequencies accordingly.

Most resampling techniques are described as changing the sampling rate but would be equally applicable to this. So if the resampling technique describes changing the sampling rate by $M$, we can equally run the new samples at the same sampling rate and time would have changed by the same $M$.


So you think that the energy of the windowed signal in the frequency domain equals to the original signal multiplied by a certain value related to RMS of the window function, which is wrong. According to Parseval's theorem, the total energy of a signal can be calculated by summing power-per-sample across time or spectral power across frequency:
$$
\sum_{n=0}^{N-1} |x[n]|^2 = \frac{1}{N}\sum_{k=0}^{N-1}|X[k]|^2
$$
Hence I'm going to prove your statement wrong in the time domain.
The energy of a windowed signal $w[n]x[n]$ is
$$
\sum_{n=0}^{N-1} |w[n]x[n]|^2 = \sum_{n=0}^{N-1} |w[n]|^2|x[n]|^2 \neq C\sum_{n=0}^{N-1} |x[n]|^2
$$
where $C$ is a constant. The above inequality holds true as long as the window function $w[n]$ is not a rectangular window.