By doing PCA, the principal components you will get will not correspond to a single recording, but rather to a mix between them. PCA is a feature extraction method, whereas what you are looking for, seems to me, as a feature selection problem.

Also, if you have so many simultaneous recordings, all affected by the same sources of noise, why not perform active noise cancelation?


The use of LS estimation with Vector OFDM

Sorry for late reply .. I was little bit busy.

you have mistakes in your code. Although mythology is right, you have mistakes in some parameters.

Check this paper "Low-Complexity Equalization of Orthogonal Signal-Division Multiplexing in Double-Selective Channels" Then,  Modify your code following part II in that paper. It's with details there, if you got any issue, you'd let me know.
Regarding the channel estimation, you follow the same way of conventional OFDM scheme except for V-OFDM you use Vectors pilots instead of symbol pilots.  I means in that mentioned paper, you will have vector of size $M$ used as pilots.

Good luck


Hi, If I really want to focus on the mathematics of signal procesing, Should I do a Masters in EE or Math?

As stated in the comments, this is more about the content than the label. The best is to play with both fields if you want to invest in the future: signal processing is borrowing  increasing from mathematical tools, for instance convex and non-convex optimization, learning, topological and geometrical tools (for dimension reduction), graphs, etc., that might  not be taught in traditional EE. Such techniques can help with the growing quantity of data with little models and knowledge.

As a personal anecdote, I attended a lecture on "math for EE engineers" for the beginning of a signal processing PhD program in a US university. It was about Euclidean vector spaces (in finite dimension of course). A young PhD student asked whether the basis vectors did belong to the space. I was quite surprised by the level of the question. But not as much as when I hear the teacher say, after a long pause: "mmmm, I'd say... yes, let me think about it" (such a question is indeed interesting in infinite dimension spaces). There are similar  stories from colleagues attending other MsC programs.

So after picking the right topics, I strongly suggest  that you either complement an EE program with maths, to get better insight and intuition about the theoretical    tools, or dive in maths, while practicing a lot on "real DSP problems" aside. The gap is large between mastering  Fourier theory and being able to properly estimate and interpret  of a power spectrum of a noisy signal.
There are plenty of online resources to assist you in that task.

I really appreciate the concept  behind the law of excessive learning in mathematics, stated by Alexander Borovik:

  To be able to use maths at certain level it is necessary to learn it
  at the next level


Analyze and reproduce sonic screwdriver sound

An approximation

Actually, the signal is not that complicated in my opinion. However, I am not a sounds guy... Just playing around with it a bit.

Interestingly, I have a very similar signal that I have produced (see attached sound file and spectrograms. Spectrogram above is my signal (somewhat tuned to your signal's parameters) and below is the sonic screw-driver (of which Doctor by the way? ;) )

I cannot post the code here, since as I said the parameters are a bit tuned to your signal. What I basically did is to create a signal that is...

A sum of 10 signals. Center frequency f_mid = 3.7 kHz + 20Hz*N, N = 0, 1, ..., 9, just to have some kind of a wobbly jitter in there
The center frequency is modulated by a "circle" (df: maximum change of the center frequency, t_circ: duration of a circle, t_mid: starting point of the circle in time) with df = 20 Hz, t_circ = 0.1s. The phase of the signal is generated from the constant frequency plus the dynamic modulation and a very simple numerical integration (phase multiplied by an increasing dt-vector and cumsuming over it). To get to a signal's frequency, you would differentiate its phase wrt to time - so this is the inverse process.

T = 1.0/44100.0
f_circ = df*(1.0 - ((t_circ - t_mid)/t_mid)**2.0)**0.5
phase1 = np.cumsum((f_mid + f_circ)Tnp.ones(t_circ.shape))
phase2 = np.cumsum((f_mid - f_circ)Tnp.ones(t_circ.shape))
circ = np.sin(2.0*pi*phase1) + np.sin(2.0*pi*phase2)

Added some somewhat "Brownian" Noise, however it is much below the noise in your signal.

Things to explore

The "duration" of the "circles" seems so vary slightly. Probably sinusoidal variation is also fine, so the period could be jittered.
Apparently, there is additionally a very broad noise-band. Maybe try some uniform noise with a time-varying windowing?
Subtract a bunch of short-term, constant frequency signals from your signal. There is a nearly some kind of a grid visible in the spectrogram.

Some code samples

I am sorry, in my previous description, I apparently have missed some details. Please find here a code snipped that should be able to produce a basic signal that exhibits the basic "circular" frequency change:

import matplotlib
matplotlib.use('MacOSX')
import matplotlib.pyplot as pp
import numpy as np
import scipy.io.wavfile
import scipy.signal

# Export function for convenience
def export_wave2(t, signal_l, signal_r, name = "test_sp.wav"):
    sample_rate = 44100
    data = np.stack((signal_l, signal_r))
    data = np.transpose(data)
    scipy.io.wavfile.write(name, sample_rate, data.astype("float32"))

# Sampling frequency
fs = 44100
T = 1.0/fs

# Circle parameters
t_circ_max = 0.2 # duration of one "circle"
t_circ     = np.arange(0, t_circ_max, 1.0/fs)
t_mid      = t_circ.max()/2.0
f_mid      = 3700 # middle frequency of the sound signal
df         = 400     # amplitude of the frequency modulation, ie frequency will vary between f_mid - fd and f_mid + fd

# Create ONE circle
f_circ = df*(1.0 - ((t_circ - t_mid)/t_mid)**2.0)**0.5

phase1 = np.cumsum((f_mid + f_circ)*T*np.ones(t_circ.size))
phase2 = np.cumsum((f_mid - f_circ)*T*np.ones(t_circ.size))
circ = np.sin(2.0*np.pi*phase1) + np.sin(2.0*np.pi*phase2)

# Repeate the circle N times, ie the signal will be N*t_circ_max long
N = 10
circs = np.hstack([circ for k in range(0,N)])
circs = circs / np.abs(circs).max()
circs_time = np.arange(0, N*t_circ_max, 1/fs)

pp.plot(circs_time, circs)
pp.show()

export_wave2(circs_time, circs, circs, name = "test_sp.wav")

The frequency behaviour is of course not visible directly in the plots, but when a spectrogram is calculated from circs. I have loaded the test_sp.wav into Audacity for this.

Sound example


How can I generate audio wave files of harmonic series/partials for a specific frequency in Matlab

use a for loop to save each frequency into a different named wav file, use string manipulation command num2str() of Matlab. The following will help:

clc; clear all; close all;

Fs = 44100;         % Sampling frequency
F0 = 200;            % base frequency
t = [0:1/Fs:0.1];   % time index
k = 0;              % harmonic index

for fsine=F0:F0:Fs/2
    x = sin(2*pi*fsine*t);
    filename = ['wav' num2str(k)];
    wavwrite(x, Fs,16, filename);
    k = k+1;
 end


How to locally quantify the 'sharpness' of an image?

The recent works I am aware of make use of tools that go beyond mere gradients. Here are a few references that could be starting points:

S3: A Spectral and Spatial Measure of Local Perceived Sharpness in Natural Images, 2012, with examples of sharpness maps and Matlab code (that could  be converted to Python)

  This paper presents an algorithm designed to measure the local perceived sharpness in an image. Our method utilizes both spectral and spatial properties
  of the image: For each block, we measure the slope of the magnitude
  spectrum and the total spatial variation. These measures are then
  adjusted to account for visual perception, and then, the adjusted
  measures are combined via a weighted geometric mean. The resulting
  measure, i.e., S3 (spectral and spatial sharpness), yields a perceived
  sharpness map in which greater values denote perceptually sharper
  regions

Image sharpness assessment based on local phase coherence, 2013, with examples  and (BROKEN yet) code (local copy of the code provided by the authors)

  Sharpness is an important determinant in visual assessment of image
  quality. The human visual system is able to effortlessly detect blur
  and evaluate sharpness of visual images, but the underlying mechanism
  is not fully understood. Existing blur/sharpness evaluation algorithms
  are mostly based on edge width, local gradient, or energy reduction of
  global/local high frequency content. Here we understand the subject
  from a different perspective, where sharpness is identified as strong
  local phase coherence (LPC) near distinctive image features evaluated
  in the complex wavelet transform domain. Previous LPC computation is
  restricted to be applied to complex coefficients spread in three
  consecutive dyadic scales in the scale-space. Here we propose a
  flexible framework that allows for LPC computation in arbitrary
  fractional scales.

The given examples and comparisons across different  could provide you with  some hints toward your goal.


How to find out the farrow coefficients if FIR coefficients are given?

There is no standard way to determine the coefficients of the Farrow structure for a given FIR filter. The Farrow structure is an implementation of a whole class of FIR filters with transfer function(s) $H_{\mu}(z)$ with a continuous control parameter $\mu$. Often this parameter $\mu$ determines a fractional delay, but it could as well be any other adjustable filter property.

The transfer function $H_{\mu}(z)$ of the Farrow structure is given by

$$H_{\mu}(z)=\sum_{k=0}^{K-1}C_k(z)\mu^k\tag{1}$$

where $C_k(z)$ are FIR transfer functions. From $(1)$, the impulse response is

$$h_{\mu}[n]=\sum_{k=0}^{K-1}c_k[n]\mu^k\tag{2}$$

So each filter coefficient (of the corresponding transversal filter structure) $h_{\mu}[n]$ is implemented as a weighted sum of coefficients $c_k[n]$. This is why there is no one-to-one mapping from given FIR filter coefficients to the coefficients of the Farrow structure.

You need to specify which filter property the parameter $\mu$ is supposed to control, and then you need to design the coefficients $c_k[n]$ of the Farrow structure accordingly.


How can I test an implementation of DFT for correctness?

A very quick and easy test that catches all relevant structural mistakes follows from a fundamental property of the Fourier transform (and also its inverse): It is a square root of the sequence order reversal operator. That means applying your implementation twice on a data vector will reverse the sequence order and applying it four times will give your original data back.

Note that this requires the so called unitary scaling convention, where the Fourier transform and its inverse share the same scaling factor. Otherwise you will get back a scaled version of the original signal.


When you multiply the image matrix with the rotation matrix, you are essentially multiplying an integer matrix with a real number matrix. The result is a real matrix, which has to be quantized back to a discrete matrix.

What is being quantized? The coordinates themselves. The integer coordinates of the previous position of the pixel is now a real valued coordinate, which must be converted back to an integer in order to be placed on the screen.

If you were asked to rotate it back to the original angle you could have been right in questioning the futility of interpolation. However, the user might not want to rotate it back, so you have to perform interpolation or quantization to be able to use the rotated image. If you are concerned about image fidelity, you might want to save a copy of the original real valued matrix in memory or file until or opertions are finished and the user exits program.

Here is a useful link you might want to peruse.


The thing you're looking for is known as an "all-pass filter", which has a flat amplitude response, but modifies the phase relationships within a signal. This would be a separate stage from the one that limits the bandwidth.

Note that the 50-500Hz bandpass filter is going to modify phase relationships, too, so the all-pass filter will need to be designed so that the overall system response is what the client is looking for.

However, it's possible that the client just wants a time delay, but expressed it in terms of the phase angle of some fundamental frequency. This will need to be clarified.


Basic difference between Fourier transform and laplace transform?

Fourier transform is an intuitive tool that's a bridge between domain of physics and mathematics, as it quantitatively describes the periodic content of the signals and also frequency response characterisation of systems that occur in physical (and engineering) applications. The use of frequencies is quite intuitive and consistent at least for stable systems...

However for unstable systems (and signals) Fourier transforms becomes mostly awkward (if not useless) to deal with. However control engineers unavoidably must make frequent use of unstable systems in their works. And for this purpuse, Fourier transform is either insufficient or awkward, hence a generalisation of the existing Fourier transform is made into the Laplace transform which conveniently yields mathematical (complex algebric) descriptions of stable as well as unstable systems which was not possible with the Fourier. The Laplace transform, therefore, includes a region of convergence parameter into it.

Another difference between the two transforms is in the time-domain transient analysis of output of LTI systems driven under nonzero initial conditions which is successfully captured in the Laplace transform only. In the sense that LCCDE with initial conditions are straightforwardly solvable by (unilateral) Laplace transforms whereas the standard FT can only solve LCCDE with zero initial conditions (initial rest)...

For one sided and two sided differences, I think Stanley has things to say...


If MFCC vectors are speaker independent, then how can we reconstruct unique speech signals from them?

First of all there is some serious "cheating" in the MFCC reconstruction experiment you linked to: not only the MFCCs are used, but also the voiced/unvoiced bit and the pitch.

MFCC are not speaker independent. In fact, they are used for speaker identification/verification tasks!

The speaker "idiosyncracies" are both in their prosody (preserved by this reconstruction experiment because the pitch is provided as a side-information to the reconstruction process) and in the articulation/timbre (preserved by the MFCC).

Two ingredients are needed to get MFCCs to work for speaker-independent recognition:

Vocal tract length normalization. A linear transformation (matrix multiplication of the MFCC vector) can map relatively well the MFCC sequences of two speakers speaking the same sentence. So even if MFCC are not speaker-independent one can optimize for a transformation matrix that "flattens out" the speaker-specific details.
Acoustic modelling. Using a large number of gaussians (or any classifier with large capacity) for a specific acoustic unit allows it to capture all the variations.


ICA - Statistical Independence & Eigenvalues of Covariance Matrix

Signals 3 and 5 appear to be quite correlated -- they share their first harmonic. If I were given two mixtures of those, I wouldn't be able to separate those, I'd be tempted to put the common harmonic as one signal and the higher harmonics as a second signal. And I would be wrong! This might explain the missing eigenvalue.

Signals 1 and 2 do not look independent too.

A quick and dirty "sanity-check" for the independence of two series is to do a (x, y) plot of one signal against the other:

plot (sig3, sig5)

and then to do the same (x, y) plot with one signal shuffled:

indices = randperm(length(sig3))
plot(sig3(indices), sig5)

If the two plots have a dissimilar look, your signals are not independent. More generally, if the (x, y) plot of the data shows "features", disymmetries etc, it's a bad omen.

Proper tests for independence (and those are the objective functions used in the ICA optimization loop) include, for example, mutual information.

ICA is recovering the most independent signals a linear mixing of which yields your input data. It will work as a signal separation methods and recover the original signals only if those were maximally independent according to the optimization criterion used in your ICA implementation.


How can I extract a sinusoid from a very noisy signal for computing the phase shift?

You are correct. Estimating the instantaneous phase of a noisy sinusoid is NOT easy.  I suggest you design a narrow bandpass filter such that your sinusoid-of-interest is in the filter's passband.  (The better the filter the more noise that will be eliminated.) Pass your two signals through the bandpass filter to generate filtered signals $x_1[n]$ and $x_2[n]$.  Next, pass your $x_1[n]$ and $x_2[n]$ signals through a Hilbert transformer to generate $\hat{x}_1[n]$ and $\hat{x}_2[n]$.  Create two analytic (complex) signals as:

$$z_1[n] = x_1[n] + j \, \hat{x}_1[n],$$

and  $$z_2[n] = x_2[n] + j \, \hat{x}_2[n]$$

where

$$ \begin{align}
\hat{x}[n] & = \mathcal{H}\{ x[n] \} \\
           & = \sum\limits_{m=-\infty}^{+\infty} \frac{1 - (-1)^{m}}{\pi \, m}   x[n-m]  \\
\end{align} $$

Next, compute two instantaneous phase sequences:

$$\phi_1[n] = \arg\{z_1[n]\}$$

and $$\phi_2[n] = \arg\{z_2[n]\}.$$

Finally, compare the instantaneous phase difference between the $\phi_1[n]$ and $\phi_2[n]$ sequences.


How to choose the size of a Laplacian of Gaussian kernel for filtering images fast?

you can approximate LoG with Difference of Gaussians. Apply 2 Gaussians (separable) with different sigmas. (for what I know with sigma ratio 1.1) and take the difference.


What are bits in the context of channel capacity units?

in Information Theory, the amount of information in a message is proportional to the logarithm of the probability of occurrence of the message:

$$ I(m) = -C \log\left( P(m) \right) $$

the only difference between different bases of logarithms is a scaling constant.  if $ C = \frac{1}{\log(2)} $, that is

$$ \begin{align}
I(m) & = -\frac{1}{\log(2)} \log\left( P(m) \right) \\
     &= - \log_2\left( P(m) \right) \\
 \end{align} $$

then we say that the $I(m)$ is in units of bits.  the reason why is if the message $m$ had only two possibilities and both had equal likelihood, then $P(m) = \frac12 $ and $I(m) = 1$.  there are two messages, two possibilities of equal likelihood, like heads and tails of a coin.  and you need only one bit, 0 or 1, to represent the two possibilities.


Am I supposed to normalize FFT in Python?

Your decision to normalize or not does not change the accuracy of your answer, as it is simply a scaling factor. If you use the common scaling of $1/N$, then the output for each DFT bin will represent the average of the portion of the input signal that is at the frequency defined by that bin, scaled to the same units as the input. So that is convenient and gives a certain meaning to what the magnitude of the output represents, but if you didn't divide by N it doesn't make the answer incorrect, as long as your approach is consistent.

For example, it is also common to use $1/\sqrt{N}$ as this would make the DFT and IDFT completely symmetric as the same scaling can be used for both.

To understand what is happening, consider the simplest case of the computation of the first bin of the DFT, which represents DC, and consider a DC signal that is 1 for all $N$ samples given. If this DC signal represents a magnitude in voltage, then this is samples of $f(t) = 1$ volt, given as $f[n]=1$, for all $N$ samples, $n = 0$ to $N-1$.

The first bin of the DFT without any normalization is simply the sum over all N samples:

$$F[k=0] = \sum_{n=0}^{N-1}v[n] = N$$

So if we were to scale this by $1/N$, then the result would represent the average for the frequency component of our signal within that bin in the same units as the input. In this case all values were at bin 0 "DC", so the average value was 1V, consistent with what we get when we scale by $1/N$. (It is clear from the formula above that when scaled by $1/N$ this is simply the average of $f[n]$).

Similarly all the other bins are the result of the average after multiplying by the time domain signal by $e^{-jk\omega_o n}$, which can be viewed as frequency translating any signal at that bin to DC and then taking the average of that result.


Negative group delay and envelope advance

Answer : No, any causal LTI system with frequency response $H(f)$ cannot produce the output $y(t)$ in advance. And, the answer lies in the causality of input signal $x(t)$ being applied to $h(t)$. Any causal input $x(t)$ which has an identifiable beginning cannot truly be Narrow-Band or Band-Limited. It will have non-zero frequency content at all frequencies.

Yes, you are right that it is not uncommon to have practical LTI systems with $+ve$ gradient of $\angle H(f)$ in parts of the response and hence making group delay $-ve$ around those parts of $H(f)$. And, if we can give a Narrow-Band input $x(t)$ such that the bandwidth of $x(t)$ is restricted in that part of $H(f)$, then you would have a time-advanced output. So, are we able to look into the future?

NO!!!! We are not. My point will get clear in a minute.

Let me take an example of a very common and practically realizable IIR filter in equivalent discrete time scenario : the Leaky Integrator.

The $H(e^{j\omega})$ of a leaky integrator is given by the following:
$$H(e^{j\omega}) = \frac{1-\lambda}{1-\lambda e^{-j\omega}},$$So,
$$|H(e^{j\omega})|^2 = \frac{(1-\lambda)^2}{1 + \lambda^2 -2\lambda cos(\omega)}, \angle{H(e^{j\omega})} = arctan \{\frac{-\lambda sin(\omega)}{1-\lambda cos(\omega)} \}$$

The shape can be plotted in MATLAB by following:

  freqz(0.1, [1 -0.9], (-pi:0.001:pi));

Now, if we can give a very narrow-band input $x[n]$ centered around $\omega = 0.6\pi$ and bandlimited within a very small $\Delta \omega$, then we would get a response as follows:

$x[n] = s[n]cos[\omega_o n]$, where s[n] is a narrowband baseband signal and $\omega_o = 0.6\pi$ and group delay of the filter is $g_d$ around $\omega_o$
$$Y(e^{j\omega}) = X(e^{j\omega}).e^{-j.g_d(\omega-\omega_o)},$$ You can work this out to get $y[n] = s[n - g_d]cos[\omega_o n]$

According to the equation above, leaky integrator is basically producing an output which is having a delayed envelope of input by $g_d$ samples. And, what happens if this $g_d$ is negative!

Check out that $g_d$ is indeed negative around $\omega_o = 0.6\pi$. Does that mean that the leaky integrator is able to produce the $s[n]$ envelope $g_d$ samples in advance?

No, it is not. The caveat is that we cannot have a perfectly bandlimited narrowband causal input $x[n]$. We cannot have a $x[n]$ which has an absolute start in time and yet it is having a bandlimited narrowband response in frequency domain.

Because we cannot have such input $x[n]$, hence we cannot have a "future seeing time machine".

In order to produce a causal input, which has an identifiable absolute start in time, the frequency response of the input will spread in frequency domain and the input $X(e^{j\omega})$ will be present at all frequencies with non-zero spectral components, and this will make the overall delay to be positive.

Indeed, if you plot the group delay response of the leaky integrator, you get the following, and check that even though the group delay is small negative number away from $\omega = 0$, it is taking high $+ve$ values around $\omega = 0$:

Hope that answers your question.


The ambiguity function is used in radar systems to get the distance and the relative speed of a moving object with respect to the transmitter. Is called ambiguity because it also tells about the ability to distinguish objects that are close between them and with a similar vector velocity.

The ideal ambiguity function is a Delta function in both domains. However there is no waveform with this ambiguity function. With your Barker-Code, you could have troubles to identify different objects that are at a distance +/- <10us, and whose Doppler is +/- <10 kHz. But you could distinguish two objects that are at 5us, one with 5 KHz and the other with 27 kHz.

In order to estimate the Doppler and the delay, you first compute the ambiguity function between the Barker-Code and the reflected signal. You will get a noisy version of the above figure but shifted in both domains, depending on the delay and Doppler of the target. If there are more targets, you will obtain overlapped ambiguity functions, each one shifted in a different place on the correlation plane.

You then just find the position of the maximum. From the obtained Doppler you can get the speed with the equation you wrote.


What is the rate of expansion of an image as it aproaches the human eye?

The magnification of an object at distance $g$ from a thin lens with focal length $f$ is given as

$$\beta=\frac{f}{g-f}$$

which is generally negative for objects at distances greater than the focal length. This means the image is upside down.

The derivative of the magnification with respect to the object distance is

$$\frac{\partial \beta}{\partial g}=-\frac{f}{(g-f)^2}$$

So for constant $f$ the magnification diverges like $x^{-2}$ where $x=(g-f)$ is the difference between focal length and object distance.

So yes, the magnification grows in an accelerated way and even approaches infinity.

However, the relationship is only strictly true if the image at the image plane is sharp and the aperture is unlimited. So for your eye, at some point the image will become blurry and the magnification will grow in an asymptotic way that is much harder to describe.


Given a noise-free channel with limited bandwidth, why is the channel capacity infinite?

Why is this intuitively true?

It's a purely hypothetical mathematical formula result. So it doesn't even need to be intuitive, at all.

It's very normal for humans to look for intuition in abstract things, but an observable noise-free channel doesn't exist, so meh.

But if you insist: in a noise free channel, a single symbol could take one of $N$ values, with $N$ becoming arbitrarily large, because, due to the absence of noise, you could tell arbitrarily fine voltage differences apart. With that, a single (equally probable) symbol $x$ would have Shannon info $I = -\log_2 p_x = -\log_2 \frac1N = -(\log_2 1 - \log_2 N) = \log N\rightarrow \infty$ information.

EDIT: I saw you just got the same answer not an hour ago:

  But in terms of pure math, there is no ($\aleph_0$) finite limit on the amount of information that a single real number (one sample) can contain.


I know this answer is very late but I think this question worth an answer

The algorithm you used goes as follows:

Blur the original image using the Gaussian filter with given Mask Size & Sigma
Subtract the blurred image from the original (result is called Mask) to eliminate background and get the edges regions
Add a weighted portion from the mask to the original image by multiplying the Mask (the edges only) by K to enhance edges regions

K: User input (If K = 1 Unsharp, If K > 1 Highboost)

Mask(x,y) = Orig(x,y) – Blurred(x,y)

Result(x,y) = Orig(x,y) + K × Mask(x,y)

Here is the code:

public static Image<Gray, byte> Sharpen(Image<Gray, byte> image, int w, int h, double sigma1, double sigma2, int k)
    {
        w = (w % 2 == 0) ? w - 1 : w;
        h = (h % 2 == 0) ? h - 1 : h;
        //apply gaussian smoothing using w, h and sigma
        var gaussianSmooth = image.SmoothGaussian(w, h, sigma1, sigma2);
        //obtain the mask by subtracting the gaussian smoothed image from the original one
        var mask = image - gaussianSmooth;
        //add a weighted value k to the obtained mask
        mask *= k;
        //sum with the original image
        image += mask;
        return image;
    }

I think your problem is that you didn't subtract the smoothed image from the original image, anyway I hope this help


The signal $x[n]=3^n$ for $n\le -3$ (and zero otherwise, I assume) can be written as

$$x[n]=3^nu[-n-3]\tag{1}$$

To compute its $\mathcal{Z}$-transform you simply need to use the formula:

$$\begin{align}X(z)&=\sum_{n=-\infty}^{\infty}x[n]z^{-n}\\&=\sum_{n=-\infty}^{-3}3^nz^{-n}\\&=\sum_{n=-\infty}^03^{(n-3)}z^{-(n-3)}\\&=3^{-3}z^{3}\sum_{n=0}^{\infty}\left(\frac{z}{3}\right)^n\\&=3^{-3}z^{3}\frac{1}{1-\frac13 z},\quad |z|<3\tag{2}\end{align}$$

It's important to specify the region of convergence ($|z|<3$), otherwise the $\mathcal{Z}$-transform does not uniquely identify the corresponding sequence $x[n]$.


How to detect the edges of periodic signal?

The autocorrelation function of a sliding window of data may change by a significant amount when the window is slid past the first or last cycle of your burst of periodic events.  You may have to play with the width of the autocorrelation window depending on the estimated frequency of your cyclic/periodic events.


These terms refer to the general duality properties of these two domains. It should be obvious that higher frequency components mean faster changes in time and lower rise/fall times. In most natural systems, some of these fundamental properties oppose each other:
Shorter rise time means more overshoot and higher settling times and vice versa; lower overshoot means a higher settling time as well as a higher damping factor resulting in a higher rise time. In general, achieving better frequency domain properties (in the sense of time domain duals) means worsening other time domain properties. Apart from that, in specific applications such as digital transmission, cleaner pulses (lower rise/fall time) means broader spectrum and frequency content, which is undesirable in the sense of bandwidth efficiency. But again there is no unique answer to what you have asked. Each application must consider its own tradeoffs.


what does musical noise sound like and is there a sample clip for it?

Here's one example: http://www.izotope.com/tech/aes_suppr/

Googling "Musical Noise" certainly turns up a lot of noise music, so I'd suggest looking for academic papers that have web links to sound examples.


Why is scale space (DoG) needed to detect scale invariant features?

As you know, you find interest point (SIFT points) by finding the local maxima in scale-space, it mean the response of the detector must be maximum regarding the coordinates and also the scale. So knowing the scale you select features from that scale. Now consider the same process for same interest point in the same picture but with different zoom, the detector maxima will be in the same scale, so it will collect the same feature. So in the matching these two same feature will be matched together.

So the scale invariant is a result of finding maxima with respect to scale.


What is the standard way to use z-score with EEG?

Need a bit of clarity in question. There can be numerous normalization techniques for EEG depending on how they are used. That can be before processing (i.e. raw) or even after processing (for example, with FFT or other transforms).

Z-transform is ONE method. Period. There is nothing else called z-transform. You can have a look on Wikipedia entry of z-transform here.

Z-score, as used in your article, is a statistical measure, also called the standard score. More info here, again.

To sum up, z-transform and z-score are NOT linked.


How to understand FFT results of scipy.fftpack?

One rarely takes an FFT of a whole audio file, or a half. Usually only a sequence of shorter segments some number of milliseconds in duration.  The length of the "list of frequencies" will be half of the FFT length (for strictly real data) and the frequencies related to each y[i] will then be proportional to the sample rate and inversely proportional to your chosen FFT length.  Thus y[1] will only represent the output of a leaky narrow-band filter centered at 1 Hz if the FFT length is equal to the sample rate.

Decibels are more related to the log of the abs() or log magnitude of the FFT results. But without a calibrated microphone and DAC, the resulting decibels will be relative to some arbitrary level or offset. Normalization merely changes the arbitrary offset.


Has this 2D filter for enhancing circular dots in images a name?

This is a matched filter, where the shape of the filter coincides with the shape of the signal to be detected.

In this case, whenever the filter convolution overlaps a circular dot in the image, a maximum will be present in the output image. Thus, searching for the single maximum pixels in the output image will give the center position of dots in the original image.

Note that your linked source mentions that it is a matched filter.


Typically, for analyses like this, you do overlapping windows with a hop size that is fixed fraction of the window size.  So, for example, you might start by doing analyses of 1024 frames per window, and move the window forward 256 samples per analysis.  Depending on how the data looks then, you may or may not want to apply a window function (like Hamming or the like) to each window.

But since you're basically trying to remove noise, it might be simpler to take an example of the noise itself, and then try to reduce that noise from the original sample.  Take the FFT of the noise, subtract the frequency bins of the noise from the frequency bins of the sample, and run an inverse FFT to produce the final filtered result.  If you take this approach, make sure to only do FFTs at half the sample rate.


There are a couple of interesting aspects of "reconstruction to unity". First, there are two ways of combining two filters: parallel and in series. For a parallel topology it is ALWAYS possible to find a complimentary filter so that the pairs add to unity. It's easy enough, actually. Simply do $\tilde{H}(\omega) = 1-H(\omega)$. In the time domain that means that the impulse response of the complimentary filter is simply the negative of the original impulse response with 1 added to the first sample. So all the "ringy" stuff cancels out. Now the shape of this complimentary filter is not always what one would expect. For a 1st order low pass it actually is a first order high pass but for higher order filters it tends to have over/under swings in the cutoff region. However, it always exists as a stable causal filter.

Series (or cascade) "reconstruction to unity" is bit more complicated. Obviously the filters would have to be the inverse of each others, i.e. $\tilde{H}(\omega) = \frac{1}{H(\omega)}$. In general this can be done for any minimum phase filter. The inverse of a minimum phase filter is minimum phase as well and both are causal and stable.

So this leaves us with the question of how to interpret group delay in these cases. The cascade case is actually the more interesting. Since the filters are inverse of each other, the phase, and hence the group delay, of one is the negative of the other. So at frequencies were one filter has positive group delay, the other has negative group delay. An easy example would be a low shelf with +6dB of gain and a low shelf with 6dB of cut.
So negative group delays are very real and certainly not a violation of causality. In practice, these show up in areas of the filter that are fairly "non-flat" so the traditional interpretation of "delay of the envelope" doesn't quite apply since there is a fair amount of amplitude distortion as well.

If you Google "negative group delay", you can find a few IEEE articles that have tackled the subject.


inverse discrete FFT in python, multiple times?

The FFT and IFFT are almost identical operation: other than scaling (depending on convention) and  complex conjugation they are the same.

Ignoring the scaling (or assuming $1/\sqrt N$). We have

Apply one time: spectrum
Apply second time: time inverse original
Apply third time: conjugate spectrum
Apply 4th time: original signal

For inverse inverse FFT: reverse the order of the steps.


Can time-invariance be determined from a given a transfer function?

The question can only be answered if it is clarified what it is that the function $G(z)$ describes. If $G(z)$ is the system's transfer function then we're done immediately, because only linear time-invariant (LTI) systems can be characterized by a transfer function of that form.
If $G(z)$ is the transfer function, the output sequence $y[n]$ is given by the convolution of the input sequence $x[n]$ and the impulse response $g[n]$, which is the inverse $\mathcal{Z}$-transform of $G(z)$. Note that in order to compute the inverse $\mathcal{Z}$-transform you must know the region of convergence (ROC), or you must be given the additional information that the system is causal.
Having computed the impulse response $g[n]$, the output is given by the following convolution sum:
$$y[n]=\sum_{k=-\infty}^{\infty}g[k]x[n-k]\tag{1}$$
It's trivial to show that the response to a shifted input $x[n-n_0]$ is given by $y[n-n_0]$. Hence, any system whose input-output relation is described by a convolution $(1)$ is time-invariant (and, obviously, also linear).


How can i derive step response in terms of impulse response from the convolution sum?

In an LTI system, any linear operation on inputs, is directly imposed on the outputs, that is, if an LTI system responses to $\delta[n]$ as $h[n]$, then it responses to $u[n]=\sum_{k=-\infty}^{n}\delta[k]$ as $s[n]=\sum_{k=-\infty}^{n}h[k]$.

This is also clear from the convolution since $$s[n]=u[n]*h[n]=\sum_{k=-\infty}^\infty h[k]u[n-k]=\sum_{k=-\infty}^n h[k]$$


The impulse response of two filters in series is the convolution of each filters impulse response. You can multiply two polynomials by convolving their coefficients. Likewise you can factor or combine filters by factoring or multiplying their transfer functions described as polynomials.  For example, a fourth order filter can be factored into two 2nd order filters in series, or two 2nd order filters can be combined into one fourth order filter- it’s simply polynomial multiplication or factoring.

So in your case since you already factored the Butterworth filter into individual filters that you know how to implement, you would then simply cascade those filters.

As @MattL points out in the comments, you must pay careful attention to impedances when cascading analog filter sections. In the analog such structures are quite common with active filter design where it is quite simple to isolate the sections with appropriate design techniques, but not so common with filters constructed with passive components.


Sampling Theorem and Signals Explained to a Mathematician

"Highest frequency" is shorthand for the highest k/T for which a and b sub k are non-zero.

Considering f, one reason it might consist of the sum of multiple sinusoidal components is that this is a solution to many ordinary linear low-order differential equations.


Yes I and Q are in quadrature and therefore orthogonal so they are naturally isolated. For example you could consider QPSK as two independent BPSK channels transmitted on the same carrier frequency (or you could consider it having twice the spectral efficiency as BPSK since you can send twice as much data at once within the same bandwidth, and at the same SNR level).

To maintain isolation in the SDR you need to have good IQ balance or correct for IQ quadrature and amplitude error. Correcting for analog errors with software as depicted in the diagram below is easy to do: if the data content is equally distributed in the long term then the data can be used directly, and this would be a reason to "whiten" the actual data with scrambling prior to transmission. Otherwise training sequences could be used for such calibration of the receiver hardware- to the extent that is needed to be sufficient for audio purposes. Amplitude is adjusted by comparing the rms magnitudes of I and Q, and phase is adjusted by correlating I and Q (average of the product of I and Q over all samples) and scaling by the amplitude (to get phase).


Assume low-pass signals
throughout.

Since $X(f)$ is usually complex-valued, using the power spectrum $|X(f)|^2$
is probably a better idea, especially if you want to take square roots
etc. afterwards. Thus, $m_k$ is defined as
$$m_k = \int_{-\infty}^\infty f^k |X(f)|^2 \mathrm df.$$
Note in particular that $m_0$ is the power in the signal, and $m_1 = 0$
Now, the Gabor bandwidth $G$ of a signal is given by
$$G = \sqrt{\frac{\displaystyle \int_{-\infty}^\infty f^2 |X(f)|^2 \mathrm df}{\displaystyle\int_{-\infty}^\infty |X(f)|^2 \mathrm df}} = \sqrt{\frac{m_2}{m_0}}.$$
To put this in a slightly different perspective, $|X(f)|^2$ is a nonnegative function,
and the "area under the curve $|X(f)|^2$," viz. $m_0$, is the power in the signal.  Therefore, $|X(f)|^2/m_0$ is effectively a probability density function of
a zero-mean random variable whose variance is
$$\sigma^2 = \int_{-\infty}^\infty f^2 \frac{|X(f)|^2}{m_0} \mathrm df
= \frac{\displaystyle\int_{-\infty}^\infty f^2 |X(f)|^2 \mathrm df}{\displaystyle\int_{-\infty}^\infty |X(f)|^2 \mathrm df} = G^2$$.

A sinusoid of frequency $G$ Hz has $2G = 2\sqrt{\frac{m_2}{m_0}}$
zero crossings per second.  Since Mohammad is reading a legacy book,
it may well be doing all this in radian frequency $\omega$, and thus
if $G$ is the Gabor bandwidth in radians per second, we need to divide
by $2\pi$ giving
$$N_0 = \frac{1}{\pi} \sqrt{\frac{m_2}{m_0}} ~ \text{zero crossings per second.}$$

Turning to extrema, the derivative of $x(t)$ has Fourier transform
$j2\pi f X(f)$ and power spectrum $|2\pi f X(f)|^2$.  Its Gabor
bandwidth is
$$\begin{align*}G^\prime &= \sqrt{\frac{\displaystyle\int_{-\infty}^\infty f^2 |2\pi f X(f)|^2 \mathrm df}{\displaystyle\int_{-\infty}^\infty |2\pi f X(f)|^2 \mathrm df}}\\
&= \sqrt{\frac{\displaystyle\int_{-\infty}^\infty f^4 |X(f)|^2 \mathrm df}{\displaystyle\int_{-\infty}^\infty f^2 |X(f)|^2 \mathrm df}}\\
&= \sqrt{\frac{m_4}{m_2}}.
\end{align*}$$
Using the same arguments as before (two zero-crossings of the derivative
per period is the same as two extrema per period), radian versus
Hertzian frequency, we get
$$N_e = \frac{1}{\pi} \sqrt{\frac{m_4}{m_2}} ~ \text{extrema per second.}$$


Does MATLAB code retain its accuracy when is executed in dsp kit?

Yes, C\C++ code generated from Matlab/Simulink code should behave the same as the original code.
However there are a few caveats, here are 3 problems I've experienced in real-life.
1 - Simulink\Matlab performs an implicit reset just before starting the simulation/script. Make sure you perform this reset in your C\C++ code, otherwise your results won't be the same and won't be comparable. Simulink coder generates a reset function for you. I've seen a corner case where the reset function would not reset all global variables, so be careful.
2 - Variable types. Matlab/Simulink typically use double-precision floating point variables. If your code uses doubles then great. The behaviour should be the same. However, for performance reasons, you might choose to use single-precision floating point variable aka floats. Then the behaviour could change slightly. There are ways to mitigates the effect of quantization, some algorithms are really sensitive to precision, others less so. Even if you use doubles there are some difference between platforms. So if you use doubles on an ARM µProcessor, the behaviour could be slightly different in corner cases.
3 - Integer saturation. I made the mistake in the past of using a 32-bit integer in Matlab as a free-running counter. It worked well, but after like 10 days the code would freeze. The counter would be stuck at 2^32-1 instead of wrapping back to 0. Obviously, you won't see that behaviour in simulation. So be careful, test your code thoroughly.


Time-bandwidth product and uncertainty on harmonics

If you a-priori assume that $k \cdot f$ is a stationary harmonic of periodic waveform of frequency $f$, that implies that you don't have any separable time locality of that harmonic.  But any tighter bandwidth of a frequency estimation might come from the a-priori assumptions regarding the stationarity of the harmonics, not from the data.


Reverse Engineering Bessel LPF (sample frequency needed)

I don't think Deve's answer is correct. The coefficient determine the transfer function in normalized frequency. You can determine at what relative frequency the transfer function hits -3dB. Let's say it's 0.01. Since you  know this is corresponds to 10 Hz you can infer that the sample rate is indeed 1000Hz.
This code is difficult to read and there isn't a single line of comment in it, nor is there any explanation of what the state variables so it's difficult exactly which number is what coefficient (b0,b1,b2, a0,a1, a2 etc) but it is in principle possible to infer the sample rate.


What is convolution of two sine waves (tones)?

Convolution in the time domain is equivalent to multiplication in the frequency domain.

If you window two sinusoids in the time domain to get finite length waveforms, and the two sinusoids are exactly integer periodic in the window width, then the DFT will be impulses.  If the frequencies are different, the impulses in the two DFT will be disjoint, and the multiplication of the two spectrum will result in zero output.

If one or both of the two sinusoids are not exactly integer periodic in your window width, then the FT will result in a Sinc function.  Then the convolution will be equivalent to multiplying a Sinc function against an impulse or another Sinc function, which will produce a non-zero result.


I seems to me like wavelets can offer a neat solution to your problem for noise filtering. Just incase you are not aware a wavelet transform uses two filters (also called mother and father wavelet). These two filters result in a sub-signal (also called sub-band) which is of half the size of the original signal. What I mean to say is say your signal length is 100. Upon wavelet transform you get 2 signals each of size 50. One of these signal is a approximation of original signal and the another one is detail signal. Usually the noise component is present in the detail section.

So as to get rid of the noise, what is usually done is to take a wavelet transform. Implementation of it is available in matlab (dwt() ). This will result in two outcomes cA and cD.

[cA,cD] = dwt(X,'haar')

Then I would suggest you to do usual noise removal algorithms on cD and with cA and cD' (ie. filtered cD) reconstruct the original signal.

cD_1 = someProcessing( cD )

X = idwt(cA, cD_1,'haar')

Please have a look at the following links to get more information on noise filtering with wavelets
1. PhD Thesis titled "Signal Recovery and Noise Reduction with Wavelets" download from HERE.


How to Plot Summation of Shifted impulses in MATLAB?

According to the features of the unit impulse function $\delta(t)$; at t = 0, the function goes toward infinity. Actually, the unit impulse function is meaningful with its area which is 1.
Additionally, MATLAB don’t show the famous arrow notation of the unit impulse function. Instead, it shows an invisible line which can be noticed where the function increases into infinity, i.e. the time value on the time axis where the function increases into infinity is blank. This can be seen from the figure.
On the other hand, the unit sample function $\delta[n]$ value is defined with the area of the unit impulse function which is 1.
We can obtain h[n] function by sampling h(t) function. But, what about those inifinities?
So, let’s do some gamble!
If we code h(t) function in MATLAB, it will actually be a row vector. In this row vector, we can change those infinities with ones in order to sample the function correctly. This process will require a for loop and an if statement.
The MATLAB source code for the complete process is given below:
t = -10: 1: 10; % Time range

Ts = 1; % Sampling period

syms k; % Symbol assingment for summation

function_1 = (t + 1) .* symsum(dirac(t - k), k, 1, 5); % Summation

subplot(2, 1, 1);
plot(t, function_1, "r", "LineWidth", 3);
xlabel("Time (s)");
ylabel("h(t)");
title("h(t) Time Function");

for a = 1: length(function_1) - 1 % For loop-if statement duo for replacing infinity results with unit sample function
    if function_1(a) == Inf
       function_1(a) = 1;
    end
end

function_2 = (t + 1) .* function_1; % Complete function to be sampled

n = t / Ts; % Sample range

subplot(2, 1, 2);
stem(n, function_2, "ro", "LineWidth", 3);
xlabel("Samples [n]");
ylabel("h[n]");
title("h[n] Sample Function");

According to the figure, sampling process has gone well and it matches with analysis that I’ve done on the paper.

Actually, there may be other approaches for the replacements of the infinities with ones but this one does the trick well.


Two type of calculating R in steepest-descent modeling algorithm

So let's have a look at your calculations:

x1(:,i)=x(N-(i-1):end-(i-1));

Let's say N = 10 and sample = 100 so that end = 110.

For i=1 that means we sum from 10 to 110.

For i=2 that means we sum from 9 to 109.

etc.

Comparing that with:

R(k)=h(k:m)*h(1:m-k+1)';

where m is the same as sample:

For k=1 that means we sum from 1 to 100.

For k=2 that means we sum from 2 to 99.

I think all that's happening is your first calculation is assuming x is zero-padded and the second one is assuming that it's not.

That there are differences between the two is hardly surprising, given that they are summing different (possibly random) values.


Conditions for which the Hilbert transform returns a correct phase

A single instantaneous phase estimate may or may not make any sense if there is more than one frequency peak in the signal's local spectrum.  So, to get a better single frequency and phase estimate, you may want to first high-pass filter your signal to remove the spectrum of any D.C. bias or slow trend, and also band-pass filter the signal to remove any outside spectral peaks due to noise or interfering signals.

So my off-hand guess for a condition is that the local spectrum has a single clear frequency peak with a monotonicly decreasing spectral magnitude away from that one peak.  All the way to zero at DC for real signals.


FFT of audio, packet sampling and others

As some comments have pointed out, there are a couple of practical ways to find out the sampling rate:

Read the sampling rate information from the WAV file header
Look for mains hum in the recorded audio, which you'd see as a peak in the FFT response in the frequency bin corresponding to 60 Hz (or maybe 50 Hz).

But it sounds like you're looking for a mathematical property of the FFT that will allow you to back out the sampling frequency. Such a thing doesn't exist.

The only inputs to the FFT are the original sampled data. And this sample data itself has no information about the real sample rate contained in it. So there's no possibility that the spectral result of the FFT can tell you anything about the sample rate.

For example,  there is no way to tell the difference between a 10 Hz signal sampled at 100 Sa/s from a 320 MHz signal sampled at 3.2 GSa/s.

In order to know the sample rate you must record it as meta-data, for example in the header of a WAV file.

  If I have a .wav file with N datapoints then split into ten N/10 segments and then fft on each segment. But this then means that my sampling frequency is essentially reduced to N/10 Fs?

No, the sample rate is still Fs. But the sample duration has been reduced to N/(10 * Fs).

This means that the spectrum will still cover the same frequency band (0 to Fs*(N-1)/N), but it will cover it in coarser steps. The bins of the DFT will be spaced by 10*Fs/N instead of Fs/N.


Non Coherent BPSK with no Channel State Information at the receiver

In any PSK modulation, the information is only in the phase.
An unknown $h$ randomly rotates your phase. You know nothing about the phase that $x$ had when you observe $y$. Hence, you have zero mutual information.
One could write that down as a very short formal proof based on conditional entropies, but I don't think that'll help you; the more intuitive understanding "without knowing $h$, the output phase is fully random" is sufficient as proof.


First, you should note that the FFT refers to an algorithm for the fast computation of the DFT taking advantage of its symmetry.

The signals you have are:

where $\mathbf{a}_y$ and $\mathbf{a}_x$ can both have negative a positive values depending on the direction along the $x$ and $y$-axis. These signals are the output of your accelerometers.

What you call "a piece-wise resultant acceleration", calculated like
$$r(i) = \sqrt{x(i)^2+y(i)^2},$$ is just $\rVert \mathbf a\lVert$, the magnitude of a vector quantity $\mathbf{a}$, which doesn't show the direction.

From NASA webpage, about vector components:

  Why go to all this trouble? Because, in aerospace, we are often
  dealing with forces and forces are vectors. Breaking a single vector
  force into several components allows us to study the resulting motion
  much more easily.

In conclusion, performing a DFT on your vector quantities $\mathbf{a}_x$ and $\mathbf{a}_y$ will give you frequency information both for about magnitude and the direction together, while performing a DFT on your scalar vector $\rVert \mathbf a\lVert$ will just give you the frequency domain of the magnitude, so yes, you are loosing information here.


Why are binary signals drawn continuous in time despite being a digital signal

You are right and wrong!

If you would like to treat a binary signal as a discrete-time sequence having values of $0$s and $1$s, then you would plot them as all other discrete-time sequences, using the lollipop display.

However, for those electronic circuits that implement binary logic operators using physical devices (gates etc.), the associated voltage and current waveforms are analog in nature. What's binary there is the information that the voltage or current waveform carries. But the physical existance is analog. So if you want to plot them, then you would do it just like a logic analyser (oscilloscope) would draw them on their phosphorus crt screen...


Multi-path is a property of the communication channel (mostly the wireless channel). The received signal $r(t)$ is the sum of delayed and weighted versions of the transmit signal $x(t)$:
$$
r(t) = \sum_{i} a_i x(t-\tau_i)
$$

Transmit diversity is a property of the transmitter. It describes a technique using multiple transmit antennas, thus creating multiple communication channels that are different because the antennas are placed at different locations in space. The idea is to transmit the same information at the same time at different antennas.

For example, let's assume one receive antenna and two transmit antennas (2x1 multiple-input single-output (MISO)) the received signal is a superposition of both multi-path channels:
$$
r(t) = \sum_{i} a_i x(t-\tau_i)+\sum_{j} b_j x(t-\tau_j)
$$
The hope is that if transmission on one channel fails you will still receive a useful signal through the other channel thereby increasing the robustness of the system.

There is also the concept of receive diversity where multiple receive antennas are used (SIMO). And of course you can use multiple antennas on both receiver and transmitter side (MIMO).

When using the scheme of transmitting the same data at the same time, an FIR filter can be used no matter if transmit diversity is applied or not, because to the receiver it will just appear like an ordinary multi-path channel. Advanced equalizer structures include the fractionally spaced equalizer and the decision-feedback equalizer. If you apply some space-time coding a the transmitter an according space-time decoder is required at the receiver.


How to do this inverse Z-transform? Is this correct?

A rational transfer function

$$H(z)=\frac{Y(z)}{X(z)}=\frac{\sum_{n=0}^{N}b[n]z^{-n}}{1+\sum_{n=1}^{N}a[n]z^{-n}}\tag{1}$$

corresponds to the following difference equation:

$$y[n]=b[0]x[n]+b[1]x[n-1]+\ldots +b[N]x[n-N]-\\-a[1]y[n-1]-\ldots -a[N]y[n-N]\tag{2}$$

So the current output sample $y[n]$ can be computed from the current input sample $x[n]$ and $N$ past input and output samples. $N$ is the order of the difference equation $(2)$.

In your example you have

$$H(z)=\frac{Y(z)}{X(z)}=-c\frac{1-z^{-1}}{1-cz^{-1}},\qquad c=\frac{1}{1+\frac{2RT}{m}}\tag{3}$$

and, consequently,

$$y[n]=c\big(-x[n]+x[n-1]+y[n-1]\big)\tag{4}$$


Why Is the DFT of a Signal Symmetric About Its Central Bin?

First of all it is only conjugate symmetric if the input signal is Real Signal.

The reason for that is the DFT is built by Complex Exponential which can be decomposed into Sine and Cosine.
Since the Sine governs the Complex part and is $ N $ periodic one could see it creates conjugate symmetry of the DFT.

You just need to follow the definition of the DFT.


500 kHz (watch the capitalization: Khz is wrong!) is a very low sample rate for the N210. In fact, the N210 has to decimate by a factor of 200 from its ADC's sampling rate, and that might introduce unexpected roll-off and aliasing, because that over-stresses the CIC decimator (if this doesn't make too much sense to you: don't worry, it's an implementation detail. Just: rather go for 1 or 2 MS/s sampling rate and do the rest of the decimation in software).
Also make sure you're actually seeing the signal you think you're seeing: your USRP AND your LoRa transmitter might both have frequency offsets. But: 15 MHz difference at 915 MHz? That would be an error of 1/61 ; the USRP has a frequency accuracy of ca 10 ppm (that's 1/100000), your LoRa transmitter maybe of 50 ppm. No way you're observing a 915 MHz signal at 930 MHz. Either this is something else, or you're misinterpreting something and your transmitter isn't trying to transmit at 915 MHz!
It's hard to tell from your slightly unusual choice of color scheme whether your spectral broadening happens from an overdrive of the receive amplifier. Reduce gain – if it gets better from that, that's happening.
Also make sure you're not actually seeing what you want to see – after all, LoRa is chirp-spread spectrum, and you're sampling these chirps at more or less irregular intervals. How does your spectrum change when you change the FFT size?


Normally it is 94 dB of Sound Pressure Level, where 0 dB is 20 µPa of acoustic pressure. This is the SPL you will find in most of the acoustic calibrators for measurement microphones.


The DTFT is independent of Nyquist sampling. The DTFT is applied to a sequence of numbers regardless of any analog waveform that may be sampled. Economic sequences are often inherently discrete time sequences.  The DTFT applies to aliased sequences as well.

The summation is in general over $n$ from $-\infty$ to $\infty$.  The summation limits can be truncated to the nonzero range of $x[n]$.

The nominal convention as presented in Oppenheim and Schaefer, which they introduce prior to their chapter on Nyquist sampling, is that the DTFT is a continuous valued, periodic function defined uniquely on $-\pi$ to $\pi$ which has a length of $2\pi$.

For real valued sequences, the interval $-\pi$ to $0$ is conjugate mirrored of the spectrum over $0$ to $\pi$.

In general for complex valued sequences, the mirrored conjugate symmetry can’t be assumed.

$0$ to $\pi$ is only one half of a length $2\pi$ period.

edit

$-\pi$ to 0 has an image at $\pi$ to $2\pi$ and $3\pi$ to $4\pi$ and ......

0 to $\pi$ has an image from $2\pi$ to $3\pi$ and $4\pi$ to $5\pi$ and ....

similarly for the negative frequencies.


The MFCC summary you link seems to leave out the typical windowing function applied before each FFT.  Segmenting longer data into shorter finite length FFT inputs does an implicit rectangular windowing, which causes the energy of the frequency of any signal not exactly integer periodic in the FFT length to be "spattered" into other distant frequencies of the FFT result.  A non-rectangular window function is applied to help reduce these rectangular window artifacts (sometimes called spectral "leakage").

A typical window function, such as Von Hann or Hamming, can be quite lossy at both edges of each window.  Overlapping windows adds a window (or more) centered closer to where other windows are lossy, thus helping preserve information about the original signal that would otherwise be lost or degraded.

Another reason is for time resolution.  The length of each FFT window might be chosen to give sufficient frequency resolution, however, due to the time-space resolution trade-off of FFT analysis, the resulting time resolution may not be adequate, and thus miss fast audio transients not being centered or isolated in one window.  Shorter FFT windows may not provide sufficient frequency resolution.  A greater density of FFT windows (using overlap instead of shorter windows) improves the capture and isolation of these shorter time domain events.  This greater density (more windows per unit time) can be accomplished by overlapping the windows.

If you don't overlap, then events near or at the boundaries will be severely degraded, and the probability of transients (such a consonants) being between windows, or not being isolated in a single window, is increased, thus reducing how well the MFCC can categorize different inputs.


How do I differentiate between signals from an unknown program?

Just a few pointers, since a complete answer would be way too long.

The link you give for the key fob lists its FCC ID: CWTWB1U331. Using a search engine, you can find its entry in the FCC database: https://fccid.io/CWTWB1U331

The FCC database is a goldmine because it lists many of the operational parameters of the device. For example, you can verify that the frequency range you're looking at is correct.

Among other things, you can learn that the modulation is ASK, and that the line encoding is Rolling/Manchester. In particular, the "Test Report" pdf has a couple of extremely helpful figures, which will help you narrow down the transmitted frequency and the line encoding format.

For example, from Figure 2 in that document, it seems like there are 10 bits in 5 milliseconds, for a data rate of 2000 bits per second.

Another thing to note is that essentially all key fobs will transmit a different word every time the key is pressed, so the job of identifying the key will imply finding each key's patterns, which will not be easy -- otherwise, car theft would be a piece of cake.

Finally, there are several resources online for reverse-engineering wireless signals, such as remote controls and garage door openers. I recommend that you look for those.


FMCW SAR chirp return from target modelling

I am not sure you really understand the processing you do so I will try to explain what you measure.
When you use the returned chirp mixed against the reference chirp you get an instantaneous phase measurement. When you do the FFT the frequency domain represents distance - since the phase differences are changing faster as the distance is greater.
The inverse FFT of the Distance will give you (assuming high signal to noise ratio) a periodic signal that is the function of the distance you measured.
In the graph you show, it is not clear what you show since the phase is not presented.


Getting the $\mathcal Z$-transform of the moving average summation

You can use the linearity property of the $\mathcal Z$-transform:

$$a_1x_1[n] + a_2x_2[n] \overset{z}\longleftrightarrow a_1X_1(z) + a_2X_2(z)$$

together with the time shifting property:

$$x[n-k] \overset{z}\longleftrightarrow z^{-k}X(z)$$

Expand the summation like below and apply the two properties.

$$\sum_{k=0}^M b_k  x[n-k]=b_0x[n]+b_1x[n-1]+\ldots+b_Mx[n-M]$$


When doing a Hilbert-transformation, why not simply multiplying by an exponential?

In contrast to Jason R's answer I claim that the Hilbert transform is a phase shift by $-\pi/2$ for real-valued signals. By definition, a phase shifter shifts the phase of a sinusoidal signal by some given phase $\phi$:

$$x(t)=\cos(\omega_0t)\quad\Longrightarrow\quad y(t)=\cos(\omega_0t+\phi)\tag{1}$$

Since

$$\cos(\omega_0t)=\frac{e^{j\omega_0t}+e^{-j\omega_0t}}{2}\tag{2}$$

and

$$\cos(\omega_0t+\phi)=\frac{e^{j\omega_0t+\phi}+e^{-j\omega_0t-\phi}}{2}\tag{3}$$

it is clear that in order for the signal $x(t)$ to experience a phase shift by $\phi$, the positive frequencies need to be shifted by $\phi$ and the negative frequencies by $-\phi$. For $\phi=-\pi/2$ that's exactly what a Hilbert transformer does.

A constant transfer function like $H(f)=e^{-j\pi/2}=-j$ doesn't really do much apart from scaling the input signal, in this case by the imaginary unit. So the signal $x(t)=\cos(\omega_0t)$ is "transformed" into $y(t)=-j\cos(\omega_0t)$. Of course, for the special case of the complex-valued input signal $x(t)=e^{j\omega_0t}$ multiplication by $-j$ achieves a phase shift by $-\pi/2$, but usually we want to shift the phase of real-valued signals, and in that case we need a frequency-dependent transfer function.


This may be a good candidate for a two step approach:

Step 1 would be a pretty coarse running cross correlation with a threshold detector to identify what parts in the signals are matching.

Step 2 would than determined the actual delay. There are various ways to get sub sample resolution:

Upsample to desired resolution and then cross correlate
Short term fourier transform and match a linear phase difference with a weighted least squares approach
Delayed lock loop with a fractional delay filter
Match the the two signals in with an adaptive filter. Then calculate the fourier transform of the filter impulse response and calculate the delay through a linear phase match


Beyond "mid-side": decompose an audio stereo recording (coincident microphones) into Mid + 45° + 90°

There is a generic class of algorithms called "upmixers" that indeed de-compose a stereo signal into multiple direction and/or it's constituent directional channels. They typically work by chopping the signal into frames, evaluating correlation/phase/directional differences in multiple frequency bands, steering the origial signal on rules based on the analysis and then splicing them all together again.
It's tricky business and typically involves quite bit of trade off between the effectiveness of the spatial decomposition and spectral and temporal artifacts.
However, I'm not convinced that this would help in your case. You have a far field recording done in a reverberant space: the energy of your direct sound component (which is the only thing localizable) is probably tiny as compared to the energy of the early reflections and reverb. Trying to pull out directional components out of the reverberant field is hard and chances are they are decorrelated anyway, so applying time delays to them is unlikely to make much of a difference.
If you just want a little more space in the recording, you could try to boost the difference channel a bit with a high boost shelf filter. No boost for low frequencies, some for the highs. You can play around with the cutoff frequency, but I would vary the cutoff from 500Hz to 2kHz and do try boost of 2-3 dB or thereabouts.
If you are aiming for "more natural over headphones", your best shot is to match the inter-channel correlation as a function of frequency to that of a human (or dummy) head in a diffuse field. At low frequencies the two ear signals are highly correlated, at high frequencies they tend be very uncorrelated and they are partially correlated in between. You can analyze the correlation of your recording and compare to a binaural target curve and then increase/decrease correlation at specific frequencies as needed. There are various ways of doing that, but the easiest way is to adjust the sum/difference balance: more sum = more correlation, more difference = less correlation.


Though acoustic speech features (MFCCs, etc) & lexical words are correlated in daily language, traditionally in speech recognition software, acoustic and linguistic information are assumed independent.

You can see this in popular open-source speech recognition packages like Sphinx3/4 (Java/C using HMMs vs DTW) or Kaldi (C++/neural nets), one must train both an acoustic model & a language model. For the English language, there are a few off-the-shelf models available, but in most domain-specific cases, it is best to train a new one or adapt an existing model. I'm not sure about the state of acoustic/language models for the Hindi language, but here is a paper you may have seen.

I recommend looking for existing Hindi models to use within Sphinx or Kaldi.


Pulse shaping with root-raised cosine filter for a non-Gaussian channel

There are many related questions on this site; make sure to review them. To summarize:

Correct.
The roll-off factor $\beta \in [0,1]$. If $R_p$ is the pulse (or symbol) rate, then the pulse bandwidth is $B=(1+\beta)R_p/2$. When $\beta=0$, the pulse is a sinc and the bandwidth is the smallest possible. For $\beta>0$, the factor $1+\beta$ is the "excess" (not "sacrificial") bandwidth.
The AWGN channel does not distort the input signal, only white Gaussian noise is added to it. The density of the noise is Gaussian (normal), not the channels' impulse response. The noise, being a signal, cannot be linear or non-linear (only systems have this property).
RRC and other Nyquist pulses can avoid ISI in many different channels. For example, in  a frequency selective channel, after linear equalization the received signal has no ISI. RRC pulses remain the most popular in practice, in all kinds of channels.


Comparision between SSIM and MAD Image Quality Assessment Algorithms

A natural practice to convert, in a monotonous way, data in $[0\,\infty)$ to $[0\,1)$ is an increasing function, like those called sigmoid functions. They behave like $x\mapsto x$ close to $0$, and flatten to $1$ when $x\to \infty$. Classical examples are:

$$  x\mapsto \frac{x}{1+x} $$
$$  x\mapsto \frac{x}{\sqrt{1+x^2}} $$
$$  x\mapsto \frac{2}{\pi}\arctan{x} $$

There are many more versions, since they are used a lot in statsitics, articial neural networks, etc.


Calculating original signal from Discerete Fourier Transform

At first glance your first paragraph seemed wildly confusing to me. (The notion of "calculating" an equation makes no sense.) But now I see what your professor is asking of you. He wants you to determine the coefficients of an equation. He gave you a continuous $f(t)$ equation of a signal containing a DC component (1cos(0$\pi$$t$/4) and three sinusoids. And the 2nd and 3rd sinusoids have frequencies that are integer multiples of the 1st sinusoid.

So, your prof wants you to select the values and length ($N$) of a $t$ time sequence that will make those sinusoids reside exactly at an $N$-point DFT's bin centers. When that happens the bin centers' frequencies will give you the frequencies of the three sinusoids in your $f(t)$ equation. The peak values of the non-zero DFT spectral components (DC plus the three sinusoids) will enable you to establish the [1, 0.75, 0.2, 0.5] peak amplitude coefficients in your $f(t)$ equation.  Your prof wants you to learn how a sinewave whose peak amplitude is $A$ produces a DFT spectral component magnitude value of $M$ given the $N$ length of your DFT.  Your job is to see how $A$, $M$, and $N$ are related to each other. This is a good homework problem.

Your $t$ = [0,1,2,3] sequence's period is one which means your sample rate is 1 sample/second. And as Ollie said, that's too slow. Set your $t$ sequence to [0, 0.5, 1, 1.5, 2, ..., 6.5, 7, 7.5], compute a 16-point DFT, compute the spectral magnitudes, and see what happens
.


I managed to find a solution that is working quite well. Im not fully optimized yet but I think its just tweaking filters, as ill explain later on.

Initially I tried to set up opencv but it was very time consuming and a steep learning curve but it did give me an idea. The key to my problem is really detecting the characters within the image and ignoring the background, which was basically just noise. OCR was designed exactly for this purpose.

I found the free library tesseract (https://github.com/ldiqual/tesseract-ios-lib) easy to use and with plenty of customizability. At first the results were very random but applying sharpening and monochromatic filter and a color invert worked well to clean up the text. Next a marked out a target area on the ui and used that to cut out the rectangle of image to process. The speed of processing is slow on large images and this cut it dramatically. The OCR filter allowed me to restrict allowable characters and as the plaques follow a standard configuration this narrowed down the accuracy.

So far its been successful with the grey background plaques but I havent found the correct filter for the red and white editions. My goal will be to add color detection and remove the need to feed in the data type.


Fourier transform 4 times = original function (2D and higher)

The n-dimensional (discrete) Fourier transform is separable, that means the dimensions can be treated independently. So the property that the FT is a fourth root of the identity applies to n-dimensional transforms too. And the property is exact, so apart from a tiny numerical round off error you won't see any deviations.


Time scaled version of discrete time signals

Now if i plot u[n/3] i know there will be two zeros added between each previous samples(samples before scaling).

No.
You are confusing up-sampling with time scaling. The time scaled version of a unit step is a unit step. You would only insert 2 zeros if you up-sample by a factor of 3.


What is done to minimize distortion due to the hold operation?

What you are describing is the distortion introduced by an ideal digital-to-analog converter (DAC) in the analog domain.  Two things are typically done to reduce this distortion:

Analog filtering
Oversampling

As you note in your question, the distortion is modeled, in the frequency domain, by a sinc rolloff.  Increasing the sample rate before converting the digital signal to analog reduces $T$/2, which make the rolloff slower, or "flatter".

Analog filtering can range from simple low-pass filters to compensation filters that counteract the sinc rolloff.  These filters do not, I believe, fix the phase issues, but they can greatly improve the frequency amplitude response.

CD players are a good example of these techniques.  Audio CD's are sampled at 44.1 kHz, but before the DAC the CD players will typically increase the sample rate to a few hundred kHz.  I'm not sure how sophisticated their analog filtering is.  My guess would be that the cheap stuff does simple low-pass filtering and the expensive ones do sinc rolloff compensation.


Running Integral of sine and cosine functions

I would not teach integrals this way.

Issues might arise form the definition of terms. In standard calculus, a (note: "a", not "the) primitive function (also called antiderivative or indefinite integral) of a continuous function $f$ is a kind of converse of the concept of derivation. SO:

  If $F$ is "a" differentiable function  whose derivative is equal to
  $f$, then $F$ is a primitive (function) of $f$. There exists many such
  functions, up to a constant $c$, because: $(F(t)+c)' = (F(t))' =f(t)$.

Here, you can check that $t\mapsto \sin(t)+ c$ is differentiable, with derivative $t\mapsto \cos(t)$.

Now (if I remember well), suppose that $f(x)$ is a continuous function, AND that it is  "integrable" over $\mathbb{R}$, then $F(t) = \int_{-\infty}^t f(\tau)d\tau\;$ exists, and is a differentiable function. THEN, it is also a primitive of $f(t)$ (among others).

But here, it is not possible, namely because $t\mapsto \cos(t)$ is not "integrable" over $\mathbb{R}$, and thus the integral expression makes no sense.

I would say that one should not define "integrals" the way you wrote it, if one wants mathematical correctness. However, for causal signals, in a lousy engineer perspective, it may be intuitive. But the cosine is not causal.

[Note: there are more general definitions of primitives]


How to choose the right segment length for a given signal?

The DFT in of itself is not a consistent estimator, or in other words, the variance of your transform doesn't decrease as the transform length increases.

Another way to think of it is that the DFT is a linear invertible transform, so whatever "information" is contained in the time domain is contained in the frequency domain.  The segment that is 256 points is expected to be less smooth than a 128 point segment. If the longer segment has more variability, it's DFT will also have more variability.

To have a consistent estimator one can use the technique attributed to Welch, which means that you average a a set of DFT segments and that generally means averaging bin magnitudes (or magnitude squares).

Picking a correct segment length is related to how many segments you will average.  The number of averages is often tied to the amount of data you can average, as well as a target bin variance.

If the data has a spectrum that varies, like a speech signal, one can average too much and obscure the dynamic nature of the data.

Window functions have a large influence on smoothness, and are implicitly related to segment length.

Without knowing what your target spectrum's purpose is, any notion of  "correct length" will be subjective.   Subjectivity isn't necessarily bad.


Combining channels into single FFT left/right vs mid/side

Whether FFT(L+R) will work depends on how the stereo recording was mic'd and mixed.

If just 2 separated stereo microphones were used to record a group of spatially separated sound sources (say an orchestra), then there will be location and frequency dependent phase differences between the 2 channels than could cause cancellations or nulls in a simple (L+R) mono mix.  If a lot of individually mic'd tracks were mixed down into 2 channels, this might be less of a problem.


Difference between diversity combining and multiplexing

diversity combining such as maximal ratio combining, the same message is sent over different frequencies

That is not inherently true.

Diversity combining makes use of different channels (not in the "TV channel" sense, but in the "everything that happens to a signal" sense).

Those channels might be created by over different frequencies all carrying the same signal, but often/usually, at the receiver, you use multiple e.g. antennas to observe the same signal.

The different signals are then weighted (often with a complex factor to reflect both phase and amplitude differences), and then added up, is what most commonly is referred to as diversity combining.

You could actually do the same with a single antenna and multiple identical signals mixed to different frequencies. However, that has little to no diversity gain: You actually need $N$ times the bandwidth of the original signal, so you could as well just transmit more data, add more redundancy etc.


Evaluating different inputs into a given system

The output will be equal to the input $(y[n]=x[n])$ if no aliasing occurs as a result of downsampling. Which means $X(e^{j\omega})=0$ for $\frac{\pi}{4} \leq |\omega| \leq \pi$.
a) $X(e^{j\omega})$ has impulses at $\omega = \pm \frac{2\pi}{7}$ $=>$ aliasing exists $=>$ $y[n]\neq x[n]$
b) $X(e^{j\omega})$ has impulses at $\omega = \pm \frac{2\pi}{9}$ $=>$ no aliasing $=>$ $y[n]= x[n]$
c) $X(e^{j\omega})$ has impulses at $\omega = \pm \frac{\pi}{8}$ $=>$ no aliasing $=>$ $y[n]= x[n]$
d) You need to sketch $X(e^{j\omega})$ to see if aliasing exists  hint: the impulse response of your filter is $h[n]=\frac{sin(\pi n/4)}{\pi n}$


Purpose of eigenspace of covariance matrix of a blob?

The eigenvectors will always be orthogonal. There is no guarantee that the eigenvectors will give you the bounding box with minimal area.


Prove a property using shift theorem and duality

If $X(f)=\mathcal{F}\{x(t)\}$ you have, in analogy with Eq. $(7)$,
$$x(f)\Longleftrightarrow X(-t)\tag{1}$$
where we assume that functions with the independent variable $f$ are Fourier transforms of the corresponding functions with independent variable $t$.
Applying $(1)$ to
$$x(t-b)\Longleftrightarrow e^{-2\pi jbf} X(f)\tag{2}$$
gives
$$x(f-b)\Longleftrightarrow e^{2\pi jbt} X(-t)\tag{3}$$
According to $(1)$, $X(-t)$ is the inverse Fourier transform of $x(f)$, so with $g(t)=X(-t)$ and $G(f)=\mathcal{F}\{g(t)\}=x(f)$ you obtain the desired result:
$$G(f-b)\Longleftrightarrow e^{2\pi jbt} g(t)\tag{4}$$


What's the advantage/disadvantage of oversampling followed by decimation, verses sampling at the correct rate to begin with?

I may have originally misinterpreted your question so I'll change my answer here.

As to why we don't typically just sample at our a decimation rate (i.e. the Nyquist rate for the bandwidth of interest). It really comes down to the practicality:

If you had a 200hz signal of interest, yes you could sample at 400hz assuming there was absolute no noise in your system. However, there will be noise and you would need an anti-alias filter to block all frequencies greater than 200hz. The issue with this is that many anti-aliasing filters are slow roll-off because they're cheaper and easier to implement. Higher order analog filters are simply more expensive and difficult to implement (they're more sensitive to component tolerances).

Instead, it would make more sense for us to over-sample our signal such that we could use a lower order anti-alias filter.

Let's say your signal of interest is still 200hz and you decide to sample at 10khz this time. You put in an anti-alias filter at say 5kHz with a slow roll off. This would prevent aliasing frequencies of getting into your digital signal while being able to use a cheap/easy anti-alias filter. You would then perform your decimation operation. Note that the decimation operation includes a digital low-pass filter such that aliasing doesn't occur in the down-sampling process.

Summary: the further we separate our bandwidth of interest from the Nyquist rate, we can use lower-order anti-alias filters which are cheaper to implement.


That's perfectly normal and expected. You are removing harmonics and as a result you change the time domain wave form. There is no guarantee or expectations that the maximum amplitude will be maintained.
In fact, a square wave has the lowest possible Crest factor (ratio of peak to RMS) and ANY change will increase the Crest factor. In order to not clip you need to remove a substantial amount of power (which your lowpass probably does).


Gonzales Digital Image Processing and the Laplacian Operator

Isotropic operator is an operator that has no bias to any particular direction in the image.
If you see the below Laplacian mask, Edge detection is taking place in every direction.

If you see the below image, three filter kernels are applied. The first kernel detects edges in x direction, the second in y direction and the third is Laplacian kernel detecting edges in all directions, hence it is isotropic.


Accelerometer - coordinate system transformation

It will call your function each time any motion event gets fired. This can affect the timing of css animations responding to sensor input.  Try getting the orientation values and updating your transformations from within an interval.


Decibels conversion into amplitude ratio

You will want to work with what Wikipedia calls "field quantities."

Your engineer is telling you to work in steps of .1dB, which would be an amplitude ratio of .1=20log10(x).  To get x you do 10^(.1/20) = 1.01158

From x dB to ratio (amplitude):  10^(x/20)

From x ratio (amplitude) to dB:  20Log10(x)

If you can calculate the ratio, then you should be able to calculate the correction needed, round to the next .1dB step, and apply it in one step.


How do I filter audio with modified spectrogram?

Modify the spectrogram and then use scipy.signal.istft

scipy.signal.istft(Zxx, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, input_onesided=True, boundary=True, time_axis=-1, freq_axis=-2)
Perform the inverse Short Time Fourier transform (iSTFT).


Matched filter - do I need to know the signal I am trying to find?

As stated by others, the signal $s[n]$ has to be known to calculate the corresponding matched filter.
An alternative is to use the power of deep learning to emulate a matched filter response as shown in this paper.


Sum of Squared Intensity Difference for block matching with ORB features, c++

Your results look reasonable. However, you're likely to get a somewhat unreliable result: the good matches are located in a small portion of the image, while a large spread would help in getting overall more accurate estimates (but that's probably due to your images).
Experimenting will work here! Otherwise, you can compute the expectation of the SSD between 2 blocks that are a rightful match under the noise level of your images.

Eventually, bad matches should be eliminated by the RANSAC part, but you may have to do a trade-off between getting many but sometimes wrong matches (loose SSD threshold) or stricter ones (at the risk of not having enough matches o estimate the pose).


Visual Tracking / Object Tracking: Estimate the Direction the Whale Is Moving

I wouldn't rely on edges as they are not connected and foam/waves also generate them.

I would rely on the contrast between the colors of skin, water and foam, to get a solid blob.

Possibly a two step process: discard the more saturated pixels (the sea), then among the achromatic keep the dark ones. After experimenting I saw that erosion can be useful to get rid of unwanted residues.

In the end you should get a big blob, and determine its dominant orientation (equivalent ellipse of inertia, Feret diameter, most distant points...)

Unsaturated pixels:

Dark pixels:

Combination:

Eroded:

Blobs and ellipse of inertia:

As there is a risk of blob fragmentation, you can consider merging the largest blobs.


Image Processing - find the centre of the estimated circle of eye

Unfortunately, I dont have access to a MATLAB image processing toolbox. I've coded a solution in python, showing the idea of how to solve it:

import skimage
import skimage.io
import scipy.ndimage as ndi
import skimage.morphology as m
import skimage.measure as measure

plt.figure(figsize=(8, 12))
I = skimage.io.imread("https://i.stack.imgur.com/eCVzl.jpg")
plt.imshow(I[:,:,2])
bw = I[50:250,100:300,1] < 100  # Import and binary image. Cut the region of interest (e.g. remove the white border)

plt.subplot(321)
plt.imshow(bw, interpolation='none')

plt.subplot(322)
plt.title('Distance transform of input image')
R, D = m.medial_axis(bw, return_distance=True)  # Calculate the distance transform (Matlab: bwdist)
plt.imshow(D, interpolation='none')

plt.subplot(323)
plt.title('Convex hull of input image')
C = m.convex_hull_image(bw == 0)  # Calculate the convex hull of the pixels (i.e. the are where the maxima can be)
plt.imshow(C)

masked = C * D  # zero-out all distances that are not within the convex hull
plt.subplot(324)
plt.imshow(masked)
plt.title('Masked distances according to convex hull')

M = measure.moments(C.astype(np.uint8), order=2);
cy, cx = M[0, 1] / M[0, 0], M[1, 0] / M[0, 0]   # Calculate the convex hull center of mass
x, y = np.meshgrid(np.arange(C.shape[0])-cx, np.arange(C.shape[1])-cy);
weights = np.exp(-0.003*(x*x+y*y))  # calculate a weighting mask

plt.subplot(325)
plt.imshow(weights)
plt.title('Distance weights')

masked_weighted = masked * weights

plt.subplot(326)
max_pos = np.unravel_index(masked_weighted.argmax(), masked.shape) # detect maximum pixel and its 2d-coordinates
plt.imshow(masked_weighted)
plt.plot(max_pos[1], max_pos[0], 'kx', mew=2)
plt.title('Weigthed distances with detected center')

plt.subplot(321)
plt.plot(max_pos[1], max_pos[0], 'kx', mew=2)
plt.title('input image with detected center')

plt.tight_layout()

The idea is based on an operation called "Distance Transform" and that the center of a circle is the point that has the largest distance to its boundary. The program assumes that the center of the circle is roughly in the middle of area spanned by the edges. Furthermore, it assumes that the maximum distance pixel corresponds to the circle center (i.e. other edges are closer together).  So, in rough steps:

calculate distance transform of Image
calculate convex hull of the edge pixels
caluclate center of the convex hull to detect the rough center
weight the obtained distances from 1) by the distance from the estimated center
use the strongest pixel as the center of the circle.

The code sohuld be easy to translate into Matlab, esp. as I have plotted all intermediate images.


For finite length or time limited signals, the Nyquist rate only makes sense given some finite noise floor requirement and finite total signal energy.  You look for the frequency at or above which nothing is above that noise floor and the total remaining energy is below that noise floor.  Twice that is the minimum sampling frequency to avoid aliasing causing any noise above that given floor.

For infinite length signals, there is the possibility of a frequency domain window of finite support width.


how does this equation correspond to smoothing?

The first equation you give is the difference equation for a lowpass FIR filter, or a linear filter with an impulse response that is finite in duration. I'll write it a bit differently (so that it is expressly discrete in time and causal):

$$
f_s[n] = 0.1 f[n-2] + 0.8 f[n-1] + 0.1 f[n]
$$

$f_s[n]$ is the smoothed version of the discrete-time input sequence $f[n]$, generated by passing $f[n]$ through an FIR filter with the coefficients $[0.1, 0.8, 0.1]$. The frequency response of this filter is as follows:

As it turns out, it's not a very good lowpass filter. As the name implies, a lowpass filter should pass low-frequency content while removing higher frequencies. This provides the "smoothing" action that you're looking for, as "jagged", non-smooth features are associated with high frequencies since they change rapidly with time.

Your second equation is an example of a lowpass IIR filter, a linear filter whose impulse response is infinite in duration. The filter's difference equation is:

$$
y[n] = \alpha y[n-1] + (1-\alpha) x[n]
$$

where $x[n]$ is the filter input and $y[n]$ is the filter output. This type of filter is often used as a low-complexity lowpass filter and is often called a leaky integrator. It is favored because of its simple implementation, low computational complexity, and its tunability: its cutoff frequency depends upon the value of $\alpha$. $\alpha$ can take on values on the interval $[0,1)$. $\alpha = 0$ yields no filtering at all (the output is equal to the input); as $\alpha$ increases, the cutoff frequency of the filter decreases. You can think of $\alpha = 1$ as a boundary case where the cutoff frequency is infinitely low (the filter output is zero for all time).

As an example, if $\alpha = 0.8$, the frequency response of the filter is as follows:

which is a better filter than your FIR example; it yields much better attenuation of frequencies toward the upper end of the band. Even though it might not be obvious by looking at the difference equation (because of the feedback from the filter output back to its input), it effectively performs smoothing on the input due to its lowpass nature. I'm not sure if this description will be particularly meaningful to you for your application, but these are pretty fundamental signal processing concepts; some study of introductory DSP texts could help fill in the gaps.

Edit: By request, here's a plot that shows both responses on the same axes, illustrating the relatively poor attenuation provided by the FIR example filter:


Interleaving over OFDM symbols in same transmission

Interleaving data that can be contained in one OFDM symbol between $n$ OFDM symbols will significantly reduce your bitrate and spectral efficiency by a factor of $\frac{1}{n}$. The advantage of doing this would be to increase your resilience to short bursts of noise, but only over $n$ symbols, at the cost of reducing your bitrate by $\frac{1}{n}$. I don't think many applications (if any) would benefit from this.

Edit: I think my initial response wasn't right. Suppose you're interleaving $n=10$ symbols with each other over 10 symbols, and only after you receive the whole block of 10 symbols can you read the code. The amount of information transmitted over the channel over time is the same as when symbols aren't interleaved, there is just some additional latency for the receiver to get the message. As long as your data can wait the extra time, it's worth it for your data to be that much more resilient to noise.


I suspect that you suspect that you need to window because there is a huge literature on windowing, but all things FFT, aren’t about window.

At a minimum, we use a boxcar window so it makes sense to consider other windows but if the noise in your application is flat over the band the signal occupies, the boxcar is very likely the appropriate window.

If the noise isn’t flat, a whitening filter is  typically a better idea.

Spectral leakage is a problem when you have strong interference tones and the leakage masks a weaker set of tones that are your signal of interest.  A signal of interest leaking on itself is not usually a problem.  If you do have strong interference, an adaptive algorithm will provide an effective optimal window.

If accuracy is your main concern you really should be asking yourself what the bounds on accuracy will be and that will be governed more by FFT size, how many averages and what kind of averaging you do, which also leads to issues like how much your received signal changes.

The subject can be complicated.  There are literally thousands of papers in the IEEE library that touch on the subject.  You can dedicate yourself to thousands of hours of study, or you can just try a few windows on some data and pick one.

In many situations, a customer will not be interested in details, but some customers will and you can find yourself being asked questions about issues that you didn’t know existed.

I suggest you should know what a Cramer-Rao bound is and what it is for your problem, and Pearson Correlation.

So some answers.

If you minimally averaging, a window smooths the cross correlation. If you are minimally averaging, don’t throw away data.  There are too many factors to pick a specific type. You need to experiment.

If you window, do so prior to zero padding.


Why does the representation of the filter matter? How explain it to people outside signal processing theory?

The short answer is that each for each structure there is a different sensitivity to small perturbations of the coefficients. The locus of poles and zeros change under those perturbations.

The different forms correspond to different factorizations of the polynomials that correspond to the filters transfer function.

So with infinite precision arithmetic the different forms are exactly equivalent but the factorizations are not equivalent with respect to small perturbations.

Similar ideas are found in numerical linear algebra. A Cholesky factorization is more accurate than the unfactored matrix the factorization corresponds to.  This is why rank one updates are favored in various algorithms


Rearrange the signal into a matrix containing one period in each row and average down the columns.


How do they produce an impedance/phase plot for a speaker with a chirp that lasts 1 sec?

Quoting Wikipedia article on Electrical impedance:

  Impedance is defined as the frequency domain ratio of the voltage to the current.

You can measure instantaneous current indirectly by connecting a small resistor (say 0.1 ohm) in series with the speaker and by measuring the instantaneous voltage across it, as function of time. Because the resistor is purely resistive, by Ohm's law $I = V/R$. To know the instantaneous voltage seen by the speaker, subtract the instantaneous voltage drop across the resistor from the instantaneous amplifier output voltage, which you can also measure as function of time.

Figure 1. Impedance curve test circuit. Voltage is to be measured at the test points TP2 and TP3 referenced to TP1.

In order to do the frequency domain division, play back the test signal, measure both instantaneous voltage and instantaneous current across the speaker as function of time, digitize those two signals, do a discrete Fourier transform on both, and complex divide bin by bin in the frequency domain. As long as there are no bins that are (too close to) zero-valued, the result should be meaningful. A short chirp, even a click, would work as long as it contains all frequencies of interest. The discrete Fourier transform should be long enough to fit the portion of the current signal where "things are happening". That can be significantly longer than the chirp, if it is short.

This kind of testing relies on the assumption that the system is linear. Some causes of nonlinearities (mainly at the speaker if the amplifier has good specs) can be avoided by limiting the peak amplitude of the test signal. A chirp is better test signal than an impulse, because it carries more power for the same peak amplitude, improving the signal-to-noise ratio.