Time delay estimation of two sensor signals

Comment: As a new user I am not allowed to add a comment to your question - hence this approach.

You mention that you do get a peak from the cross-correlation but that the lag estimation is incorrect.  This implies that you know the correct result - is from a calibration signal or ..? And what is the lag roughly?

Could you explain what the signals look like (sinusoidals for example)?  Also what is the SNR, sequence length, sampling frequency? Have you filtered the signals before sampling - if so how?

If the offset is incorrect it might be related to insufficient time resolution. Using a sampling frequency $f_s$ means that your time resolution is $T_\mathrm{s}=\frac{1}{f_\mathrm{s}}$.  If this time resolution is not sufficiently small you may get an estimate that does not agree with your expectations.

There may be other techniques but this depends on the signals you have - and also if you for example know that the lag is always a fraction of a period time of a periodic signal. Some more details could help with this.


Equal loudness curves: any good numerical approximation equations?

I believe that the actual standard (access to which of course must be purchased) specifies all the mathematics. I did find this one MATLAB file in file exchange:

ISO 226 Equal-Loudness-Level Contour Signal

If you're going to use it in a commercial applications though, the distribution license won't allow you to do so. Still, a good start if you don't have the standard handy.


Camera calibration and extrinsic parameters for perspective transformation

For your application, here are the steps I would follow:

Calibrate the camera and obtain the intrinsics - You have obviously done that and corrected for the distortion already. initUndistortRectifyMap only lets you to precompute the projection. Otherwise, you might as well use projectPoints each time. This is a detail and I would not care about that initially. But it's a good practice anyways.)
The main idea is to find the laser plane (with extrinsic calibration to the camera). To find the extrinsics:

a. Place a checkerboard in front of the camera and leave enough space on the board for the laser line to be projected.

b. Image the laser line projections together with the checkerboard.

c. Now we can compute two planes with respect to the camera: The plane joining the laser line with the camera center, and the plane that checkerboard lies in (we now the checkerboard).

d. Intersect the two planes to find the line in 3D space.

e. Collect all such lines from multiple images and fit a plane to all of them. This is the equation of the laser plane.
To make measurements in 3D, simply shoot lines through the observed laser points in 2D image and intersect the rays with the laser plane. Each intersection will give a single point in 3D space.

This way we also do not explicitly need $(R,t)$. Our extrinsics are rather the equation of the laser plane.

In case we use a point-laser, the equation also remains the same, as instead of lines, we can collect sets of points and still do the plane fit.

By the way, here is a toolbox and a recent further read.


Decimating Non-Uniform large time-series data

I probably wouldn't bother interpolating at all. You want to downsample by factor 128 anyway, so from that point of view, your raw data is already a good approximation of being uniformly sampled at 128Hz. Just go with low pass filtering and downsampling and you should be good. Try it on a small fraction of the data, to get an idea of the error you introduce by this approximation, my gut feeling says, it will be acceptable, i.e. small compared to the accelerometer error. If actual time span does matter, calculate the average sampling rate and use this value instead of 128Hz.

For handling large data sets with limited RAM, I would opt for generators in python (yield instead of return) This way, there is no need to keep all the data in memory, you just "stream" the data from file to file in (overlapping and  windowed) chunks of appropriate size.


How to get started analyzing audio signals for detecting cracks in a material?

Not really a programming question!

In any case if you need to read books to understand DSP and are trying to compare two non stationary processes I suggest you start reading right away because most steps I tell you won't make sense and you need to answer a host of questions even before anyone can answer this question.

Signal processing deals with lot of statistical phenomena and you have to work with it emperically many a times. Not to say there are no fixed rules but the step-by-step to process this signals is hard enough when data is stationary. Given that you have non stationary data, get reading and once you have a hang of it ask some pointed questions on dsp.stackexchange.com


How can I correctly plot an impulse_response() of a discrete transfer function?

So finally this is the way I got to plot it right:

Gp = ct.tf(100,[1,0,100])
ts = 0.05
Gz =ct.sample_system(Gp, ts)
sym.pprint(Gz)

t = np.linspace(0,14*ts,num=15)

T, yout = ct.impulse_response(Gz, t)
yout=yout.flatten()

fig, ax = plt.subplots()
ax.step(T, yout)
ax.yaxis.set_minor_locator(MultipleLocator(0.05))
ax.yaxis.set_major_locator(MultipleLocator(0.1))
plt.margins(0.1, 0.1)
plt.title('Impulse Response c.1)')
plt.xlabel('Time [sec]')
plt.ylabel('Amplitude')
plt.grid(True)
plt.show()

and this is the output:

this would be the equivalent of matlab's impulse() in python


Bilateral Laplace transform and existence of Fourier transform

The bilateral Laplace transform converges in a vertical strip $a<\text{Re}\{p\}<b$, called the region of convergence (ROC). Compare this to the bilateral $\mathcal{Z}$-transform which converges in an annulus centered at the origin of the complex plane: $r_1<|z|<r_2$. For causal signals we have $b=\infty$ and $r_2=\infty$.

If the vertical strip $a<\text{Re}\{p\}<b$ does not include the imaginary axis, i.e., if $0<a<b$ or $a<b<0$, the bilateral Laplace transform does not converge for $p=j\omega$ because the imaginary axis is not inside the ROC. Consequently, the Fourier transform does not exist because the corresponding integral does not converge. For the $\mathcal{Z}$-transform the analogous case would be that the ROC does not include the unit circle, in which case the discrete time Fourier transform does not exist.


What is the meaning of $Ta_k$ of fourier series or transform?

The idea is that a Fourier series is only defined for periodic signals. In the discussion in the linked slides, the author is considering a rectangular pulse train with period $T$. That is, a pulse of width $2S$ repeats periodically with a spacing of $T$ between them. The pulses are therefore centered at:

$$
[\ldots, -2T, T, 0, T, 2T, \ldots ]
$$

Now, consider what happens as $T \to \infty$: in the limit, the only pulse that remains is the one centered at zero; the others are infinitely far away. When the author makes the claim that:

$$
\lim_{T\to \infty} T a_k = E(\omega)
$$

He or she is trying to show that the Fourier transform, which is defined for suitably well-formed aperiodic signals, can be thought of as the Fourier series of that signal (which typically wouldn't be defined since the signal is not periodic) in the limiting case of an infinite period. Stated a little differently, you can in some way think of an aperiodic signal as a periodic signal with infinite period.

The multiplication by $T$ in the limit is to account for the differences in definition between the Fourier series and Fourier transform: the series representation typically has a factor of $\frac{1}{T}$, while the transform does not. I don't know that there is a lot of insight to be gained via this analysis, but it shows that the series and transform representations are intimately related.


How does one know how many times to downsample in Harmonic Product Spectrum?

Take a look at this paper on page 4. It has a great graphical representation of the algorithm. I'll try to write up some pseudocode here that you can adapt to your language (looks like Java to me) as necessary.

/* Harmonic Product spectrum PSEUDOCODE
 * @param fftData Discrete Fourier transform of data
 * @param N Number of times we downsample the spectrum to get HPS
 */
int HPS( Complex[] fftData, int N )
{
    // Find magnitude of the FFT
    Real fullSpectrum[] = absOfComplex(fftData);

    // Keep only the positive frequencies (DC to Nyquist)
    Real spectrum[] = discardNegativeFrequencies(fullSpectrum);

    // Make a new array to store HPS
    Real hps[] = copyOf(spectrum);

    // Perfrom HPS:
    // Go through each downsampling factor
    for (int downsamplingFactor = 1; downsamplingFactor <= N; downsamplingFactor++)
    {
        // Go through samples of the downsampled signal and compute HPS at this iteration
        for(int idx = 0; idx < spectrum.length()/downsamplingFactor; idx++)
        {
            hps[idx] *= spectrum[idx * downsamplingFactor];
        }
    }

    return findIndexOfMax(hps);
}

This isn't bug-proof but it should get you started.


Oversampled Binary Image Sensor vs Active Pixel Sensor(CMOS)

The success of your camera depends on too many factors! and they're out of the scope of this page. If the dynamic range requirements of your application are not too demanding then you could take advantage of the fast sampling rate and try to reconstruct a low quality image (by the same principle by which sigma delta modulation works, eg reconstructing a PCM signal from the output of a 1 bit output digital microphone). There's technique called halftoning, which has been the basis for reproducing images in print media since the invention of the newspaper. It is based on sampling an image with 1 bit and diffusing the error. Determining if using a Binary Sensor instead of a CMOS is a good or bad idea depends on your intended application. You should also pay attention to the size of the objects you want to picture. CMOS technology is widely used and you could probably find a cheap sensor to use in your camera.


How do you choose a size for a windowing function in FFT?

Depending on the pitch source (voice? Instrument?), the magnitude peak from a windowed FFT may not even be good for pitch estimation.  If you are post processing FFTs as a sub-component of a more complicated pitch estimation method, then the type of post processing (parabolic log magnitude interpolation, complex Sinc interpolation, phase vocoder estimation, cepstrum, etc.) may alter the choice of more optimal windows.

A shorter FFT window provides better temporal localization for changing pitches.  A longer FFT window provides more frequency resolution, providing the spectrum stays stationary for the duration of the window width.


Signal-to-interference-plus-noise ratio (SINR) is by definition

$SINR = \frac{P}{I+N}$,

where $P$ - is useful signal power, $I$ - interference signal power and $N$ - noise power in the band of interest. Look

http://en.wikipedia.org/wiki/Signal-to-interference-plus-noise_ratio

for example.

So the first formula you've presented seems to be rigth one. Time span of calculation $SINR$ depends on your needs and interference nature. The common way is using duration of $x$. In general you have to distinguish both $x$ and $i$ at first, because $i$ is non-stationary random process so you can't think it as receiver noise. So in most common case in real life you can estimate $SINR$ only if you are able to combat such an interference in you system.

If inteference appears as short pulses in time you can do measurement with relatively short sliding window to estimate $SINR$ in dynamics. If signal power is quite stationary in short averaging you can use previously estimated $P$ in formula above to find out $SINR$ even if demodulator fails to defeat it.

You can combine measurement of $SINR$ at the input of demodulator, $MSE$ at its output and $BER$ at the decoder's output to make some sophisticated modem quality control.


The formula you gave is an estimator for the "statistical power", i.e. the expectation of the square of a continuous  underlying random process:
$$ P_X = E(X^2(t)) $$
Computing the latter theoretically would require to know the probabilistic behavior of $x(t)$, which you don't have access to in real life, especially since:

you get sampled data (not continuous),
you generally have a limited number of realizations.

This is even more complicated when you add variables, because you should know their paired probability behavior. A classical result is that, for variables $X$ and $Y$,
$$\mathrm{Var}(aX+bY) = a^2\mathrm{Var}(X)+ b^2\mathrm{Var}(Y)+ 2.a.b.\mathrm{Cov}(X,Y)$$
This quantity can take different values, depending on the last term: $\mathrm{Cov}(X,Y)$. For uncorrelated variables, the covariance vanishes, and the variables add as:
$$\mathrm{Var}(aX+bY) = a^2\mathrm{Var}(X)+ b^2\mathrm{Var}(Y)$$
So, for uncorrelated signals, powers add. And this generalizes to several signals.
Note that when  you sum over a time frame to estimate such quantities, you somehow assumes that you can estimate quantities by averaging in time, with assume that the signals are correctly sampled for a stationary and ergodic process.
For $K$ signals of  identical power $P$, the combined power can take value:

minimum: $0$ when signals sum to zero,
$KP$ when uncorrelated,
maximum: $K^2 P$ when signals are identical.

Additional pointers:

Addition of signal variance and noise variance
Variance of Uncorrelated Variables
Signal averaging
Superposing Signals


Orthogonal signal generator using integer arithmetic

I managed to find an answer to my problem.

The problem is with rounding in integer division. For example, -4/3 will be rounded to -1, just as -5/3. Because of this, the integration error is constantly accumulated. Instead of explaining, here I give a code how to fix this.

Code declaration.

// Global variables
int w = 643398L; // (2*PI)*50*2048
int i1u1=0, i1u2=0, i1u3=0;
int i2u1=0, i2u2=0, i2u3=0;
short i1y1=0, i2y1=0;

// Numerical integrator implementation
short NumInt3rd(int u1, int u2, int u3, short y1) {
    int iu = 23*u1-16*u2+5*u3;
    int iy = 240000L*y1;
    short y0 = (((iu+iy)>>12)*2237+65536)>>17;
    return y0;
}

Main function.

// Get voltage measurements (-2048 to +2048)
short vg = (short) Input(0);

// Numerical integrators, 3rd order
short vg_x = NumInt3rd(i1u1,i1u2,i1u3,i1y1);
short vg_y = NumInt3rd(i2u1,i2u2,i2u3,i2y1);

// Downsample frequency to prevent overflow
short wb = w>>5;

// Update integrator #1 states
i1u3 = i1u2;
i1u2 = i1u1;
i1u1 = ((int)((vg-vg_x)-vg_y)*wb)>>6;
i1y1 = vg_x;

// Update integrator #2 states
i2u3 = i2u2;
i2u2 = i2u1;
i2u1 = ((int)vg_x*wb)>>6;
i2y1 = vg_y;

It should be noted that instead of using integer division, which is very expensive in terms of required number of instruction cycles, I rather use bit shift operation. For example, 1/240000 can be well approximated as 2237/2^29 - the approximation accuracy is 0.00169277%.

And here is the sine wave estimation error. As you can see, the results are much better now.


Estimating variance in arbitrary, periodic signal

Do a cross correlation with your known $s$; you will get peaks of height $\|s\|A_k$ at the positions of your shifts.

You just sum up these peaks, call that sum $P$, calculate the overall energy in your signal (i.e. $\|x\|$), and calculate $\text{SNR}=\frac P{\|x\|-P}$ if you need the SNR. If you just need the noise: $\|x\|-P$ (quite directly). Done!


Rotated object detection by Viola jones Algorithm

You are correct, Haar features, as well as LBP and HOG are not rotation invariant. I have experimented with vision.CascadeObjectDetector in MATLAB's Computer Vision System Toolbox, and found that the face detector model it comes with can tolerate about 15 degrees of in-plane rotation.

So if you want to handle in-plane rotation, the easiest thing is to rotate the image by 15 degree increments and run the detector each time.  If you are talking about out-of-plane rotation, then you would need multiple detectors, such as the frontal and profile detectors that come with MATLAB.

There is another possibility. If you need to detect the faces in video, then you can handle rotations and even partial occlusions by tracking them using KLT.  See this example.


Will using window in Welch's method introduce artificial peak in estimated PSD?

It is not any artifact. When you calculate the DFT of your signal you will get following magnitude spectra:

So you can clearly see that around $25 \mathtt{Hz}$, and $50 \mathtt{Hz}$ there is some harmonic content present.

Reason why you don't see your harmonics on PSD in first case is that you are not using any windowing. Probably you are aware of side-lobes attenuation in characteristic of the window functions. For rectangular window it is approximately $13.3 \mathtt{dB} $ and for Hanning you get $31.5 \mathtt{dB} $. It is therefore possible that sidelobes will cover some harmonics of your signal.

Please take a look on your figures; for rectangular window lowest value of noise is approx. $-14$, whereas for Hanning you are getting around $-18$.

Below you can see frequency domain representation of these two windows.

For more info about windows, please refer to:

  F. J. Harris - On the Use of Windows for Harmonic Analysis with the DFT


Transfer function estimation using system identification

Let us denote the output feedback SISO plant transfer function $\frac{Y(s)}{U(s)}=P(s)$ so that the open-loop (OL) transfer function of the system loop is $CP$ where $C(s)=K_p + K_d s + \frac{K_i}{s}$ and the closed-loop (CL) transfer function is $\frac{CP}{1+CP}$. From the attempt at the solution to the problem in the OP, we can see that a second order model might be considered a good idea.
Since the controller transfer function ($K_p, K_d s, {K_i}$) and the system is using direct output feedback, the control input (or actions) may be obtained from the provided data as $u(t) = K_p y(t) + K_d \dot{y}(t) + K_i \int_0^t y(t) \; dt$. The data available is then $\{y(t_i), u(t_i)\}$ which can be directly fed into a system identification tool. Alternately one may used least squares to obtain the second order model parameters ($A, B)$ to fit the state transition model ${\Huge[} \begin{matrix} y \\ \dot{y}\end{matrix}{\Huge]} = A {\Huge[} \begin{matrix} y \\ \dot{y}\end{matrix}{\Huge]} + B {\Huge[} \begin{matrix} 0 \\ u\end{matrix}{\Huge]}$. Care must be taken not to use the CL data to fit the model instead since, as clearly seen from the model, there is no integral term in the state space representation the model is attempting to use. Finally, from a preliminary observation of the graphical representation of the data, the data can be divided into three distinct data set split by time-intervals in which the system undergoes large control actuations which seem to become saturated. These three data sets are observed to have the typical second order behavior reflecting a first or second order system controlled by a PID controller.


In statistical case OFDM (without specific PAPR reduction precoding and with most of subcarriers used) signal could be modelled as white Gaussian process. So PDF of OFDM signal is well-known one for normal distribution:

$f(x) = 1/(\sigma\cdot\sqrt{2\cdot\pi}) \cdot \exp(-(x-\mu)^2/(2\sigma^2))$

For OFDM $\mu=0$. When we calculate crest factor we need to define the number of samples to be used for measurement. We say measurement for $1e+6$ samples is valid since crest factor doesn't increase significantly, when we take more values. In another words if we have e.g. $1e+2...1e+4$ samples, crest factor will vary from one measurement to another. But if we have a million of samples or more, crest factor of white noise or OFDM will converge to one specific value, it is about 13.5 dB. You can check it easily with MATLAB.

It can be derived analitically with aid of PDF of normal distribution. Look, we need to estimate how much peak value exceeds $\sigma$ to estimate crest factor. If for example $1e+6$ samples are used the probability of peak value is $1e-6$. It is rought approximation. For normal distribution large peaks have low probabilty and it decreases with peak amplitude increase (see plot of normal distr. PDF). We need to estimate maximum peak for million samples, so it is one peak for million and its probabilty is $1e-6$.

The equation to be solved is

$P(x)=\int_x^\infty f(x)dx = 1e-6$,

where $f(x)$ is presented above.

To avoid infinite integration we can roughtly simplify (forgive me, mathematicians!) this equation to something like this

$f(x) \le P(x) = 1e-6$.

By usign $f(x)$ instead of $P(x)$ and solving this equation for different $\sigma$ we can find out that $4.5\sigma \le x \le 5\sigma$ in general. The creast factor could be find as $20\cdot\log_{10}(x/\sigma)\approx 13...14$. The result is very similar to what we can observe in practice. It can be more accurate if someone suggest the way to approximate $\int_x^\infty f(x)dx$ term in the original equation (I believe it could be done since $f(x)$ is decreasing in the interval from $x$ to $\infty$).


Z Transform - Do i always need 2 poles for every "peak"?

If your system has a real impulse response, then pole & zeros will be in complex-conjugate pairs; so a pole at $\omega = \pi/3$ implies another pole at $\omega = -\pi/3$ .

Further, by the definition of poles, as the frequency $\omega$ get closer to the pole frequency, the frequency response magnitude increases (so as to say peaks!). So for every peak in the (positive frequency axis) then you would have at least two poles at the same frequency (one for positive frequency and one for the negative frequency). Note that there could be more, but at least two is required for a real system.

However note that, having $2N$ poles does not necessarily mean that the frequency response curve has $N$ peaks (sharp rising transitions). Indeed counter-examples are easiy to find, such as an Butterworth or elliptic IIR lowpass filters having no peaks, compared to a sharp isolated rising transition.

Note that in those filters the poles are bundled together such that their indivual peaks have overlapping effects on the total frequency response (yielding a flat passband rather than isolated peaks).


Discussion - Integrating Accelerometer data to Position data from Frequency domain FFT

The point is that acceleration is the second derivative of displacement. So to get the displacement, one should integrate acceleration twice. To do this, one way is to numerically integrate. A more interesting approach is to use the time-integration property of the Fourier transform:
$$\mathcal{F}\left(\int_{-\infty}^t x(\tau)d\tau\right)=\frac{X(j\omega)}{j\omega}$$
Note that we require $X(0)=0$. So division by $-\omega^2$ means twice integration in time.
Assuming the sampling frequency is fs and a is accel. data, you can use something similar to the following code:

 l = length(a); % assume l even
 w = 2*pi*fs/l*(-l/2:l/2-1);
 A = fftshift(fft(a));
 D = A./-w.^2;
 d = ifft(ifftshift(D),'symmetric');


Difference Between `fft(signal, nfft)` and `fft(signal)`

Assuming a signal x[n] of length L, The Matlab function fft(x) computes an L-point DFT of the signal, whereas fft(x,N) computes an N point DFT of x[n]. For certain reasons N must be greater or equal to L, otherwise a wrong (aliased) DFT result would be computed, so that one cannot reconstruct x[n] back from such a DFT via inverse DFT.

The two calls would result in two DFTs which are of different lengths and which ,in general, take different values. This is because they are in fact L and N samples of DTFT $X(e^{j\omega})$ of the signal x[n].


Can you predict how long a discontinuity will take to reach its target if it is low pass filtered at a certain frequency?

Infinitely long. It's an exponential decay, so it never goes back to 0.

The time constant for a 200 Hz filter is $\frac{1}{2 \pi f} \approx 0.8ms $. About every 0.8ms it drops by about $2/3$ or every 1.84 ms it drops by a factor of 10.

If you define "almost zero" as one millionths you need to wait $6*1.8ms$ or about 11ms.


You made a minimum-phase filter but with a different magnitude response than the original linear phase filter. What you have to do to keep the magnitude the same is to reflect the zeros outside the unit circle into the circle. Since for linear phase filters the zeros are symmetrical with respect to the unit circle (or on the circle), you simply have to double each zero inside the unit circle and leave the zeros on the circle unchanged. Apart from doubling each zero inside the circle you also need to apply a scale factor $1/z_i$ for each zero that gets reflected (where $z_i$ is the location of the zero).

So each pair of zeros of the linear phase filter (assuming $|z_i|<1$)

$$(1-z_iz^{-1})(1-z^{-1}/z_i^*)$$

needs to be replaced by the following scaled double zero:

$$\frac{1}{|z_i|}(1-z_iz^{-1})^2$$

This will keep the magnitude unchanged while transforming the linear phase function to a minimum-phase function.


Help with equations for exponential ADSR envelope

I think that what confuses you is that a decreasing exponential ($e^{-x}$) does never reach 0, so an ADSR generator with truly exponential segments would stay stuck ; because it would never reach the target value. For example, if the generator is at the height of the attack phase (say $y = 1$) and has to land to a sustain value at $y = 0.5$, it can't go there with a true exponential, because the true exponential won't decay to 0.5, it'll only asymptotically go to 0.5!

If you look at an analog envelope generator (for example the 7555 based circuit everybody seems to use), you can see that during the attack phase, when the capacitor is charging, it is "aiming higher" than the threshold used to indicate the end of the attack phase. On a (7)555 based circuit powered by +15V, During the attack stage, the capacitor is charged with a +15V step, but the attack stage ends when a threshold of +10V has been reached. This is a design choice, though 2/3 is the "magic number" found in many classic envelope generator, and this might be the one musicians are familiar with.

Thus, the functions you might want to deal with are not exponentials, but shifted/truncated/scales versions of it, and you'll have to make some choices as to how "squashed" you want them to be.

I am curious anyway as to why you are trying to get such formulas - maybe it's because of the limits of the tool you are using for synthesis ; but if you are trying to implement those using a general purpose programming language (C, java, python) with some code running for each sample of the envelope, and a notion of "state", read on... Because it's always easier to express things as "such segment will go from whatever value it has just reached to 0".

My two pieces of advice on implementing envelopes.

The first one is not to try to scale all the slopes/increments so that the envelope exactly reach start and end values. For example you want an envelope that goes from 0.8 to 0.2 in 2 seconds, so you might be tempted to compute an increment of -0.3 / second. Don't do that. Instead, break it down into two steps: getting a ramp that goes from 0 to 1.0 in 2 seconds ; and then applying a linear transform that maps 0 to 0.8 and 1.0 to 0.2. There are two advantages to work this way - the first is that it simplifies any computation you'll have relative to envelope times to a ramp from 0 to 1 ; the second is that if you change the envelope parameters (increments and start/end times) midway everything will remain well-behaved. Good if you're working on a synth, since people will ask to have envelope time parameters as modulation destinations.

The second is to use pre-computed lookup table with envelope shapes. It is computationally lighter, it takes out many dirty details (for example you don't have to bother with an exponential not reaching 0 exactly - truncate it at your whim and rescale it so that it is mapped to [0, 1]), and it's dead easy to provide an option to alter envelope shapes, for each stage.

Here is the pseudo-code for the approach I describe.

render:
  counter += increment[stage]
  if counter > 1.0:
    stage = stage + 1
    start_value = value
    counter = 0
  position = interpolated_lookup(envelope_shape[stage], counter)
  value = start_value + (target_level[stage] - start_value) * position

trigger(state):
  if state = ON:
    stage = ATTACK
    value = 0  # for mono-style envelopes that are reset to 0 on new notes
    counter = 0
  else:
    counter = 0
    stage = RELEASE

initialization:
  target_level[ATTACK] = 1.0
  target_level[RELEASE] = 0.0
  target_level[END_OF_RELEASE] = 0.0
  increment[SUSTAIN] = 0.0
  increment[END_OF_RELEASE] = 0.0

configuration:
  increment[ATTACK] = ...
  increment[DECAY] = ...
  target_level[DECAY] = target_level[SUSTAIN] = ...
  increment[RELEASE] = ...
  envelope_shape[ATTACK] = lookup_table_exponential
  envelope_shape[DECAY] = lookup_table_exponential
  envelope_shape[RELEASE] = lookup_table_exponential


For option 2, you could describe the system as a cascade of your FIR filter, a moving average filter and a decimator (taking one of every 12 moving average filter outputs). The FIR filter and the moving average filter could be combined together into an equivalent filter whose impulse response is the convolution of the FIR filter impulse response and that of the moving average filter. So option 1 and option 2 are quite similar in that they both correspond to some filter followed by a decimation by 12.

If you've designed and optimized your FIR filter to meet your specifications (including choices over various parameters including noise performance, resource usage, bandwidth, etc.) then the extra moving average filter would just go against that optimization (for example by reducing the effective bandwith, increasing resource usage, etc.) So from that perspective option 2 would be at best wasteful (or even harmful with respect to meeting your specifications) and option 1 would be the better option.

That said, given the same FIR filter in both option 1 and 2, and provided that that filter's design has enough margin against the specifications to accommodate the cascade with the extra moving average filter without causing the cascade to go out of specifications, the extra filtering is likely to remove more noise (so in that case option 2 would be the better option, even though that isn't exactly a fair comparison from a resource usage perspective).


Assuming that the sensors share the same characteristics, have the same timing (acceleration signals are aligned), the model with $y_1= x + n_1$ and $y_2= x + n_2$, $n_1$ and $n_2$ being uncorrelated noises of the same power, averaging them is a way to reduce the noise. The theory that asymptotically, averaging $N$ sensors reduce the variance by a factor of $N$, hence the noise amplitude by $\sqrt{N}$. For small $N=2$, and with quantization, one should not expect too much (since the asymptotic formula is not valid), but this can be tried, you may have nice surprise. If noise powers are different, this method can be improved to take this into account.

Further, this problem can be addressed as a form of source separation, for which many methods have been developed (independent component analysis, Non-negative matrix factorization, low-rank decomposition to name a few). The choice will depend on how much more information you have, and your algorithmic constraints.


How to use matlab toolbox for 2qam dfe equalizer?

From the MATLAB documentation (https://uk.mathworks.com/help/comm/ref/comm.decisionfeedbackequalizer-system-object.html): "The equalizer uses the reference tap location to track the main energy of the channel". The equalizer object has a property called the reference tap which has default value of one. Looking at your channel coefficients, it looks like the tap with the main energy is the second tap so you should set eq.ReferenceTap = 2.

You can also work on choosing an appropriate step size. The equalizer returns the error signal back, you call it e, you can plot that and use it as a guide to pick a step size. You can also use the in-built function maxstep which gives you the maximum step size needed for a LMS-like algorithm to converge.

Last suggestion I have is unrelated but when calling awgn and the input is real, then the function does not add complex noise so your noise will only appear in the in-phase/real axis. You can get around it just by creating your own complex noise, or by specifying the input as complex with imaginary part equal to zero: r = awgn(complex(filter(h, 1, x), zeros(size(x))), SNR_dB, 'measured').

Edit

The OP is still having issues so I will add the settings that I use to get results.

eq = comm.DecisionFeedbackEqualizer;
eq.Constellation = qammod(0:MQAM-1, MQAM);
[~, eq.ReferenceTap] = max(h);
eq.StepSize = maxstep(eq, r)/2; % maxstep returns largest so choose something smaller to be safe

And here are the resulting error signal and before/after constellation plot.


how to extract telephone audio only from an mp3

There might be much discrepency in the signal's characteristics of the interviewee that depends partialy on the phones they're using, their mobile operator, etc...
It is much easier to extract the radio's speaker signal. Then, everything that does not come form the radio's host must be from one of the interviewees.
What you need to do is to analyse the signal of the radio's speaker that might contains more prominent low-pitched components for example (depends obviously on the station you recorded from). Also use a window of analysis sufficiently large to take into account that every stream must last at least 500/1000ms. That might introduce a little stability in your extraction.
If you want to split the different interviewees, you can pose the simple condition that interviewee/host must alternate.


Actually, yes. As Hilmar pointed out, it might be cheaper to use an ordinary filter, but that depends on the circumstances.

If you have two signals $A$ and $B$ with Fourier transform being $\mathcal{F}$, their convolution $A\ast B$ will be

\begin{equation}
A\ast B = \tilde{\mathcal{F}}(\mathcal{F}(A)\cdot\mathcal{F}(B))
\end{equation}

With $\tilde{\mathcal{F}}$ being the inverse transform. This is the well known "convolution lemma". There is not much of a theory to see. In the Fourier domain, you just replace convolution by (point-wise) multiplication. Your question already contained the answer, to me it just appears that you were not sure about the details.


In Kalman filters why is it necessary to transform the systems dynamics matrix to the state transition matrix?

As I said in the comments, there should be no need at all to do such a transformation if all of your equations start in discrete-time.

The equation
$$
 (s \mathbf{I} - \mathbf{F} )^{-1}
$$
is just the Laplace transform of the solution for the state equation:
$$
\dot{\vec{x_t}} = \mathbf{F} \vec{x_t} + \vec{w_t}
$$

However, I really can't see how $\mathbf{F}$ works in:
$$
\vec{x_t} = \mathbf{F_t} \vec{x_{t-1}} + \mathbf{B_t} \vec{u_t} + \vec{w_t}
$$
because this is a discrete-time equation, and the equation you're looking at generating the sinewave is continuous-time.

Can you clarify?


How to compute Discrete Cosine Transform

You don't need to split image into blocks. The DCT equation can be applied to the whole image. The block division has been chosen for JPEG standard partly because DCT was costly to compute in the past (but that's not the only reason).

You can choose any size of block (including the single block, which is the image itself), then split image into the blocks and apply DCT for every block separately.

After applying DCT on, say 8x8 block (matrix), you get another 8x8 block, but with DCT coefficients (not pixel values). You normalize them, quantize them and then encode. The coefficients in the bottom right corner of the block (I think) can be stripped off to achieve compression (the article about JPEG describes this well).

It is very similar to DFT, which uses sine and cosine bases. The DCT uses cosine bases only. The DC component is located in the top-left corner of transformed image/block and in DFT the DC component is in the middle of the image/block.


Adding Poisson noise to a time series signal

If your data $\mathbf{x}$ has variance $\text{var}(\mathbf{x})$ and you want to add Poisson noise to it with signal to noise ratio of $\text{SNR}$ dB, then you can just generate Poisson random numbers with appropriate variance. The variance of the noise $\mathbf{n}$ you need is $\text{var}(\mathbf{n})=\frac{\text{var}(\mathbf{x})}{10^{\text{SNR}/10}}$. So now you just need to use your Poisson random number generator and then add the noise to get the noisy data: $\mathbf{x}+\mathbf{n}$


Because the discrete Fourier transform is discrete in space, it's periodic in frequency

Let's think of the Fourier Transform of a continuous, periodic function: FT explains the function as a projection of the signal onto a sum of sinusoidal basis functions, i.e. correlates the signal with a bunch of complex sinusoids, each with its own frequency. Thus, any data point contributes to a sinusoid with a particular frequency - data is encoded into the frequency. Thus FT becomes discrete. Now inverse also holds: If the space is discrete, one would require a periodic FT to represent that (non periodic signal on one domain causes continuous signal on the other).

I really like this animation showing the continuous FT.

Mathematically, it is always possible to show that for DFT, $X_K = X_{K+N}$ where $N$ is the period such as in here.

And because the spatial data is real, the frequency data is symmetric. The middle of the square correspond to the Nyquist frequency (half the sampling frequency).

First let's address the symmetry:

I will greatly cite this stackexchange post and for completeness will include some shortened portion here:

Due to its definition, FT will generate symmetric responses for real data. Take the transform equation:
$H(f) = \int h(t)e^{-j2\pi ft}dt$.

Using Euler's formula $e^{jx} = \cos(x) + j*\sin(x)$, we could modify the Fourier transform as follows:
$$
H(f) = \int h(t)(\cos(2\pi ft) - j\sin(2\pi ft))dt
$$

FT for negative frequencies become:
$$
H(-f) = \int h(t)(\cos(2\pi (-f)t) - j\sin(2\pi (-f)t))dt \\
= \int h(t)(\cos(2\pi ft) + j\sin(2\pi ft))dt
$$

Comparing the negative frequency version with the positive frequency version shows that the cosine is the same while the sine has phase inversion of $90^\circ$. Magnitude is the same. Due to this orthogonality, they will both respond to real signals in the same way.

Now for the Nyquist part:

The Nyquist frequency represents the highest frequency that can be represented by the sampling rate (differs from Nyquist rate). The DFT of a real series, i.e. zero imaginary part, results in a symmetric series about the Nyquist frequency - this is where the signal folds. In the image you have given, the mid-part is exactly about where the symmetry is created. This is also where the DC component lies.

I hope these bring in some clarification.


Using the standard deviation of time lags is not a bad idea - the problem is that for very noisy signals such as consonants you won't really get a pattern with peaks. Your suggestion would be more useful in the context of musical instruments sound (for example to measure the inharmonicity of a sound, from violin to piano to bell...)

You can look at the ratio between the value taken at the highest peak in the vocal range (say 80-400 Hz), and the value of the ACF at 0 - this is a common criterion used for voiced/unvoiced classification. Another simple metric that works relatively well with clean recordings is the zero-crossing rate.


Frequency of a periodic signal with distortions

The signal looks like narrow band noise and may not be periodic at all.  The FFT will generally reveal a periodic signal that is in the presence of some noise.

In the FFT results, if you have a periodic signal, you should see peaks at the fundamental and its harmonics.  If the "periodic" signal is unstable (modulated in some way) the peaks will become less distinct.

You can try different window functions on your data to see if more defined harmonic peaks become visible in the FFT data.  For example compare results using a number of windows such as Gaussian and Blackman-Harris windows.

Also the choice of the window size (length of data) used when computing the FFT will affect your frequency resolution (bin spacing).  If you can roughly estimate the period of the signal, choose a data length that is a multiple of the estimated period so the bins line up with the signal harmonics.  In general a longer data set will provide better frequency resolution than a shorter data set, so if you don't know the period of the signal, go with a longer data set.


As hinted out, you firstly shift signal by π by multiplying it by cos(πn). Now, your low frequency components are around π. In the next step of diagram, you filter the components around frequency 0 since the filter in diagram is lowpass. But since your signal was shifted by π, you filter the high frequency components of the original signal. In the next step, you shift your original signal back to 0 frequency.

So, overall you filter the high frequency components.


How to make bode plot when output signal changes amplitude?

You are probably still seeing other frequency components from the transient of the system. I assume you that you simulated your system for a fixed number of periods of the excitation signal. If this is the case then it would also explain why you only see the transient at high frequencies. This would be because a fixed number of high frequency periods take less time than the same number of low frequency periods. Therefore there is not enough time for the transient to die out for the high frequency excitation.


How do I know what is causing the noise in my signal

Since it's home work, I can't give the full solution.

The noise seem to have been artificially added. This type of noise doesn't happen normally
The noise is basically two sine waves that really don't belong there
The best filter type to eliminate these are notch filters.


The point is that

  The output of an LTI system to an input is a weighted sum of shifted versions of the impulse response.

You can see that the input is $$x=(1)\delta[n]+(-1)\delta[n-1]+(1)\delta[n-2]$$
The response to $A\delta[n-n_0],\, \forall n_0$, is $Ah[n-n_0]$. By the principle of superposition, the output to a sum of shifted $\delta$ terms is a sum of shifted $h$ terms. Here:

$$y[n]=(1)h[n]+(-1)h[n-1]+(1)h[n-2]$$
where $h[n]$ is $$ \cdots \quad 0 \quad 1 \quad -2 \quad 3 \quad -4 \quad 2 \quad0\quad\cdots$$


Cascading first order digital filters in C++

My late night brain made foolish mistakes, for the record, if anyone needs this, the working code is as follows:
Header:
float state1;
float state2;
float state3;

Implementation:
In Constructor:
state1 = 0;
state2 = 0;
state3 = 0;

Robert's function:
// this processes one sample

float first_order_filter(float input, float pole, float zero, float *state)
{
    float new_state = input + pole*(*state);
    float output = new_state - zero*(*state);
    *state = new_state;
    return output;
}

Then in the loop, for (int i=0; i < numSamples; i++)
float first = first_order_filter(whiteNoise, 0.99572754, 0.98443604, &state1);
float second = first_order_filter(first, 0.94790649, 0.83392334, &state2);
float third = first_order_filter(second, 0.53567505, 0.07568359, &state3);
out1 = third;

Of course this works like a charm, producing Pink Noise and if the poles and zeros are swapped it produces Blue/ Azure Noise.


Is there an algorithm for finding a frequency without DFT or FFT?

FFT is actually not a great way of making a tuner. FFT has inherently a finite frequency resolution and it's not easy to detect very small frequency changes without making the time window extremely long which makes it unwieldy and sluggish.

Better solutions can be based on phase-locked loops, delay-locked loops, auto correlation, zero crossing detection and tracking, max or min detection and tracking and certainly intelligent combination of these methods.

Pre-processing always helps.


Obtaining signal energy modulation for voice activity detection

From my understanding of your question, you want to ask whether one should look for energy(of FFT coefficients) of bands directly, or one should use MFCC's (energy) for  making VAD decisions.

As MFCC based approach is to mimic the human auditory systems response, it will be better to use MFCC coefficients.

if one uses same bands as MEL Filterbank for calculating band-wise energy (after taking FFT), than the results will be similar to using MFCC'S, as MFCC's are just the transformation from FFT coefficients.


representing simple definite integral as convolution in matlab using CONV

A definite integral is a convolution with idial all '1' filter(by the theory)

No. Where does that theory come from ?

If you want to approximate a definite integral with a discrete sequence, you can simply sum up the samples in the integration interval and scale the sum by the sampling interval.

The output of a convolution of two functions is another function, not a single value.


Constellation diagram for QPSK modulation

I think you may be missing some fundamental understanding of QPSK modulation, which was likely the point of the homework exercise so let me hopefully fill in some blanks, from which you can figure out how to code up in your Matlab.

This diagram of a QPSK modulator may be helpful, taken from http://www.masoodkh.com/projects/communications/18-qpsk-modulatordemodulator:

Note the similarity to your implementation, where the data is converted to +/-1, and multiplied by the sine and cosine of the carrier.  The process of doing this takes your waveform that would otherwise be centered at frequency = 0 (DC) which is the baseband waveform and move it to a carrier frequency (where it is easier to move the signal over great distances). With such a linear modulator, everything about that signal at DC is the same at the carrier in terms of the signal's amplitude and phase, we just change the frequency the signal rides on. So we start with a "DC" carrier and move it to $f_c$ which is what the block diagram is doing and what you are doing in your code. The constellation is the baseband amplitude and phase assignments that are made at DC prior to shifting the signal up to a higher frequency carrier.

We call the baseband data IQ data as it is typically a complex data stream with real and imaginary components, "I" is in-phase and represents the real portion and "Q" is the quadrature-phase and represents the imaginary components. This is certainly the case with QPSK, we take every other bit and assign it to the real axis and the others to the imaginary axis to get I+jQ at baseband.

Your original constellation and how to plot it without the noise should now be very straight forward (plot I+jQ before it gets multiplied by the sine and cosine functions).

To plot the constellation after adding noise seems to me would be very challenging to ask for a homework question, which makes me think there is an easier way than what I proceed to describe below (and hopefully someone else can provide that), but here are my thoughts:

It is not clear if the constellation with noise should be plotted for all samples over the symbol duration, or if you are to get an average of the noise to plot the one sample that a demodulator would determine before making a decision (it would be the latter if the assignment was clear that there should be only one sample per symbol), but let me describe both cases:

Note that the "carrier" is multiplying your baseband signal by $e^{j2\pi f_c t}$ since

$$e^{j2\pi f_c t}= \cos(2\pi f_c t)+ j\sin(2\pi f_c t)$$.

We assign a $j$ to the Q input of the modulator and multiplying by this identity according to Euler matches the block diagram above. The implementation is specifically equal to taking the real part as follows:

$$Re[(I+jQ)(e^{-j2\pi f_c t})] = I\cos(2\pi f_c t)+jQ\sin(2\pi f_c t)$$

Importantly based on the 4 possible dibits on the input ({1,1},{1,-1},{-1,1},{-1,-1}) we assign 4 phases (45°, 135°, 225°,315°) to the carrier. Note that the math and seeing the result is actually much easier using the $e$ instead of sines and cosines, understanding that $Ke^{j\theta}$ is a vector on the complex plane with magnitude $K$ and angle $\theta$.

Note that multiplying a signal by $e^{j\omega_c t}$ (where $\omega_c=2\pi f_c$) simply shifts the signal in frequency in one direction, depending on the sign of $\omega_c$:

A modulated signal that is varying in amplitude and phase (and in your case of unfiltered QPSK, varying predominantly in phase only) will have the same amplitude and phase variation at baseband (DC) as at any carrier, as represented in the IQ diagrams below. It's just in the case of a carrier, the amplitude and phase is plotted with reference to the carrier. This plot shows how multiplying a signal at the carrier frequency by $e^{-j2 \pi F_o t}$ (where $F_o$ was used in the diagram but is the same as your $f_c$) will move the signal from the carrier back to baseband (DC).

Ultimately to see your constellation with the noise, you want to move your signal that is at the carrier back to DC. A further complication is that your signal as modulated is completely real (note that I took the real part after multiplying by $e^{j2 \pi f_c}$ when we did the upconversion to make it match your implementation. A complex signal can be single sided (positive frequencies only for example) as we described would happen after multiplying by $e^j2\pi f_c t$ (and note that I specifically multiplied by $e^-j2\pi f_c t$ in order to match your implementation, which means I shifted the frequency to a negative frequency carrier!). Once we take the real part of a complex frequency, what was once a single sided spectrum (if it was, doesn't have to be but will always be the case is the positive spectrum is different from the negative spectrum), will become a double sided spectrum that is complex conjugate symmetric: For all real signals the spectrum will always be complex conjugate symmetric- meaning the negative frequency spectrum will equal the positive frequency spectrum in magnitude and opposite phase.

Since your signal is real, there will be two complex conjugate symmetric spectrums, one at $+f_c$ and one at $-f_c$. (When I took the real part of what was originally just a spectrum at $-f_c$, the companion spectrum at $+f_c$ is immediately created by the operation of taking the real part).

That said, to move your signal back to DC, simply multiply by $e^{-j2\pi f_c t}$ which will shift both spectrums to the left, and you will get one at DC and another that is at $-2 f_c$. It is possible here that the spectrum will be inverted, as I did not follow the math carefully in my processing described, in that case it would be better to multiply by $e^{-j2\pi f_c t}$, or simply swap I and Q and change the sign of Q which will also invert the spectrum. You then need to filter out the higher frequency signal that will be at twice the carrier frequency in order to get your baseband constellation only. This is the part I am surprised you would have to do as a homework problem; and if you had access to the noise samples at baseband this would be much easier but given the way the noise was added it doesn't appear that you do?

If you needed to plot the constellation with noise as 1 sample per symbol (IQ dibit), then you would simply take the average over a symbol period of the signal after converting back to baseband. Actually, by taking the average, the required low pass filtering will also be done (by nature of averaging is a low pass filter). So such a result for the constellation with noise is actually easier achieved when plotting one sample per symbol. Note too that taking such an average, if synchronized with the symbol period, represents the ideal matched filter for your example of QPSK with a rectangular pulse shape. This is often depicted on diagrams of receiver implementations as an integrate and dump filter (moving average where you accumulate or integrate over the symbol period, dump the result and reset the moving average).


On the symmetry of a $2$-dimensional discrete-time signal

What you have displayed is true for which it's odd symmetric for only one of the variables $n1,n2$ and even-symmetric for th other.

For 1D signals, an even-symmetric signal has the property that
$$ h[-n] = h[n] $$

and an odd-symmetric signal has the property as:
$$ h[-n] = - h[n] $$

For 2D signals, like $h[n_1,n_2]$, more symmetries are possible:

For even signals, the signal with the property:
$$ h[-n_1,n_2] = h[n_1,n_2]$$
is even-symmetric only in $n_1$ variable, whereas the signal with the property
$$ h[n_1,-n_2] = h[n_1,-n_2]$$
is even-symmetric only in $n_2$ variable. Consequently the signal withe property that
$$ h[-n_1,-n_2] = h[n_1,-n_2] = h[-n_1,n_2] = h[n_1,n_2]$$
is even-symmetric in both $n_1$ and $n_2$.

For odd signals, the signal with the property:
$$ h[-n_1,n_2] = -h[n_1,n_2]$$
is odd-symmetric only in $n_1$ variable, whereas the signal with the property
$$ h[n_1,-n_2] = -h[n_1,n_2]$$
is even-symmetric in $n_2$ variable. Consequently
$$ h[-n_1,-n_2] = -h[n_1,-n_2] = -h[-n_1,n_2] = h[n_1,n_2]$$
is odd-symmetric in both $n_1$ and $n_2$.

So the property you have displayed
$$ h[-n_1,-n_2] = -h[n_1,n_2]$$

implies an even-symmetry on $n_1$ and odd-symmetry on $n_2$ OR an odd-symmetry on $n_1$ and even-symmetry on $n_2$


Echo cancellation supporting long delays and without frequency domain processing

I believe you logic is reasonable. The fact that there are companies offering algorithms for long echo tail (even above 1 second) with reasonable performance indicates that an algorithm similar to your suggested logic can work.


Why motion cannot be estimated in DCT domain using block match?

Block matching will only work well if the pixel offset is (close to) a multiple of the block size, in both dimensional axises.  Odds against that might be (8*8-1):1 or (16*16-1):1.


Your reasoning is correct. The signal $x(t)=|\sin(2\pi f_0t)|$ is periodic with period $T/2$, where $T=1/f_0$. Consequently, it has components at frequencies $2k/T$, $k=0,1,\ldots$, i.e., at $0$Hz, $2000$Hz, $4000$Hz, etc. Since the ideal low pass filter has a cut-off frequency of 1500Hz, all harmonics are completely suppressed, and only the DC component appears at the output.

I'm sure that you know how to compute the DC component of $x(t)$ (which, if positive, equals its RMS value).

Your concerns about actually implemented filters and capacitors etc. are irrelevant to the question, which explicitly states that the low pass filter  can be assumed to be ideal. But even a non-ideal low pass filter of appropriate order would suppress all harmonics sufficiently, such that they would be buried in noise.


What methods can be used to compute homography

I've never actually done Panorama Stitching, but I was talking to people who were, and I was involved in a project where we had to estimate the homography for fine localization for visual navigation, so here it goes.

Basic approach

Rough outline of the approach using RANSAC (you can read more of the approach in this Master thesis, page 12,13 - unfortunately, in Croatian):

estimate pairs of corresponding points for two images $A$ and $B$ (lot of them) and get rid of the obviously bad correspondences

(detailed summary of a feature matching approach in my Master Thesis, pages 13-20 - In English)
pick any four points for the initial estimate of homography $H$

(most details on the numerical methods used to estimate the homography can be probably be found in my Bachelor Thesis (in Croatian), but there should also be some details in this Master Thesis (Cro) and this one (Eng))
project all points from the image $A$ onto the image $B$ using the first estimation of the homography $H$
compare the projected points from $A$ to the actual points present in $B$ (using Euclidian distance)
discard all the point pairs where the distance between projected point from $A$ and corresponding point from $B$ is greater than a threshold $t$ -- those points are outliers
(this can be repeated several times, possibly with a more constrictive threshold $t$ every time, to get rid of more outliers)
the points that are left are our population points
using only the population points, do the final estimation of homography $H_{final}$ just from those points

Some additional info

Coincidentally, if you have the initial rough matching results, the rest of this is implemented in OpenCV, it's a function called findHomography and there is a good StackOverflow post on that. It includes RANSAC (and from what I can read from that post, Least Squares method should have a similar goal as RANSAC, but RANSAC gives better results).

If you need better estimates, you can find references for methods other than RANSAC in this Master thesis (Eng) - pure RANSAC approach is compared to something called MLESAC, and RANSAC coupled with Least Median of Squares method.

Suggestions on an efficient implementation of Homography estimate from points (presuming that you do not want to use OpenCV or some other available implementation) can be found in my Bachelor Thesis (Cro).

To estimate the initial (not-very-precise) pairs of correspondences, you should use matching methods. Those methods are based on:

detecting feature points

(there's a lot about that in this DSP question, and especially in this answer you will find a link to the current widely used framework for comparing feature detectors)
using appropriate feature descriptors

(nice info in this DSP question)
doing matching between those points and possibly some kind of filtering to further better the matches obtained

Matching a point from an image $A$, $x_{A}$ with a descriptor $d(x_{A})$ is usually done just by finding the point in the image $B$, $x_{B}$ with the most similar descriptor $d(x_{B})$ according to Euclidian distance.

The most simple (and surprisingly efficient) filtering method is to compare the distance from $d(x_{A})$ to the closest descriptor $d(x_{B_1})$ and second closest descriptor $d(x_{B_2})$. If $\frac{D_{eucl}(d(x_{A}), d(x_{B_1}))}{D_{eucl}(d(x_{A}),d(x_{B_2})}$ is greater than some coefficient $k \in [0,1]$ (typically $k = 0.6$), the match is rejected. First introduced by D.G. Lowe: Distinctive image features from scale-invariant keypoints.

Additional filtering can be done according to spatial arrangement of the corresponding points. One such interesting method was described in J. Sivic, A. Zisserman: Video Google: A Text Retrieval Approach to Object Matching in Videos.


I think you should plot something like:

t = [0:0.05:1];      %20Hz sampling
a = sin(2*pi*2*t);   %2Hz sine wave
b = sin(2*pi*18*t);  %18Hz sine wave
plot(t, a, 'bo');
hold on;
plot(t, b, 'ro');

T = [0:0.001:1];     %1000Hz sampling frequency
A = sin(2*pi*2*T);
B = sin(2*pi*18*T);
% plot for 1000Hz sampling frequency
plot(T, A, 'b');
plot(T, B, 'r');

legend({'2Hz @ 20 Hz f_s', '18Hz @ 20 Hz f_s', '2Hz @ 1000 Hz f_s', '18Hz @ 1000 Hz f_s'})

That allows to understand which exact samples of your sinusoids you are taken at a given sampling frequency:


What is an example of using aliasing to your advantage when recovering an input signal?

you can successfully recover (or reconstruct) the signal so long as $f_2<F_s/2$

Nope. The sampling theorem says that you need 2 samples per Hz of bandwidth, so in this case you'd need

$$Fs > 2*(f_2-f_1)$$

For more info on how this works google "bandpass sampling" or just read https://en.wikipedia.org/wiki/Undersampling.

The basic idea is that sampling creates a periodic repetition of the original spectrum with a period of $F_s$. As long as none of the repeated spectra falls on the same frequencies as the original spectrum, you don't get aliasing.


How can I upsample 22 kHz speech audio recording to 44 kHz, maybe using AI?

What you want to do is called "Super Resolution" in the field of imaging.
This is an ill posed problem.
Hence in order to solve it you need some prior / model about your audio data.
For example you may look at the model in Speech Super Resolution Using Parallel WaveNet.


The Wiener Filter for the Time Derivative of the Signal

The key point here is the derivation / integration is a Linear Operator.
Since the Wiener Filter is also a linear operator it makes things easier.
Write the problem using the Wiener Filter Model:

Neglecting the integrator we know the Wiener Filter for $ \hat{H} \left( w \right) $ (See Wikipedia - Wiener Deconvolution):
$$ G(f) = \frac{1}{\hat{H}(f)} \left[ \frac{ 1 }{ 1 + 1/(|\hat{H}(f)|^2 \mathrm{SNR}(f))} \right] $$
By defining $ H \left( w \right) = \frac{1}{j w} \hat{H} \left( w \right) $ as the fourier operation of integration (This is not exactly accurate, see The Integration Property of the Fourier Transform)
$$ G(f) = \frac{1}{{H}(f)} \left[ \frac{ 1 }{ 1 + 1/(|H(f)|^2 \mathrm{SNR}(f))} \right] $$
Basically, we defined the integration as a filter and used the known building blocks.
The interpretation is, with respect to the non integral case, is apply the derivative operator on the estimator and have some regularization on the SNR estimation as we estimate a derivative.
Resources

How Is the Formula for the Wiener Deconvolution Derived?
Question on Wiener Filtering.
1D Deconvolution with Gaussian Kernel (MATLAB).


Here is a pretty clear paper explaining secrecy rate and uses over fading channels.

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.105.2844&rep=rep1&type=pdf

See equation 3 for the formal definition. One way to think about it is like a tug-of-war between the main channel and wiretapping channel. If the main channel has a higher capacity then the secrecy rate is positive, otherwise it is negative.


How to visualize and understand Gabor filter using MATLAB

The following is my implementation, of the above problem :

In my_gabor_filter.m, I have the following code :

function psi = my_gabor_filter(x,y,mu,nu,sigma)
    phi = pi*mu/8;
    f = sqrt(2);
    k_max = pi/2;
    k_nu = k_max/(f^nu);

    % This is the wave vector
    k_vec = [k_nu*cos(phi),k_nu*sin(phi)]';
    z_vec = [x,y]';

    k_vec_norm = norm(k_vec);
    z_vec_norm = norm(z_vec);

    exp1 = (k_vec_norm/(sigma^2));
    exp2 = (exp(-((k_vec_norm^2)*(z_vec_norm^2)/(2*sigma^2))));
    exp3 = (exp(complex(0,k_vec'*z_vec))-exp(-(sigma^2/2)));

    psi = exp1*exp2*exp3;
end

and, in plot_gabor_filter.m I did the following :

function plot_gabor_filter()
    mu = 0;
    nu = 0;
    sigma = 1;

    [X,Y] = meshgrid(-2:0.1:2,-2:0.1:2);

    Z = zeros(41,41);

    for i = 1:1:41
        for j = 1:1:41
            Z(i,j) = imag(my_gabor_filter(X(i,j),Y(i,j),mu,nu,sigma));
        end
    end

    surf(X,Y,Z,'EdgeColor','black')

    xlabel('x axis');
    ylabel('y axis');
    zlabel('g(x,y)');
end

and when I run plot_gabor_filter, I get the following output, which helped me in visualizing the gabor function :


How does the Magnitude of an FFT relate back to the Time Series wave it was created from?

You can calculate the coherent or DC gain of a window with:

coherent_gain = sum(window)/n;

with a Hann window, and n = 40, this comes to 0.4875, as hotpaw2 mentioned.

You can therefore eliminate the effect of the coherent gain of the window on your spectral estimate by replacing windowed_time_signal = time_signal.*window; with windowed_time_signal = time_signal.*window/coherent_gain;.

If your main priority is very accurate amplitude - and you're unconcerned with spectral leakage as Sebastian mentioned, then the flattop window
is the way to go. Here's the matlab docs.

Just replace window = hann(n)'; with window = flattopwin(n)';.

Also, if you have the ability to collect more samples, that will likely improve the spectral estimate from your FFT (assuming the underlying process in your application is somewhat stationary/ergodic). You can also take FFTs of successive (possibly overlapping) slices and average them.

Furthermore, as hotpaw2 mentioned, selecting a sample rate depends very strongly on your application, and how many samples you intend to collect. That said, I've often used fs = 2.56 * bw, with good results when it comes to spectral analysis or digital filtering.


What happens when you wrap a 3x3 kernel at the corners of an image?

There are three main attitudes when dealing with borders of images and filter kernels.

The simplest is to only process the pixels that have enough neighboring pixels. As a result, some pixels on the sides are left unprocessed, and the filtered image can be cropped.

The shape of the  filter can be adapted (notably in size) to fit the neighborhood of of the running pixel.
The image is extended to the appropriate size. There are many extension methods, from filling with zeroes or constant values to clever extrapolation.

Traditionally, symmetric or antisymmetric extensions are used in many cases, since this is quite consistent with the symmetries or antisymmetries of filter kernels themselves, and this can be performed efficiently. Wraparound is also performed, and sometimes implicit when using FFT filtering or discrete wavelet restoration. Due to discontinuities (the far ends of the image can be quite different), a preprocessing is advised (windowing, smoothing, background removal).


Hough transform is an edge based operation. For this reason, if you could extract the edges of the outer boundary (together with many noisy pixels) and adjust the radius parameter of the Hough transform, you will be finding the circle of the iris, but not anywhere else. You might want to check:

[centers,radii] = imfindcircles(A,radiusRange)

radiusRange is the parameter you should get right.

However, sometimes for some noisy images, Circular Hough Transform might not be the solution. In such cases you might try Hough ellipses or other methods for locating geometric primitives.


There are no doubt plenty of ways you could define one of the measures in the linked question, but the question remains, after doing so, is the measure in any way meaningful?

This all depends on what you intend to do with the metric after computing it: If it's a perceptual estimate (quality of a Codec,  subjective similarity and so on), you will need a set of measured perceptual data to which you then attempt to fit (mathematically) your measure.


Echo cancelling using autocorrelation function

Haven't you already got it there?

$$
\hat{k} = \arg \max_{l > \frac{k}{10}} r_{xx}(l)
$$
and
$$
\hat{\mu} = \sqrt{\frac{r_{xx}(0) }{ r_{ss}(0) } - 1}
$$

The R code below outputs the figure and:

  "k Estimate: 613  vs  613  mu Estimate :  0.747619585689531  vs  0.768768"

R Code

# 26617

T <- 10000
mu <- 0.768768
k <- 613

s <- filter(runif(T,-1,1), rep(10/k,k/10), circular = TRUE)

x <- s + mu*c(rep(0,k),s[1:(T-k)])

par(mfrow=c(2,1))
r_xx = acf(x, lag.max=1000, type="covariance")

min_idx <- floor(k/10)
mx <- which.max(r_xx$acf[(min_idx+1):10000])

khat <- min_idx+mx-1 # R indices start from 1, not 0

points(min_idx+mx,r_xx$acf[khat], col="red", pch=19)

r_ss = acf(s, lag.max=1000, type="covariance")

muhat <- sqrt( r_xx$acf[1] / r_ss$acf[1] - 1 )
print(paste("k Estimate:" , khat, " vs " , k, " mu Estimate : ", muhat, " vs ", mu ))


As far as I understood, by image derivation you mean extracting edges. I would recommend to filter the image by a relatively large Gaussian filter. If computational cost of image derivation is uncritical to your work, I would recommend using canny edge detector. It is less sensitive to noise and does not fool by noise, and finds weak edges along with strong edges. Matlab instruction for that:

   [MinThresh MaxThresh]=[-0.3 0.5];
   EDGE_No_SMOKE=edge(im,'canny',[MinThresh MaxThresh]);

and the result is (I know it might not be the results you're seeking, however playing with Threshold variables and filter size will bring you desirable results):

Note that you do not see smoky effect anymore. Also about those wrong edges, you can remove them using image opening and closing techniques.


Estimating the standard deviation of Gaussian filter from smoothed image

You can exploit the convolution theorem: Convolving with a Gaussian corresponds to Multiplication with a Gaussian of inverse $\sigma$ in frequency. Let x be the input signal and y be the output signal. Further, Let g be the Gaussian filter.

$$
y = g*x
$$

Or, in frequency domain (using captial letters)
$$
Y=GX
$$

Hence, if both Y and X is known, you can do
$$
G=Y/X
$$
to get your Gaussian Kernel. However, due to boundary conditions etc this will not yield the perfect Kernel. You have to fit a Gaussian function to it to get its standard deviation.

Additionally, if you have no knowledge of the input signal, but know its white noise, you can also fit a Gaussian Curve on the FFT of the output (since the input is white noise and hence flat):

N = 16*128;
img = randn(N,N);
s = 2;

G = fspecial('gaussian',[25 25],s);
Ig = imfilter(img,G,'replicate');
subplot(3,2,1);
imagesc(img);colormap(gray); axis off; title('Noise image')
subplot(3,2,2);
imagesc(Ig);colormap(gray); axis off; title('filtered image')

% Perform FFT
Img = fft2(img);
GImg = fft2(Ig);

subplot(3,2,3);
imagesc(abs(Img));
subplot(3,2,4)
imagesc(abs(GImg));

subplot(3,2,5);
fy = -N/2:(N/2-1);
hold off;
plot(fy, cumsum(fftshift(abs(GImg(:,1)))));
hold on;
plot(fy, normcdf(fy, 0, N/(2*pi*s))*sum(abs(GImg(:,1))), 'r');
hold off;

subplot(3,2,6);
div = abs(GImg ./ Img);
plot(fftshift(abs(div(1,:))));
title('Division of known input with output');


It seems that you've correctly found the two poles: $z_{\infty,0}=2$ and $z_{\infty,1}=-\frac13$. The only other thing you need to know is that the possible ROCs are limited by the pole radii, and that there can be no poles inside an ROC.

Note that $z<2$ doesn't make much sense because $z$ is a complex number. If you mean $|z|<2$ then the pole at $z=-\frac13$ would be inside that ROC, so $|z|<2$ is not a valid ROC for the given transfer function.

The three possible ROCs for your example have the following forms:

$|z|<a$
$a<|z|<b$
$|z|>b$

The numbers $a$ and $b$ are of course related to the poles of the transfer function. I trust that you can take it from here.


Proving that the H1 estimator ignores noise on the output:
If we assume that the output is the input processed by the transfer function plus some output noise, we have $Y = H(f)X(f) + N$, where everything is in the freq domain.
So, $S_{xy}(f) = E[X(f)^*Y(f)] = E[X(f)^*(H(f)X(f) + N(f))]$. where E[] denotes expectation and * denotes complex conjugate.
If you multiply this out, you get
$S_{xy}(f) = H(f)E[X(f)^*X(f)] + E[X(f)^*N(f)]  = H(f)S_{xx}(f)$
Crucially, we assume that the input is uncorrelated with the noise so the expectation of the noise term cancels to zero.
Rearranging, we get the H1 estimator, which is not affected by noise on the output:
$H1 = \frac{S_{xy}}{S_{xx}}$
If you try this with noise on the input, the answer isn't as clean - this is because the H1 estimator is not immune to noise on the input. If you repeat a similar proof for the H2 estimator you will be pleasantly surprised :-)


If by feature you mean a group of pixels (like Goofy) I suggest to try SIFT + SVD.
 (http://weitz.de/sift/)


Library to differentiate people by their voice timbre

The task of taking a long contiguous audio recording and splitting it up in chunks in which only one speaker is speaking - without any prior knowledge about the voice characteristics of each speaker - is called "Speaker diarization". You can find links to research code on the wikipedia page.

If you have prior recordings of each voice, and would rather do classification, this is a slightly different problem (Speaker recognition or Speaker identification). Software tools for that are available here (note that general purposes speech recognition packages like Sphinx or HTK are flexible enough to be coaxed into doing that).


Your answer is correct. Each pole gives 90 degrees phase lag. Furthermore, here can you download a file lnyquist1 which is able to give a more accurate result when you have a pole at zero.

num = 1;
den = [1 0 0 0];
lnyquist1(num,den);

In addition, WolframAlpha is also able to give an accurate result, see here.


Like @sansuiso said, compressed sensing is a way of acquiring signals that happens to be efficient if the signals are sparse or compressible.

Compressed Sensing is efficient because signals are multiplexed, hence the number of multiplexed samples (called measurements) is smaller than the number of samples required by Shannon-Nyquist where there are no strong assumptions on the signal.

In the noiseless case, it can be shown that compressive sensing reconstruction solver can recover an exact solution.

In the compressible case, as opposed to the strictly sparse case, it can be shown that the reconstruction error is bounded.

And yes most signals, including ultrasounds are somehow either sparse or compressible. It generally comes down to figuring out the dictionary where the signal is sparse. Domain experts generally know these things.

The interesting question you have is: Imagine you have a non sparse signal and then add zeros to make it sparse and then use compressed sensing to sample that signal, wouldn't it better than directly sampling the full signal?

The answer is no.

It turns out that the sampling requirements for which CS work require more infromation than just performing a full sampling of the original (full/non-zero) signal. In other words, the number of CS measurements required would be higher than the number of non-zero elements in the signals. By sparsifying the signal, you are "losing" on purpose the information about where the signal is supported (i.e. non-zero). The hard part of Compressive Sensing and attendant reconstruction solvers is to find that location where those non zero elements of the signal live: If you know beforehand the locations of those non zero elements, then there is no need to go to a less efficient method of sampling that signal. Indeed, finding the location of the non-zero elements of a signal is the reason we talk about compressive sensing being NP-Hard, BPP and so forth..... Up until 2004, we thought it was hard to do.

Let me put it some other way: Let us assume a signal has K non zero components. If you know the location of these K elements, then you need only K information to know your signal. If you add zeros anywhere in the signal and make that signal of size N, you now need to sample the signal N times through traditional sampling or  O(Klog(K/N)) times with a compressive sensing approach. Since O(Klog(K/N) > K, losing the information about the location of the non zeros elements has yielded a larger set of samples/measurements.

You might be interested in reading my small blog on the subject:
http://nuit-blanche.blogspot.com/search/label/CS
And the following resource:
http://nuit-blanche.blogspot.com/p/teaching-compressed-sensing.html


What you are trying to do has been tried over and over by hundreds of researchers and there is quite a large body of work about this. Check the proceedings of the ISMIR conference. Even if it is not up to date, read Elias Pampalk's thesis : http://www.ofai.at/~elias.pampalk/publications/pampalk06thesis.pdf

To quickly orient you on the right track:

Music can be similar according to many dimensions: a) timbre/texture/genre ; b) rhythmic pattern ; c) melody/chord progression... and more! From your message it's not clear what you want to measure!

If you are interested in a) the features you might want to look at are MFCC (Mel Frequency Cepstrum Coefficients), since they somehow capture the way human hearing works (frequency warping, log scale), since they are decorrelated (making modeling easier), and since they have lower dimensionality (13 coefficients vs 2048).
If you are interested in b), look at a feature called "Fluctuation Patterns" (Pampalk, in short autocorrelation of the signal in the 0.1 .. 10 Hz range over a few bands) ; or Whitman's "Penny" features (FFT of the MFCC along the time axis).
If you are interested in c), look at chromagrams. Start with Ellis' chromagram code (http://labrosa.ee.columbia.edu/matlab/chroma-ansyn/) then move up to Mauch's implementation if you need something more robust (http://isophonics.net/nnls-chroma).

That's for the features. Now you'll have to think of a better way to compare your songs once they have been represented as a sequence of those features. Computing pairwise differences between sequences is not very smart - eg: comparing a song and the same song offset by some silence will yield a difference while it is exactly the same! You'd rather compare the distribution of those features ; for example compute the mean / standard deviation of the features over song A and the mean / standard deviation of the features over song B and then take a probabilistic distance (KL, Bhattacharyya over those).

Last point, but which will matter later: computing the distance between a song and the rest of the corpus to find the nearest matches is quite inefficient. When dealing with large collections, techniques like LSH or Ball trees allow such nearest neighbors queries to be performed without explicit comparison with the whole corpus.

As an aside, tempo detection is an entirely different matter. If you want to look into it, best performance / accessibility paper on the topic is Ellis' Beat Tracking by Dynamic Programming. http://www.ee.columbia.edu/~dpwe/pubs/Ellis07-beattrack.pdf . It's incredibly simple but is close to state of the art algorithms.


How to prove that this system is an invertible system or not?

To invert your system you take the derivative in $t$
$$ \dot y(t) = x(t) - \int_{-\infty}^{t} e^{-(t-\tau)} x(\tau) d \tau$$
Then $$x(t) = \dot y(t) +  \int_{-\infty}^{t} e^{-(t-\tau)} x(\tau) d \tau$$
The problem is that the integral reduces to $y(t)$ as well.
If you look at the response, the transfer function is $1/(s+1)$ the inverse response would be $s + 1$, a system with more zeros than poles. Does it mean it is not invertible in the context you are asking it.


How to calculate the entropy of random normally distributed samples?

It's (fitted) parameters are mean = 35 units, std.dev. = 8 units and I read the values as 8 bit integers.

The choice of parameters hint at a normal distribution. Is this the case? You might want to obtain a histogram and perhaps substitute mean and standard deviation with modal mean and range. This will also help with an estimation of entropy.

  I'm trying to determine a theoretical rate of entropy generation per reading, so for example 3.4 bits /reading or the like.

A first step would be to estimate the Shannon Entropy. In this case, your alphabet is the range of 256 symbols that are possible to be formed with an 8-bit integer and you would be estimating the entropy of a source that generates these symbols with probabilities dictated by the histogram of a long enough sequence of symbols.

If these 256 symbols appeared with equal probabilities in your sequence then you would need all 8-bits to represent them. However, if a subset of symbols appeared with higher probabilities than others then, on average, it would be possible to represent this sequence with less bits per symbol. To determine Shannon's Entropy, obtain the histogram of a long enough sequence, divide each entry by the sum of all entries so that it is now a probability density function (PDF), i.e. it has a total sum of 1.0, and then calculate:

$$H(X) = -\sum_{k=0}^{N} x_i \log_2(x_i)$$

Where $X$ is your PDF (derived above) and $N$ is the length of the histogram.

  PS. As entropy is actually generated by the observer rather than by the process itself, I'm particularly interested in the significance of the bit depth of the individual readings.

Entropy is a characteristic of the source. This particular way of deriving entropy applies to a particular kind of source called ergodic and this is why I am also emphasising the "long enough sequence" (For more information please see this and specifically section 5.)

I am not sure what you mean by "significance of the bit depth". The idea is this: Suppose that you have an alphabet with 256 symbols. Obviously, you need 8-bits to represent them. But, if you were to discover that your source is actually generating numbers in the range $0..32$ for 98% of the time and a single 256 value for the rest 2% of the total length of a given sequence, then, for 98% of the time you would be wasting 3 bits, because $2^{5} = 32$ and you would need all 8 bits for only 2% of the time. Imagine these bits on a sequence:

00000011  00000100  00001110  00001001  00001100  10001101  00011001  00001110

If you could somehow squeeze those 3 bits out of your stream, you could save time of transmission (therefore fit more data through a channel) or storage space.

This is why figures such as 2.5 bits make sense. You can't define half a bit BUT over two symbols (i.e. 3 bits + 3 bits = 6 bits) you would be saving 1 bit (because 2.5 bits + 2.5 bits = 5bits).

For more information on encodings that allow you to effect exactly that "squeezing", please see variable length coding.

Hope this helps.

EDIT:
If you are absolutely sure that the source is normally distributed, then substitute $X$ for an expression of the normal PDF and the sum for an integral over the sufficiently probable figures


What is the difference between normalized peak of correlation, versus peak of correlation divided by its average?

To add some perspective, you could go back to the probabilistic interpretation of correlation.
Remember that the correlation is used because it measures the degree of similarity under some linear generative model of the signal knowing the template (additive gaussian noise). The autocorrelation gives a measure of the log-Probability.

To go back to your question, there are a number of free parameters, notably the variance of noise, and this relates to your choice of normalisation.


Software to extract frequency spectrum data from an audio file

An important thing to consider is that the pitch of musical notes does not directly correspond to frequency peaks in a spectrogram.  This is because the energy of typical musical sounds are partially to mostly in the overtones and harmonics, at frequencies well above the note pitch (depending on the source of the music).  So it would help to investigate algorithms for pitch detection/estimation, and the difference between those and the frequency spectrum of music.


How to identify useful elements of audio recording and ignore lulls?

EDIT: I don't know why i thought this question was matlab based, but my solution is written with matlab in mind. a quick search showed my the python scikit has otsu's method for automatic thresholding (this is what matlab uses). Other than that I think most of the code is relatively python safe and can easily be translated. My apologies for my error

from skimage import filter
val = filter.threshold_otsu(camera)

END EDIT

Could you threshold the amplitude of the recording? say this is our recording

where we are just plotting

plot(abs(sound))

can we say the chirping is probably above 2? and assume everything below magnitude=2 is just noise? If so, and assuming you have the imageprocessing toolbox, you could

normalize the audiofile amplitude (between 0-1) needed for next step
get automatic thresholding of file imageprocessing toolbox function
find indicies that exceed threshold and say these are bird chirps

the only problem is with sound even with a strong signal, it will usually cross 0 so you would have discontinuities even in the actual bird chirping. Maybe you can add some sort of hysteresis to the masking operation. Like
if the current sample is below the threshold BUT a large percentage of the previous few samples were above the threshold, we will assume this sample is also an interest point

roughly your code would be something like this (im not near a comp with matlab, so this may have slight syntax errors). This is also not very optimized, but maybe it's a springboard for you or others

%% step 1
%we are using magnitude only
my_min = min(abs(soundfile));
my_max = max(abs(soundfile));

%gets sound file magnitude between 0-1 only
normalized_sound_mag = 1/(my_max-my_min) * (abs(soundfile) - my_min);

%% step 2
%get a threshold, we needed the magnitude between 0-1 for this function to work
sound_thresh = graythresh(normalized_sound_mag);

%% step 3
%this vector will store all sample indexes that are interesting
idx_sound = 0;

num_for_hysteresis = 10;        %the number of samples to use for hysteresis
hysteresis_percentage = 0.8;    %percentage of samples that must be above threshold
                                %so say in our set of 10, 2 samples are below threshold
                                %since .8 are above it, we say they are all of interest

%because the way we implement hysteresis we have to skip the first few samples
for (i=num_for_hysteresis:1:length(normalized_sound_mag))
    %if we are above thresh, without a doubt save the index
    if (normalized_sound_mag(i) > sound_thresh)
        idx_sound = [idx_sound, i];
    else
        %hysteresis check prev samples, creates a boolean vector. 1 means value above thresh, 0 means it was below
        samples_above_thresh = normalized_sound_mag(i-num_for_hysteresis:i) > sound_thresh;

        %nnz, gets the number of nonzero elements in a matrix
        num_prev_samples_above_thresh = nnz(samples_above_thresh);

        %if the number of samples in the prev window met our criteria, this current sample
        %is probably of interest as well
        if (num_prev_samples_above_thresh > hysteresis_percentage * num_for_hysteresis)
            idx_sound = [idx_sound, i];
        end
    end
end

%idx_sound should now have all the indicies f interest, these can now be used
%directly on the original soundclip

we could also use the frequency power and threshold, rather than the sound amplitude, the basic outline would still be the same.


You are probably interested in the real discrete transforms. Here are the basis functions for a transform of length 6. The inverse transforms are weighted sums of these basis functions, and the forward transform figures out the weights given the time domain sequence. A reconstruction of the time domain sequence by the inverse transform is exact. The discrete transforms don't actually involve the continuous functions shown below but rather their samples at integer time indices (here 0, 1, 2, 3, 4, 5).

As a service to those who prefer to see equations, the discrete Hartley transform $H_k$ and discrete Fourier transform $X_k$ of $x_k$, where $k = 0, \dots, N-1$ are defined by:

$$H_k = \sum_{n=0}^{N-1} x_n \left[ \cos \left( \frac{2 \pi}{N} n k \right) + \sin \left( \frac{2 \pi}{N} n k \right) \right]\\
X_k = \sum_{n=0}^{N-1} x_n \left[ \cos \left( \frac{2 \pi}{N} n k \right) - i \sin \left( \frac{2 \pi}{N} n k \right) \right],$$

where $i$ is the imaginary unit.

Hartley

Figure 1. Hartley transform basis functions are harmonic sinusoids. Each basis function is a sum of sine and cosine.

Being harmonic means that the frequency of each is an integer multiple of the lowest non-zero frequency. Some of the sinusoids have frequency greater than the Nyquist frequency.

Or as discrete sequences: [1, 1, 1, 1, 1, 1] (yellow), [1, 1.366, 0.366, -1, -1.366, -0.366] (black), [1, 0.366, -1.366, 1, 0.366, -1.366] (green), and so on.

Figure 2. Hartley basis functions, constructed from the discrete basis functions by perfect bandlimited interpolation assuming that the sequences repeat forever. The discrete Hartley transform could alternatively be expressed directly with these functions, colored the same as in Fig 1.

Except for the special cases of frequency 0 and π, the aliased Hartley basis functions come in pairs of harmonic equal frequency sinusoids that are 90 degree apart, just like in the real discrete Fourier transform below. Except for scaling and phase shifts (with one scaling constant and one phase shift constant for each pair), the basis functions are identical to those of the real discrete Fourier transform.

Fourier

Figure 3. Real discrete Fourier transform basis functions are harmonic cosines and sines, here colored the same as the equivalent scaled and phase-shifted Hartley basis functions in Fig. 2.

For frequencies 0 and π there are no sine basis functions because their samples would be zero-valued and would never contribute to the weighted sum in the inverse tramsform. Unlike with the Fig. 1 formulation of Hartley transform, the discrete Fourier transform basis functions stay within the Nyquist frequency limit.

Popularity of Hartley

The fast Hartley transform algorithm came out in 1984, but there has been no significant increase in its popularity since.

Figure 4. Google Ngram of Hartley transform counts the relative frequency of the word pair in literature as it developed over the years.


Is it possible to design noncausal filters in Matlab?

The short answer is: yes, of course it's possible to design non-causal filters. Consider a simple first order IIR filter:

$$y[n]=ay[n-1]+bx[n]\tag{1}$$

If implemented in the way as indicated in Eq. (1), this filter is causal. It is also stable if $|a|<1$. Now rewrite (1) as

$$y[n-1]=\frac{1}{a}y[n]-\frac{b}{a}x[n]\tag{2}$$

Implemented in this way, you compute $y[n-1]$ from the future samples $y[n]$ and $x[n]$, so the system is anti-causal. Note that it is stable if $|a|>1$. Similar considerations apply to higher-order filters.

One way of realizing non-causal IIR filters is implemented in Matlab's filtfilt function. How this works is described in this answer.

Non-causal FIR filters can always be implemented by appropriate buffering of the input signal and delaying the output signal.


EDIT: I took this answer down for a time because I realized a lot depends on whether there are any idiosyncrasies in how images are mapped to pixels by actual devices, and, not being an image guy, I don't know a lot about that.  I decided to bring it back, with the caveat that the answer may be insufficient given said idiosyncrasies.

  The image should be upsampled before supersampling. But how much?

You are just trying to soften quantization (pixel resolution) error by collecting information from neighboring pixels.  Thus, oversampling by a factor of 2 should be plenty.  In other words, once you oversample until your sample resolution is twice that of the pixel resolution, nothing can be gained by oversampling more.

  Is there any special filter to be used for upsampling prior to supersampling?

I might be misunderstanding the terms, but it appears to me that upsampling is part of supersampling.  No, you don't need to use any special filters.  Insert 0's in between your samples and then low-pass filter.  You'll probably need to adjust the gain of your low-pass filter by your upsampling rate to get the pixel intensity to come out right.

  What if I want to resize image to e.g. 80% of its size with supersampling?

Then factor that into your upsample/decimation calculations.  After you upsample you have to decimate down to the pixel resolution.  By resizing the image you are, in effect, changing the pixel resolution.

  How much I have to upsample it and how to treat non-integral size ratio?

How much you have to upsample it depends on the details of the image, pixel resolution, etc.

  How much I have to upsample it and how to treat non-integral size ratio?

You can deal with non-integral ratios with fractional resampling.  Essentially you upsample by an integer and decimate by another integer, giving you a fractional sample change.


How would I find the function given the magnitude plot and the phase response?

There is more than one way to write the same thing. You can use $rect()$ functions or define it piece-wise, but in the end it is the same thing, since the $rect()$ is just a shorthand for a piece-wise definition.

If you do it in sections, you don't need to use $rect()$ function at all, you could simply write is as

$$X(j \omega) = \begin{cases}
 &  j \text{ if } -\pi < \omega < 0\\
 & -j \text{ if } 0 < \omega < \pi \\
 & 0  \text{ otherwise}
\end{cases}
$$

The answer just puts it together as two rectangle functions: one for negative frequencies and one for the positive ones.

Either way would be correct and I don't think one is "better" than the other unless there are some clear criteria what "better" means.


How to perform spectral inversion in the frequency domain to convert a low-pass filter into a high-pass filter?

The flaw in your logic is assuming that the frequency response of your lowpass filter in its passband is identically $1$. You're inherently assuming this when thinking that you can form the highpass filter by simply calculating $H_{hp}(f) = 1 - H_{lp}(f)$.

However, the lowpass filter has group delay because the peak of the sinc function doesn't occur at $n = 0$. Instead, as a symmetric FIR filter, it will have a delay of $\frac{N}{2}$ samples, where $N$ is the filter order. This means that the frequency response of the filter across the passband may have unity magnitude, but it will have linear phase! Recall the time shift property of the discrete-time Fourier transform:

$$
x[n-k] \Leftrightarrow X(f)e^{-j2\pi fk}
$$

So, instead of your prototype lowpass filter having a response of $1$ in the passband, it actually has a response of $e^{-j2\pi fk}$, which imparts a phase shift that is linearly proportional to frequency. If you would like to transform this filter to a highpass structure with equivalent group delay, you would need to take this into account when synthesizing your highpass filter in the frequency domain:

$$
H_{hp}(f)= e^{-j2\pi f\frac{N}{2}} - H_{lp}(f)
$$


How is wavelet center frequency computed?

Short version: DFT's bin indices are input length-dependent; "center frequency" is measured relative to the function generating the wavelet. Generated length can vary, so must be accounted for to yield 'correct' $f_c$. This requries injecting information from knowledge of the function, as DFT is blind to "absolute duration"; PyWavelets does this via domain, the range of input values to the function. domain=16, so DFT's bin indices are x16 shifted; undo by dividing.

The answer's in DFT. Take Real Morlet, with the $\sigma >=5$ approximation:
$$
 \psi(t) = e^{-t^2/2}\cos(\sigma t) \tag{1}
$$
Suppose we evaluate it for $\sigma=5$, $t\in[-8,8]$; what's it look like?

So what is its "center frequency", or frequency of a sinusoid with which the wavelet correlates most strongly? Try tweaking $\beta$ here. Spoiler: sigma; the exponential is almost ignored, just slightly shifting optima locations (but strongest correlate still is $\sigma$):

As there's no $2 \pi$ in $(1)$, 5 is the angular frequency, so $f=5 / 2\pi=0.796$, or $f^{-1}=1.26$; these match the red cosine, which cycles once every 1.26 units, or 0.796 times every 1 unit.
But in the discrete realm, we've no notion of "time" unless explicitly labeled; all we have is samples. What happens if we take the DFT of red? Do we peak at $5$, or $5 / 2 \pi$? Neither; both of these numbers are "per second", and DFT doesn't understand that. All it knows is that a cosine completing one full cycle from left to right ... completes one full cycle from left to right. It calls this $k=1$.
Then we seek $k$ such that the cosine matches most closely with red. Since red makes 0.796 cycles every 1 unit, and there are 16 units, then it'll make 0.796 * 16 = (1/f) * domain = k cycles in total. That's the index of max(abs(DFT)); we just swap f and k in the equation, inputting k to estimate f.
All said of red above also applies to blue; their peaks match. Putting it all together:
f_center = (t[-1] - t[0]) / k_max = domain / k_max

But this looks different from PyWavelets
Because PyWavelets looked ahead. First, we never expect our wavelet to be flat, so we eliminate the possibility of k_max=0 in case the wavelet averages nonzero:
argmax(abs(fft(psi)[1:]))  # _k_max - 1

But this shifts all k, such that original k=1 is now k=0; so we do +1:
index = argmax(abs(fft(psi_pad))[1:])) + 1  # _k_max

What if $k > N/2$, where N == len(psi)? Negative frequency. This is fine, except we seek the counts of cycles, which are non-negative, so we just find its positive. Suppose N=101; then, there is no Nyquist frequency, only max positive and negative frequencies (equal) at k=50, 51.

_k_max = 100. This is negative k=1. To obtain positive: N - _k_max = 101 - 100 = 1.
_k_max = 50. This is the positive max; do nothing.
_k_max = 51. This is the negative max; make it (+): N - _k_max = 101 - 51 = 50.

This explains
if index > len(psi) / 2:
    index = len(psi) - index  # now we get actual `k_max`

And that's all. Except it's not; PyWavelets does domain / (index - 1), and +2 instead of +1, and -index + 2. The behavior is identical for most, but not all, possible len(psi) and _k_max; grid search. Consider (len_psi, _k_max) = (101, 50). pywt finds k_max=51, which is a negative frequency.

How accurate is it?
Not even within two significant figures. Authors may have mistaken higher precision for higher precision, when in reality k_max gets scaled down by the exact amount psi gets scaled up. (Exception for precision < 5, for which estimation greatly degrades):

Can we do better? -- Yes; above right is the result. See here. Opened a PR, also on k_max.

More rigorously, but briefly: center frequency isn't necessarily at max frequency component.
From Mallat 4.1, if $\psi$ is a wavelet, then $|\psi|^2$ is interpreted as a probability distribution per
$$
\int_{-\infty}^{\infty} |\psi (t)|^2 dt = 1
$$
centered in time at
$$
u = \int_{-\infty}^{\infty} t |\psi (t)|^2 dt
$$
(i.e. expected value, $\mathbb{E}[|\psi(t)|^2]$) and in frequency ("center frequency") at
$$
\omega_c = \frac{1}{2\pi} \int_{-\infty}^{\infty} \omega |\hat{\psi} (\omega)|^2 d\omega
$$
Corresponding variance relations are given, from which we define Heisenberg boxes for time-frequency resolution:

If a wavelet's FT happens to be symmetric about the peak (as with most used wavelets), then center is at peak's location.

Addendum: 'relative to what' is $f_c$? How to interpret?
Relative to "the function", surely, but how to interpret this in terms of signal data we work with? Like with DFT, to find "physical" values, we need to know the sampling frequency, $f_s$. Physical $f_c$ has units of Hertz, so if we know the number of samples that comprise one second, we can relate "frequency relative to samples" (DFT) to "frequency relative to seconds" (more here). Thus:
$$
f_c \left[ \frac{\text{cycles}}{\text{second}} \right] = \left( f_c \left[ \frac{\text{cycles}}{\text{samples}} \right] \right) \cdot \left( f_s \left[ \frac{\text{samples}}{\text{second}} \right] \right)
$$