Wiener Filter for Image Noise Reduction (Image Denoising)

There is another Wikipedia entry on Wiener filtering more applicable to image processing.

To summarize (and convert to 2D), given a system:
$$
y(n,m) = h(n,m) * x(n,m)  + v(n,m)
$$
where

$*$ denotes convolution,
$x$ is the (unknown) true image,
$h$ is the impulse response of a linear, time-invariant filter,
$v$ is additive unknown noise independent of $x$, and
$y$ is the observed image.

We want to find a deconvolution filter $g$ so that we can estimate $x$ as follows:
$$
\hat{x}(n,m) = g(n,m) * y(n,m)
$$
where $\hat{x}$ is an estimate of $x$ that minimizes the mean square error.

In the frequency domain, the transfer function of $g$, $G$ is:
$$
G(\omega_1, \omega_2) = \frac{H^*(\omega_1, \omega_2)S(\omega_1, \omega_2)}{|H(\omega_1, \omega_2)|^2S(\omega_1, \omega_2) + N(\omega_1, \omega_2)}
$$
where

$G$ is the Fourier transform of $g$,
$H$ is the Fourier transform of $h$,
$S$ is the mean power spectral density of the input $x$, and
$N$ is the mean power spectral density of the noise $v$.

The equation for $G$ can be re-written as:
$$
G(\omega_1, \omega_2) = \frac{1}{H(\omega_1, \omega_2)}
\left[ \frac{|H(\omega_1, \omega_2)|^2}{H(\omega_1, \omega_2)|^2 + \frac{N(\omega_1, \omega_2)}{S(\omega_1, \omega_2)}}
\right]
$$
So the Wiener filter has the inverse filter for $H$, but also a frequency-dependent term that attenuates the gain based on the signal to noise ratio.


If you have a signal of 9 samples, namely $ {\left\{ {x}_{n} \right\}}_{n = 0}^{8} $ its DFT is given by $ {\left\{ {X}_{k} \right\}}_{k = 0}^{8} $

Now if the sampling rate is $ {F}_{s} $ the bin frequency interval is given by $ \frac{{F}_{s}}{N} $.

Now, you have 9 samples, 1 is the DC ($ k = 0 $) and the rest spans $ \left[ \frac{-{F}_{s}}{2}, \frac{{F}_{s}}{2} \right] $.

You can look on the first 5 which are $ \left\{ 0, \frac{{F}_{s}}{9}, \frac{2 {F}_{s}}{9}, \frac{3 {F}_{s}}{9}, \frac{4 {F}_{s}}{9} \right\} $.

Now, the 7 bin would have to be symmetric and hence would be $ \frac{-3 {F}_{s}}{9} $.


Identifying impulse based on experimental data

Compute the position function from the accelerometer data using a scaled aggregator of some kind presumably, with initial conditions about resting position and velocity of the sound bowl. If you don't want to model the hammer's impact exactly, and want to think of it as an impulse, then the position function is the impulse response. Otherwise, go to 2.
Model the impact of a hammer as a non-ideal impulse function in time. Deconvolve the position function and the hammer model to get the impulse response. You could try to use the 3d fft on both these functions and find their ratio, and the compute its inverse 3d fft.


The answer to your suplemental question: watch Openheim's free youtube digital signal processing MIT lectures, especially the lectures number 2,3,4, and you will get answers on the first question as well. The lectures are well grounded in math, and I am sure that this is the most efficient way to learn DSP on your own, from the point of view of a mathematician.


I just need to learn the basic nuts and bolts of signals without having to go in depth into the math of deriving the Fourier Transform or signal design etc. Is that a possibility, and what is a good resource for this?

Sorry, I don't think it works this way. Without a proper understanding of the fundamentals, it's unlikely what you will every get useful results of your analysis. Check https://ccrma.stanford.edu/~jos/filters/ or http://www.dspguide.com/ for starters. This is good too, if you are up for it https://ocw.mit.edu/resources/res-6-007-signals-and-systems-spring-2011/introduction/


Understanding sine wave generation in Python with linspace

Change the following line :

  time1 = np.linspace(0, capture_size1 * timestep1, capture_size1)

To the following:

  time1 = np.linspace(0, capture_size1 * timestep1, capture_size1, endpoint=false)

You will see correct results. Your original time instances is not what you intend because Python will create 2048 equally spaced point between 0 and 2048*Ts. What you want is equally spaced 2048 points starting from 0 and 160usec apart.

You can also use the following line if you dont want to use 'endpoint=false' :

  time1 = np.linspace(0, (capture_size1-1) * timestep1, capture_size1)

At your current sampling period of $160\mu sec$, it will mean that you are taking 1 sample of the sinusoidal every 2048 periods. There will be aliasing, but you are not bothered about that because you want to see fixed amplitude across discrete time. Basically, you want aliasing to happen.


Direction of arrival and distance to source

Your statements are a little confusing. First you say that the sensors and the source all are on the same XY plane and then you mention that the source signal $s(t)$ is transmitted from a point source in far field. I will go with the latter.

In array signal processing, the source's or the sensor's position is usually estimated using the concept of relative time delay. For a given source at location $\boldsymbol{\beta} = \left[\phi, \theta \right]$ where $\phi$ is the azimuth angle and $\theta$ is the elevation angle, the signal received by the $k^{\text{th}}$ sensor at position $\boldsymbol{p}^{\text{T}} = \left[x_k,y_k,z_k=0 \right]$ will experience a time delay of

$$ \tau_k(\boldsymbol{\beta}) = \frac{1}{c}\boldsymbol{p}^{\text{T}} \begin{bmatrix}{\cos\phi \sin\theta \\ \sin\phi \sin\theta \\ \cos\theta} \end{bmatrix} $$

To simplify this setting, let us only work with the elevation angles. Let the sensor $a$ is at location $\left[x_a, y_a \right]$ and the sensor $b$ is at location $\left[x_b, y_b \right]$. So the distance between them will be $d = \sqrt{\left(x_a-x_b\right)^2 + \left(y_a-y_b\right)^2}$. Let the plane wave reaches sensor $a$ first. Solving a basic geometry will tell you that the extra distance the plane wave has to travel to reach $b$ is $\delta = d\tan\theta$. So the time delay will be

$$\tau = \frac{\delta}{c} = \frac{d\tan\theta}{c} $$

So to answer your first question, Yes! It would be possible to estimate the direction of arrival at sensor $a$ and/or $b$. You do not need to get it for both the sensors as they are receiving the signal from the same source. You would need to know the time delay which you can get by cross-correlating the signals received by the two sensors. Once you know the time delay, you can use it to get the direction of source i.e., $\theta$.
To answer your next question, Yes! Once you find out the position of the source. Find out the time delay between the sensors $a$ and $c$ and the sensors $b$ and $c$ and then apply the same equation. You should pay attention to the sign of the delay i.e., between sensors $a$ and $c$ whether is reaches $a$ first or $c$ first and same with $b$ and $c$. Since the formulation only has the distance $d$ and not the specific coordinates, the absolute value of the delay will be the same in the circle of radius $d$ with origin being the sensor $c$.
I have included the formulation of the delay in terms of the 3-D coordinates. There are a lot of good resources available online to help you understand this in detail. A basic google search will go a long way.

There are a lot of other ways to solve the direction of arrival estimation problem. You should find out more and select the one based on your needs.

Hope this helps!


Clear examples of Fixed Point Algorithm Specification Block Diagram for FPGA Implementer

I designed a lot of DSP algorithms and I would usually be the one implementing them in an FPGA. I rarely had to explain the internals of an algorithm to other people. That beind said, it is a good idea to represent the number of bits (I usually used the s:m:f notation from Sony Playstation 2) in a block diagram/design document for every components. You should also specify whether or not the FPGA designer must truncate or round.
Otherwise if you leave the number of bits used for the internal calculations up to the FPGA designer, you're looking for trouble. Using an insufficient number of bits can render your algorithm useless or unstable.
It is also a good idea to create a fixed-point model of your algorithm in Matlab/Python/Octave/C++ for 2 reasons. First, you will be able to validate the number of bits needed so that your fixed-point algorithm behaves similarly to the floating-point arithmetic algorithm. Secondly, you will be able to compare the FPGA implementation to your fixed-point model. The FPGA implementation should behave exactly the same as your fixed-point model.


How do you normalize a sample-by-sample processed signal?

Calculate the right audible half volumes in dB and add those?

Can't do that: dB is a logarithmic thing. Adding logarithms means multiplying the values. You want to build a linear filter, not a self-intermodulation generator, don't you?
If your filter coefficients are really just [1 1], well, then multiplying by 0.5 exactly retains your original signal power for low-frequency signals. Since this is a low-pass filter, of course higher frequencies than DC are attenuated, which is not a fault of the 0.5, but exactly the result of your filtering.
A usual way to avoid clipping that only happens rarely  would be a soft-clip function that lets through any "benign" amplitudes and starts saturating "softly" towards the maximum amplitudes.
If you want to do something more "time-smooth", that would for example be an adaptive gain control (AGC), which reduces gain gradually as your signal works towards higher peaks. However, that's not a "sample-by-sample" thing, but something with at the very least memory (and often even implemented in a somewhat acausal manner).


Finding shapes / lines is a good idea. Threshold tweaking is almost certainly a bad idea.

Hough transform typically expects your lines to be (1) in a binary image and (2) of one pixel width. Try to apply an edge detector (Matlab: edge) on your filtered image in order to get the boundaries of your white regions as one-pixel lines. Then Hough transform should find something.

I would also consider doing edge detection on your original image rather than your feature image. Then you can try to tweak the edge detector, rather than first running your "certain operations" and then detecting edges.

Once you have some lines you will face the problem of finding the good ones. You could look at (i.e. google for) "vanishing point detection", and read about RANSAC.


Low samples numbers machine vision training for object detection

If the camera has a shallow depth of focus and the background is mostly blurry and the foreground is only moving at one focal plane, why not use a simpler method, such as Template Matching?

Since the printing head is changing its angle of view with respect to the lens as it moves in the $x,y$ dimensions, you might have to use more than one templates but overall the process is straightforward: Apply the cross correlation in a sliding window across the image and stop where the cross correlation peaks.

This example is in MATLAB but easily transferable to other platforms.

This sliding window might seem slow but is unavoidable. Neural networks do it implicitly at the lowest level of inputs. The good news here is that the printing head image does not move in the $z$ dimension, in which case the scale would change too. This would mean that you would have to repeat this process for the scales of interest.

You might be able to derive "hints" of where the printing head is roughly in the image, by tracking its silhouette. If you threshold the image for example, you might want to look for a distinct pattern of large rectangular areas in black, that roughly match the size and shape of the arm. This would help in making the search area where the normalised cross correlation is applied, smaller. But overall, I would not expect a modern Desktop or even SBC to struggle with an HD image and a patch even as big as 128x128 pixels.

Hope this helps

EDIT:

In terms of a pesudo-code instruction:

Take a representative image of the nozzle (don't crop it too much around the nozzle but make sure the bit you want to recognise is clearly visible).
Given an image from the monitoring camera, run a normalised cross correlation across the image. How you are going to do this depends on the platform you are using but MATLAB, Octave, Python all have means of running a normalised cross correlation, fast, in this way
The output of this will be a Real two dimensional array that will have a single peak exacty at the point where the template had the  best match.
Threshold that match and take the centroid of the region (if it spans more than one pixel) to find out where the nozzle is in the camera image.


Isolate g(t) starting from $\int_{-\infty}^{+\infty}g^{2}(t)e^{-i2\pi f t}dt$

The integral is the Fourier transform of $g^2(t)$, so an inverse Fourier transform can recover $g^2(t)$. But without any further information you can't recover $g(t)$ from $g^2(t)$. So the problem is not the integral transform but simply the squaring operation.


In an attempt to solve the question on why to normalize, and implicitly how to normalize:

$F$ is your reference patch, $I$ is the patch under inspection.

So make $I$ only consist of the maximum possible value in your image format. Then, your $SSD$ becomes very large, although $I$ isn't "similar" to $F$ at all.

Now, for comparison, set $I=F$. In theory, this should give a large (if making any sense at all, the largest) $SSD$ possible, right?

So, obviously, your $SSD$ formula as is isn't useful at all for comparing things, because an $I$ with a high average amplitude will just come out as a "winner", always. So, you need to get the energy of the $I$ out of your formula first. That's easy: Just divide $I$ by the sum of the squares of all its values.


Estimating device height from 3-axis accelerometer data (possibly by supervised learning)

Any tiny offset error will cause a double integration to fly off quadratically in time.  Even if you manage to calibrate out this error for one sample test setup, any drift in offset due to change in temperature, orientation, age, power fluctuations, etc., will eventually render that calibration useless.

Continuous recalibration (airspeed, compass, optical fix, etc.) is likely required for a more useful result over a sufficiently long time period.


Clipping is a non-linear operation. Consider a system that clips all input whose absolute value is larger than 1: $$ y(t) = \begin{cases}x(t)\text{ if }|x(t)|\leq 1 \\ 1 \text{ if } x(t) > 1 \\ -1 \text{ if } x(t) < -1 \end{cases}$$ While this system is linear for small inputs, it becomes non-linear as soon as the input becomes large enough.

Non-linear systems with no memory can be modeled by a power series: \begin{align}y(t) &= \sum_{n=0}^\infty a_n x^n \\ &= a_0 + a_1x(t) + a_2x^2(t) + a_3 x^3(t) + \cdots \end{align} Usually the coefficients $a_n$ become small as $n$ increases. Let's focus on $n=2$. If the input is a sine wave, then the clipped output contains the square of the sine wave, which has a DC term plus a term of twice the frequency of the input. It is the first distortion-causing harmonic.

If you play a bit with the equations you'll see that the terms for $n=3,4,\ldots$ each introduces new harmonics. When the input is more complex, such as a sum of sine waves of different frequencies, the harmonics will be more interesting.

Another, maybe more intuitive way to see it is as follows. A clipped sine wave looks a bit like a square wave -- the more clipping there is, the more similar it is. It is natural, then, that its Fourier series will resemble that of a square wave, which has plenty of harmonics.

If you have access to a spectrum analyzer, you can easily experiment with this by producing clipped sound in your computer, and seeing the output in the spectrum analyzer.


generating spectrograms in python with less noise

None of these spectrograms show music. I would validate your algorithms by substituting a known signal, e.g. a 440 Hz pure tone. That should be a single line in your spectrogram. Audacity can generate that as a .wav file, but you can also numerically generate the signal in Python. Next, mix two known tones and check if your spectrogram is lineair (i.e. showing two lines).

Audacity can also show you what the spectrogram should be for your inputs.


FMCW Radar Signal Processing Flowchart / Ambiguity Function

The Ambiguity Function is just a name for 2D Correlator.

It is known that given a shifted 1D Signal the optimal estimator for the shift (Range in RADAR / LIDAR, etc...) it the correlation with the signal.

What if our signal had 2 parameters and both are shifted?
In the case of RADAR we have a shift of the signal due to the range (Time Shift) and due to velocity (Frequency Shift - Doppler Effect).

So how would you solve this?
Well, correlation with the transmitted signal shifted on both parameters.
The result of this operation is the Ambiguity Function which is a grid (Discrete) saying the correlation value between the received signal and a shifted (In both parameters) transmitted signal.


How to recover complex values (downconvert) from FM modulation (composite signal)?

I am not quite sure if you are asking how to perform the demodulation or how to down-mix a signal at some intermediate frequency to baseband prior to demodulation.  I will assume the prior.

If you have the modulated signal it is quite easy to perform demodulation and recover the original signal.

Simply take the angle (atan2) of the output of a delay discriminator.  If you are concerned with efficiency you can get fancier with your implementation.  Example python code would look like the following:

# Calculate the complex vector between two adjacent data points
vectorDiff = data[1::1] * np.conjugate(data[0:-1:1]);
# Record the angle of the complex difference vectors
fmDemod = np.angle(tmp);

I recently wrote a blog post detailing an FM demodulator I wrote using Python.  Here is a link.  I used this project as an excuse to learn python and numpy.

** EDIT ** Quick and dirty way to perform down-conversion using the hilbert transform prior to FM demodulation (in python)

# Import the plotting library
from matplotlib import pyplot as plot
import scipy.signal
import numpy as np

# Define the setup
fCarrier = 10;
fAudio = 1;
fs = 1000;
timeEnd = 1;
time = np.linspace(0,2,fs*timeEnd);

# Create the signals
carrier = np.sin(2*np.pi*fCarrier*time);
audio = np.sin(2*np.pi*fAudio*time);
audioInt = -np.cos(2*np.pi*fAudio*time);
freqMod = np.sin(2*np.pi*fCarrier*time + 2*np.pi*1*audioInt);

# Downconvert
analyticSignal = scipy.signal.hilbert(freqMod); # wikipedia analytic signal
baseband = analyticSignal * np.exp(-2*np.pi*fCarrier*time*1j); # complex mixing
audioDemod = np.angle( baseband[1::1] * np.conjugate(baseband[0:-1:1]) ); # fm demod


how to take a second order transfer function without any zeros and represent it using two first order transfer functions

For $|\zeta| \le 1$, let $\zeta= \cos\theta$, so $\theta=\mathrm{arccos}\,\zeta$
$$\begin{align*}g(s) &= K\frac{\omega_n^2}{s^2 + 2\zeta \omega_n s + \omega_n^2}\\
\\
&= K\frac{\omega_n^2}{s^2 + 2\omega_n s \cos\theta+ \omega_n^2(\cos^2\theta +\sin^2\theta)}\\
\\
&= K\frac{\omega_n^2}{(s + \omega_n \cos\theta)^2+ \omega_n^2\sin^2\theta}\\
\\
&= \sqrt{K}\frac{\omega_n}{s + \omega_n \cos\theta + j\omega_n\sin\theta} \cdot \sqrt{K}\frac{\omega_n}{s + \omega_n \cos\theta  -j\omega_n\sin\theta} \\
\\
&= \sqrt{K}\frac{\omega_n}{s + \omega_n e^{j\theta}} \cdot \sqrt{K}\frac{\omega_n}{s + \omega_n e^{-j\theta}} \\
\\
\end{align*}$$
For $|\zeta| \ge 1$, since $\cosh\theta \ge 1$, let $\zeta= \pm\cosh\theta$, so $\theta=\mathrm{acosh}\,\left(\pm\zeta\right)$
$$\begin{align*}g(s) &= K\frac{\omega_n^2}{s^2 + 2\zeta \omega_n s + \omega_n^2}\\
\\
&= K\frac{\omega_n^2}{s^2 \pm 2\omega_n s \cosh\theta+ \omega_n^2(\cosh^2\theta -\sinh^2\theta)}\\
\\
&= K\frac{\omega_n^2}{(s \pm \omega_n \cosh\theta)^2- \omega_n^2\sinh^2\theta}\\
\\
&= \sqrt{K}\frac{\omega_n}{s \pm \omega_n \cosh\theta + \omega_n\sinh\theta} \cdot \sqrt{K}\frac{\omega_n}{s \pm \omega_n \cosh\theta  -\omega_n\sinh\theta} \\
\\
&= \sqrt{K}\frac{\omega_n}{s \pm \omega_n e^{\theta}} \cdot \sqrt{K}\frac{\omega_n}{s \pm \omega_n e^{-\theta}} \\
\\
\end{align*}$$


Is DSP a branch of Computational Science?

Nope. Its not . Although , computational sciences covers very wide range of subjects but to be able to implement DSP , one must have a thorough understanding of signal fundamentals . That is not a part of computational sciences .


Difference equation system identification

Note that for a small sampling interval $T$, $\big(d[k+1]-d[k]\big)/T$ is a good approximation for the velocity. So if you fit $au[k]+b$ to a given set of measurements $v[k]$, it is valid to conclude
$$d[k+1]=d[k]+T\big(au[k]+b\big)\tag{1}$$
In the text you refer to they might have normalized $T$, so it changes the units without changing the values of $a$ and $b$.


Still not making any sense to me.  This is what I have parsed from what you are saying:
Red Line = Noise + Gaussian
Blue Line = Red Line + Lower Gaussian
Blue Line = Noise + Gaussian + Lower Gaussian
Blue Line = Noise
Red Line = Original
Blue Line = What you got
Clear as mud.  What you got isn't what is original.  You arbitrarily chose a spot to add a contrived Gaussian to the red line to get the blue line, yet that is what you started with.
You call the red line a gaussian + noise, but it doesn't include the gaussian you added.
Let's start with this.  Show the data set you are starting with.  The original data, the time series as you will measure it.
I'll take another look in the morning, and probably delete this answer, but it was way too long for a comment.
In the meanwhile, maybe a look at this:  Estimate peak width from a vector that is a superposition of unknown number of identical Gaussian peaks with different heights?


How do I compare 2 audio files, returning a percentage of difference

Could you please define the "percentage of difference"?

For example, please tell us which value do you expect to get:

Between two versions of the same file encoded with different settings?
Between two versions of the same file, except that one is slightly faster than the other?
Between two recordings of the same audio input, but through different microphone/recording equipment?
Between two recordings of the same person saying the same word?
Between two recordings of the same person saying the same word, with a different prosody (rhythm/melody)?
Between two recordings of different persons saying the same word?
Between two recordings of the same person saying the same sentence except for a few words?
Between a recording of a person and that of a cow?
Between a recording of an airplane and music?

What do you want to measure? The speech contents (words)? The melody, rhythm? The overall audio similarity? Encoding/recording equipment differences?


Addition of signal variance and noise variance

If the signal and noise are uncorrelated then the variance of the sum of the two equals the sum of their variances. So, yes, you can add variances (of uncorrelated signals). However, you cannot add dB-values! First of all, a variance in dB does not make sense because dB quantifies the ratio between two values. What you probably mean is dBm or dBW. So before adding variances you have to compute their actual (linear) values from the given values in dB(m/W).


The following depicts a QPSK Carrier Recovery loop that works on the baseband I and Q samples following proper decisions from a Timing Recovery Loop. For more information on this implementation see

High modulation index PSK - carrier recovery

And for more information on a Timing Recovery implementation that can precede this see:

Symbol timing synchronization using a high sampling rate

(Copying the links here since the question asked would not lead you easily to the prior titles)


FFT of two sets of samples vs FFT of sum of samples

Fourier Transform is linear, hence if $\mathcal F[x(t)]=X(f)$ and $a$, $b$ are complex numbers, then:

$$\mathcal{F}[ax(t)+by(t)] = a X(f) + bY(f) $$

So in your case, simply sum the results together.


I may have fond the answer. The $Q$ matrix, made from the sum of the Hankel and Toeplitz forms, stems from:

$$cos(k\omega)cos(n\omega) = \frac{1}{2}\left[cos((k-n)\omega)+cos((k+n)\omega)\right]$$

but with the additional $0.5$ for even orders, it comes out as:

$$cos((k+0.5)\omega)cos((n+0.5)\omega)=\frac{1}{2}\left[cos((k-n)\omega)+cos((k+n+1)\omega)\right]$$

and since the $k+n$ term makes the Hankel matrix, it needs an increment in its indices, and the $q$ vector longer by 1. The $b$ vector only needs the extra $0.5$. Now the relevant code lines look like this:

q = [wp + K*(1 - ws), wp*sinc(wp.*n) - K*ws*sinc(ws.*n)];
Q = (toeplitz(q(m + 1)) + hankel(q(m + 1 + oddN), q(m + M + 1 + oddN)));
b = wp*sinc(wp.*(m + 0.5*(~oddN)))';

and this is a comparison between a 30th and a 31st orders, with f=[0 0.3 0.4 1], A=[1 1 0 0], K=[1 10]:

What I don't understand now is why the $q$ vector is not made with the $0.5$ sampling offset.


A FIR filter is (similar to, or) just a weighted moving average filter, where (unlike a typical equally weighted moving average filter) the weights of each delay tap are not constrained to be identical or even of the same sign.  By changing various values in the array of weights (the impulse response, or time shifted and sampled version of the same), the frequency response of a FIR filter can be completely changed.


Reconstruction of Audio Signal from its Absolute Spectrogram

One thing commonly done (for example in the source separation community) is to use the phase data of the original signal (before transformation where applied to it) - the result is much better than null or random phase, and not so far from algorithms aiming at reconstructing the phase information from scratch.

A classic reconstruction algorithm is Griﬃn&Lim's, described in the paper "Signal estimation from modiﬁed short-time Fourier transform". This is an iterative algorithm, each iteration requires a full STFT / inverse STFT, which makes it quite costly.

This problem is indeed an active area of research, a search for STFT + reconstruction + magnitude will yield plenty of papers aiming at improving on Griffin&Lim in terms of signal quality and/or computational efficiency.


You can, if you increase phase between samples slowly enough, using unwrap(angle(signal)). "Slowly enough" means the phase doesn't jump by more than $2 \pi$; unwrap works by tracking "total phase". Python equivalent implementation here (note you can configure discont there for jumps greater than $\pi$ (or TOL) but it's pointless with angle since it outputs within $[-\pi, \pi]$).
Otherwise, can't, as described by others. For an infinite number of inputs, you get the same output, like asking "what's $N$ in $N\ \%\ 2 = 0$" (division by 2 remainder).


Why does an anti-symmetric function has zero amplitude at the center of an even length window

A real symmetric function has the following property $$f(t) = f(-t)$$

whereas a real anti-symmetric function has the following property $$f(t) = -f(-t)$$

So at $t=0$ you have $$f(0) = -f(0)$$ which implies that $$f(0) = 0$$

A complex conjugate-symmetric function has the following property $$f(t) = f^*(-t)$$

whereas a complex conjugate-antisymmetric function has the following property $$f(t) = -f^*(-t)$$

So at $t=0$ you have $$f(0) = -f^*(0)$$ which implies that $f(0)$ is real and further $$f(0) = 0$$


I assume that "small event" means "short (in time) event".  Fourier transforms are not good instruments for detecting short events.  You would be better off looking for such events in the time domain, or with a wavelet transform, perhaps, that is set up specifically to detect the events you are interested in.

That being said, if for some reason the events have to be detected in the frequency domain then I would keep the transforms short.  The reason for this is if the event is short with respect to the transform length it will look like a non-ideal impulse, which basically just creates a lot of high frequency content.  Because of the small amount of energy (due to the short time) it will probably be difficult to pick the signal out of the noise.

If the transform is short, on the other hand, then its energy is a significant part of the total transform energy.  That will make it possible to pick the signal out of the noise, or even enable detection by measuring the total energy.


Why does my GAIN remain constant after a few cycles?

$P$ is the state covariance matrix, not the Kalman gain, though the Kalman gain $K$ is usually calculated using it.

It is normal in simple Kalman filters (ones that do not vary with time) for the state covariance matrix to stabilize after a while so that it is constant.

See, for example, the plots in Figure 3(b) and 3(c) in this PDF which show how both the Kalman gain $K$ and the state covariance $P$ go to a constant value after a relatively short space of time.

Because this happens, the Kalman filter evolution is often ignored in practice and the steady state version is what is implemented.


Filtering data so that only rising edge is left

Try looking into posts on "onset detection", where typically "flux" (or first order difference, typically of spectral magnitude) is used to compute the "rate of change" of a signal. The high rate of change is the most salient feature of the images you drew above.

However, your third image seems to translate the maximum of the curve (to the left) is that your intent?


Removing halo effect in image composition: Google Meet's solution

Based on what I read in your reference I would say that
The key part for this solution is the model to compute the mask, if you have one that is a good starting point.
The second important feature is that it uses a position dependent filter, and the filter is parameterized with the length of the value of the mask. In their word avoiding the foreground to bleed in the foreground. They don't know how exactly, but I would say that from the mask a parameter $r_{i,j}$, that is large for the foreground and small for the background.
The possible realization of separable filter would be
$$ v_{i,j} = \frac{ \sum_k f(k\cdot r_{i,j}) x_{i-k,j}}{\sum_k f(k\cdot r_{i,j})} $$
$$ y_{i,j} = \frac{ \sum_k f(k\cdot r_{i,j}) v_{i,j-k}}{\sum_k f(k\cdot r_{i,j})} $$
And you are free to play with the normalization factors.
Using separable filter means that you use two 1D filters rather than small 2D filters for the image in multiple resolutions. This may be better in terms of memory. It can be implemented with two copies of the image, $x$ is written to $v$, and the second equation update the image $x_{i,j}$ in place.
Other possible advantage, since the filter is computed for every pixel separately. A separable filter of size $N \times M$, would require $N + M$ evaluations of the filter, while a 2D filter would require $M \times N$ evaluations, and in the realization above since the horizontal and vertical filters are the same it will require only $N$ evaluations.
Having the mask dependent filter is important to pass the foreground (the faces) with little distortion.


What is the Compression Algorithm used by Facebook?

From the images you posted its quite clear that the image has been downsampled and re-compressed with lower quality jpeg settings.  If you look round the mouth you can clearly see JPEG-like artefacts.


Rate: the actual rate you transmit data at.

Throughput: the percentage of signal data that is actually delivered correctly to the receiver in the overall data (signal data + coding + overhead ... etc).

Capacity: maximum rate that a channel allows for a given SNR such that the received signals can be reliably decoded at the receiver.


Reconstructing the song using its dominant frequencies

Short answer: You can't.

Sound is more than its dominant frequencies, and you threw away anything but that.


The figures are correct. You can see that there are samples at indices $0,1,\ldots,N-1$. There are no samples shown at index $N$, neither in the time domain nor in the frequency domain. The value $N$ is only shown on the abscissa because in the frequency domain it corresponds to the sampling frequency, and in the time domain it corresponds to the end of the continuous-time signal represented by $N$ samples (if we assume that each sample represents a portion of length $\Delta T$).
For a very detailed discussion on the definition of the actual duration of a discrete-time sequence have a look at the answers to this question.


How can a linear operator on DFT vector produce the same vector using only half of the DFT vector?

I've suggested in the comments above (and I think Matt L is trying to say the same thing also) that you should be able to just discard the negative frequency elements entirely for your analysis. They don't seem to be relevant to what you're doing (which makes sense, since they don't carry any information, as the time-domain signal is real). Look at your equation again:

$$
\mathbf{y} = \mathbf{Tx}
$$

(I changed the casing so that lower-case symbols are vectors and upper-case are matrices, as is the typical convention)

Assume for the time being that $\mathbf{x}$ and $\mathbf{y}$ are length-$N$ and thus $\mathbf{T}$ is $N$-by-$N$. This would correspond to the general case, where the time-domain signals corresponding to $\mathbf{x}$ and $\mathbf{y}$ need not be real. If $\mathbf{y}$ is a column vector, then you can think of the matrix operation instead as:

$$
\mathbf{y} = \sum_{k=1}^N x_k\mathbf{t_k}
$$

where $x_k$ is the $k$-th element of the vector $\mathbf{x}$ and $\mathbf{t_k}$ is the $k$-th column of the matrix $\mathbf{T}$. My question at this point would be, instead of using the full $N$-length vector and $N$-by-$N$ matrix, why not just use the $N/2$ rows that contain information? Then, the calculation would turn into:

$$
\mathbf{y^+} = \sum_{k=1}^{N/2} x_k\mathbf{t_k}
$$

(neglecting the details of the boundary cases for the zeroth and Nyquist bins, which depend upon whether $N$ is even or odd)

The resulting vector $\mathbf{y^+}$ contains the positive frequencies for the desired result. You can then just symmetrically extend it to get a full length-$N$ vector that has the conjugate symmetry that you seek. This should make intuitive sense in the following ways:

If there is no information in the negative frequencies in the input vector $\mathbf{x}$ (which there isn't, since the entire vector is defined by its positive frequency components), then you shouldn't need to use them in your calculations.
If there is no information in the negative frequencies in the output vector $\mathbf{y}$ (which there isn't, since the entire vector is defined by its positive frequency components), then you shouldn't need to waste time calculating what their values would be.

This should free you from worrying about any bizarre constraints on the properties of the matrix $\mathbf{T}$. I think you're trying to think of the action of the precoder as too coupled with the OFDM modulation. Instead, think of it just as a preprocessing operation that happens on the encoded symbols before you assign them to frequency bins in the OFDM modulator's IDFT.


How do I pick an optimal window size for overlap add/save method algorithmically

The optimal size is likely not just algorithmic, but depends on the cache sizes and levels and the compute-to-memory-bandwidth ratio for the particular processors and systems(s) on which you plan to run this code.  Also, your applications latency requirements (which filling big FFTs add).


How to understand the de-correlation property of DCT? What does de-correlation mean?

The term "decorrelation" generally accounts for processing that reduce autocorrelation within single, or cross-correlation within a group of  signals or  images. Doing that, it should preserve important features in the  data.

In other terms, images often contains simpler "objects" (in terms of  morphology: bumps, edges, textures) that are mixed, or distorted. Generally, unmixing them via decorrelation somehow simplify further processing.

Generally yes, unless you have uncorrelated noise in your images. A classical and simple model is an order-one Markov or autoregressive process. Basically a pixel somehow  depends (with a parameter $\rho$) on the past pixel, with uncertainty. If $\rho$ is close to $0$, you are already decorrelated. If $\rho$ is close to $1$ (typically $0.90-0.95$), the  data is quite correlated.

The DCT was meant for diagonalizing the  resulting autocovariance matrices with  Toeplitz structure, to give fast estimates of their eigenvectors. Transforming a fat autocovariance matrix into a thin, close to a diagonal, matrix is an instance of decorrelation.


Your denominator is wrong. It should be on this form

$$
(1-z_0z^{-1})(1-z_0^*z^{-1}) = 1-2R\cos(\phi)z^{-1}+R^2z^{-2}
$$
where $z_0 = Re^{j\phi}$.

If your goal is to block DC and create an anti-aliasing filter in a simple way, I would recommend not to use complex poles. Instead use the procedure described in detail here: http://www.dsprelated.com/showmessage/172787/1.php.

It even provides a C implementation.


Recovering a Differential Equation From the Transfer Function of a Cascaded System

You have discovered pole-zero cancellation, which is one of the reasons that state-space systems description was invented.  Another reason is because when you go from a differential equation to a transfer function, you lose information about the structure of a system that you can choose to preserve in a state-space description.

In general, you should be very careful when you see a pole-zero cancellation, because even if the pole is stable, it can still affect the behavior of the system on startup, if there is a transient that pushes the system into nonlinear behavior and excites the extra pole, or if the pole and the zero don't actually match up perfectly in the real world.

And yes, if you have an unstable pole that's cancelled by a zero then Bad Things will happen to the system that simply aren't described by the transfer function.


System Identification: smoothing in E-step of EM algorithm

The result of the $M$-step is expressed in terms of the matrices $\mathbf{P}_n^2$ and the vectors $\mathbf{x}_n^s$. These quantities are defined as expectations. Smoothing is used to compute these expectations, because in practice expectations are usually approximated by time averages (hence 'smoothing'). If you have another method for computing these expectations, the results of the $M$-step will of course be different, but the expressions for the results of the $M$-step remain valid, you just plug in different values for the respective matrices and vectors.


The statement $g[n]=f[-n]$ is meaningless unless you assign some numerical value to $n$, e.g. $g[3]=f[-3]$, and if you do say that $n=3$, then the equality
$g[3] = f[-3]$ tells you nothing about the relationship between $g[n+k]=g[3+k]$
to the value taken on $f[\cdot]$ for any choice of argument.

What you really meant to say is that

  $$\text{For all integers } n, g[n]=f[-n]$$

which asserts that the equality holds for all choices of $n$, and in
this case, it follows by definition that $g[n+k] = f[-n-k]$;
you just need to negate the argument $n+k$ of $g$ (to get $-n-k$)
and use it (it meaning $-n-k$) as the argument of $f$ exactly as
the definition says you should.  The
complicated arguments in another answer that claim that $g[n+k]=f[-n+k]$
are wrong, as discussed in the comments following that answer.


Drawing conclusions from impulse response of a discrete LTI system

You're right about your conclusions. And now a few hints to help you solve the questions:

Your textbook should state somewhere the definition of even and odd parts of a (real-valued) sequence: $$h_e[n]=\frac12\big(h[n]+h[-n]\big)\\h_o[n]=\frac12\big(h[n]-h[-n]\big)$$

In this case (FIR filter), the difference equation is simply given by the convolution of the input sequence $x[n]$ with the impulse response $h[n]$: $$y[n]=\sum_kx[k]h[n-k]$$


Image Filtering by Averaging Similar Areas of the Same Image - Patch Based Processing

Those kind of algorithms are called Non Local algorithms.
The most known algorithms of this family is the - Non Local Means which is a decent Noise Reduction (Denoising) algorithm.

Until the Deep Learning boom, this approach has been extended and usually means working in the Patch Space of image - Patch Based Models and Algorithms for Image Denoising: A Comparative Review Between Patch Based Images Denoising Methods for Additive Noise Reduction.

All those are based on a model in patch based - Average by some kind of a metric / Similarity or some kind of Low Rank model


wavelet packet transform and lifting scheme?

A lifting scheme is a method for splitting a sequence of discrete samples into downsampled subsequences, so that you can predict a subsequence from the other, and update the former later, using possibly nonlinear predict and update operators. This corresponds to a generalization of one level of discrete wavelet transform, usually obtained fully linearly by convolution (first-generation wavelet).

Lifting wavelets were called second-generation wavelets, and proved to be a generalization of first-generation wavelets.

This being said, for question 1:

yes, a wavelet packet being iterations of several one level DWT on other subbands.

For question 2:

A lifting discrete wavelet packet can be  critically sampled, like the DWT, but also can be undecimated or oversampled if needed. And they are interesting questions about what kind of  stationary lifted wavelets can be built.

For question 3:

Generally, the Heisenberg boxes are considered to be of constant area. With the STFT, titles have the same shape, whereas with wavelet they are dilated/stretched in each direction. With a uniform wavelet packet decomposition (where the time-scale plot is regularly discretized), the time-scale bins can look similar to that of  the STFT. If the packet transform is critical, the result is often poorer with wavelets packets (because of the lack of redundnacy).

At the same redundancy, aliasing and constraints on wavelet shapes often yield less degrees of freedom that full-fledged oversampled filter banks, but the filter are generally easier to optimize. Nevertheless, if you want uniform frequency bin, wavelets are rarely very beneficial.


A special case of 2 jointly Weak-Sense Stationary (WSS) stochastic processes

Yes, when $R_{X,Y}(t_1, t_2) = E[X(t_1)Y(t_2)] = 0~ \forall ~t_1, t_2$, $\{X(t)\}$ and $\{Y(t)\}$
are often called uncorrelated WSS processes,
and they are considered to be jointly WSS processes.  More nitpicky folks
agree that the processes are indeed jointly WSS when $R_{X,Y}(t_1, t_2) = 0~ \forall ~t_1, t_2$, but insist on reserving the adjective uncorrelated
for processes for which the cross-covariance function
$$C_{X,Y}(t_1,t_2) = E\left[\left(X(t_1)-E[X(t_1)]\right)\left(Y(t_2)-E[Y(t_2]\right)\right] =E[X(t_1)Y(t_2)]-E[X(t_1)]E[Y(t_2)] = R_{X,Y}(t_1,t_2)-E[X(t_1)]E[Y(t_2)]$$
equals $0$ for all $t_1, t_2$, that is, every pair of random variables $X(t_1)$, $Y(t_2)$ from the processes is a pair of uncorrelated random variables.  Of course, when at least one of the processes has zero mean
(as is usually the case in signal processing applications), then
$C_{X,Y}(t_1,t_2) = R_{X,Y}(t_1,t_2)$ and so the nit gets picked without
requiring any further consideration.


Linear Chirp parameters from a segment of the signal

Amplitude $A$ is trivial. Square and low-pass filter/average; then take the square root.

Without noise (or with good SNR), you could simply calculate a derivative:

\begin{align}\frac{\partial\, s\}{\partial t} = -A\sin\left(\omega_0t + \frac\alpha2t^2 +\phi_0\right)\cdot\left(\omega_o+ \alpha t\right)
\end{align}

and low-pass filter the result

$$\text{LPF}\left\{\frac{\partial\, s\}{\partial t}\right\} = \omega_o + \alpha t$$

From the resulting linear function, $\alpha$ can directly be extracted as the slope, and knowing the slope, $\omega_0$ is just the constant summand.

With $A$, $\omega_0$ and $\alpha$ known, you can simply compare at any fixed $t$ to find $\phi_0$.


Do all Harmonics of a periodic sound start from a 0 base

It's common to not start all harmonically related sinusoids with a phase of 0 at some initial time point.  Varying the phase relationships helps reduce the peak-to-average power ratio, which allows a higher average volume without clipping or otherwise exceeding some peak limit.


FIR Filter design, using different filter length than signal length

There's no requirement that the signal to be filtered has to have the same length as the filter - in fact, that's practically never the case.
So, there's no problem here; convolution doesn't need to "define some padding"; you'll find discrete convolution is just defined as finite sum over the available lengths.
The length difference between signal and filter isn't of relevance to "track back" a phenomen in the input signal: the time is always the time in the output signal, minus the filter's group delay. Assuming your filter is linear phase / time-symmetric, that's simply half the filter length.


Trying to plot frequency response of a filter Transfer Function in MATLAB, it looks wrong

How did you derive the second figure? I got the different result using the following code:
b = 1 - r1;
a = [1, -2*r1*cos(theta), r1^(2)];
figure; freqz(b, a)

)
The only problme I can find in your code is that
TFdb= 20*log10(TF);
should be
TFdb= 20*log10(abs(TF));
Did you take the wrong sign of $a_1$ when calculating second transfer function?
b = 1 - r1;
a = [1, 2*r1*cos(theta), r1^(2)];

The above code generates the second figure.


Determining Parameters of S Function for Image Global Contrast Enhancement

Okay, I am convinced that there is a typo in Eq. (4) of the cited paper.  It should be:

$$ S(x;a,b,c)=\begin{cases}
0 & x\leq a\\
\frac{(x-a)^2}{(b-a)(c-a)} & a<x\leq b\\
1-\frac{(x-c)^2}{(c-b)(c-a)} & b\leq x\leq c\\
1 & x>c
\end{cases} $$

This sigmoid function or smoothstep function is continuous everywhere and the first derivative is continuous everywhere, but the second derivative is discontinuous when $x=a$ or $x=b$ or $x=c$.

You must have $a < b < c$.

$a$ is defined where $S(x)$ begins to transition from 0 upward.  and $c$ is defined for when the transition is over and $S(x)=1$.  $b$, which is in between $a$ and $c$ is the point of inflection for when the curvature changes from a concave curvature to a convex curvature (or maybe it's the other way around).

There are better sigmoid or smoothstep functions than this that the author of the paper uses.  I dunno why he/she is so enamored with this definition.


I believe you should use the Nyquist bin every time. By having it, you can always recreate the full DFT of your real-valued signal from the first half of spectrum. Whereas if you discard this value there is no way back and you can recreate your time domain signal.

As in answer you pointed to, usually the Nyquist bin shouldn't contain any data when while using anti-aliasing filter. Nevertheless sometimes you might want to estimate the energy of your signal from it's DFT. You might therefore get inaccurate value. Question is, do you really care about such accuracy?

Just to wrap up. I always use Nyquist bin and I don't need to ask myself these questions.


Can I construct a linear system so that it results in a moving average?

The function you are trying to solve is equivalent to
$\frac{b}{3 \cdot a} + \frac{a}{3 \cdot b}=\frac{1}{3}$

It is also
$a^2-ab+b^2=0$

but this function does not have a real number solution.


Understanding Magnitude Spectrum of Images

Many people face the same issues. Your image is far from showing stationarity: quite edgy, with circular features, so the Fourier analysis might not be so enlightening, without showing the unwrapped phase, which might contain many interesting features for this class of images. However:

so there is a dot in the middle: standard images have integer values between $0$ and $255$, hence the average is not zero. The average is the center dot. Remove it, and the global contrast may be better
some striking vertical and horizontal lines: FFT generally uses periodic extension, and the differences between the bottom and top, and right and left, might obfuscate details.

In the following pictures, you get space data and the Fourier spectrum. The first line is your image, with an horizontal/vertical cross. When you window the image (second line), those disappear, and the most luminous  decreasing slope is, probably, due to the index finger, and the parallel lines from the two biggest gears. On the third row, the windowed image is subtracted from its  median filtered version.  One begins to see circles, related to locally oriented teeth.

So, to interpret complicated data, remove first order artifacts, and slow background features, to better interpret details.


Fun problem: Dividing a dataset into parts depending on mean value

Try total variation denoising. Here I found a piecewise constant toolbox: http://www.maxlittle.net/software/pwctools.zip


This may help: the result is the cascade of an FIR filter with coefficients $h(n)$ with the reverse conjugate filter with coefficients $h^*(N-n)$. If we can limit $H(z)$ to be the unit circle (for computing the frequency response), then
$$H(z^*)|_{z=e^{j\omega}} = H(z^{-1})|_{z=e^{j\omega}}$$
since
for $z = e^{j\omega}$, $z^* = e^{-j\omega}$
And therefore
$$H^*(1/z)|_{z=e^{j\omega}}= H^*(z^*)|_{z=e^{j\omega}}$$
And
$$H^*(1/z^*)|_{z=e^{j\omega}}= H^*(z)|_{z=e^{j\omega}}$$
So
$$Y(z)|_{z=e^{j\omega}} = H(z)H^*(z)|_{z=e^{j\omega}} = |H(z)|^2|_{z=e^{j\omega}}$$

Also from the link the OP provided to me in the comments below, it is clear that the problem is dealing with the two sided Z-transform and the full problem has a solution given as
$$H(z)H^*(1/z^*)\Phi_{xx}(z) = \frac{\sigma_v^2}{(1-az^{-1})(1-a^*z)}$$
And we see from the earlier solution in the link that
$$H(z) = \frac{1}{1-az^{-1}}$$
Therefore
$$H^*(z) = \frac{1}{1-a^*(1/z^*)}$$
And
$$H^*(1/z) = \frac{1}{1-a^*(z^*)}$$
And
$$H^*(1/z^*) = \frac{1}{1-a^*z}$$
And
$$H(z)H^*(1/z^*)=\frac{1}{(1-az^{-1})(1-a^*z)}$$
Note confirming my earlier result, when we restrict z to be the unit circle (for determining the magnitude response):
$$H^*(z)|_{z=e^{j\omega}} = \frac{1}{1-a^*(1/z^*)}|_{z=e^{j\omega}}= \frac{1}{1-a^*(1/e^{-j\omega})} = \frac{1}{1-a^*(e^{j\omega})} = \frac{1}{1-a^*z}|_{z=e^{j\omega}}$$
Also to avoid any confusion, earlier in the linked solutions where it says "from 7... we arrive at (2.55) as follows" it gives a result that indicates
$$\Phi_{xx}(z) = \Phi^*_{xx}(1/z^*)$$
But this is not generally true for any z transform, (we cannot conclude that $H(z) = H^*(1/z^*)$ for all z) but it is in this case specific to the autocorrelation function $\phi_{xx}(k)$ since that function is Hermitian:
$$\phi^*_{xx}(k) = \phi_{xx}(-k)$$


Photoshop's Auto Color is nothing more than Histogram Equalization.

It has some improvements over the classic Histogram Equalization by adding a threshold of not touching the % at the extreme values.

Try it and see for yourself:

Histogram Equalization.
Histogram Equalization (Image Processing) Part 1.
Histogram.

Good Luck.


Standard representation of FFT of an Image?

For real data, the result of an FFT is conjugate symmetric, so all but quadrant A is often ignored as redundant info.  Statement 2 seems to be about just quadrant A.


SNR in one-dimensional additive channel with uniform noise and binary antipodal input

Signal power $=E[X^2]$; Noise power $=E[N^2]$. I assume that you know how to calculate these expectations.


Why does it mean that the process/signal is not stationary when its variance varied with time?

From the Wiki page - A stationary random process is a stochastic process whose unconditional joint probability distribution does not change when shifted in time. Consequently, parameters such as mean and variance, if they are present, also do not change over time.

Additional References are: Link1 and Link2.

The implication from the definition is that the mean and variance of random process do not vary over time. Since by your definition of the variance - it varies over time, then your process is non-stationary.

A couple of notes:

Just because the mean and variance don't vary over time does not imply that the process is stationary. The definition of a stationary random process has more implications.
A random process where the mean and auto-covariance are invariant to time shifts is often called a "Wide Sense Stationary" process. There are no conditions on the higher order moments (skewness, kurtosis, etc). This is a less stringent version of stationarity that is often more applicable for real-world signals.


How to implement simple noise filter based on the frequency domain information?

Fast convolution filtering with an FFT/IFFT requires zero padding and using overlap add or overlap save methods to remove circular convolution artifacts.

You should also use a frequency domain response that has a shorter impulse response than a rectangle (which zeroing bins is).


How many bins to include when calculating SNR from FFT?

The OP's question requires further details to provide a definitive answer, but the following will give the considerations involved. If the noise is white and stationary then the answer is clear in that we can simply use a power spectral density and work with the SNR as an SNR/Hz quantity. If the noise is not white (meaning the average power across all bins is not the same) then the concern will ultimately be with the measurement bandwidth of your final system that is providing the estimate of the quantity you are looking for. You would include all noise that is within this bandwidth as that is what will give you the statistical error in your result (which is why we are concerned with SNR in the first place). Usually the choice of bandwidth comes down to a time /frequency decision and the coherence time of your signal, and for how long it can be assumed to be stationary. For bandwidth limited signals that can be assumed to be stationary for a long time, then we can use a relatively longer averaging time (reducing the bandwidth down  to that of the signal, and therefore the amount of total noise in our result, as long as the noise is a spectral density and not a "spur" that exists at one frequency).
Ultimately to answer your question, determine accurately what the bandwidth is of your measurement. Use that bandwidth in the analysis where you predict the total noise and you will get a useful SNR. (This makes sense, consider the result you will get from your measurement system with a noise-free signal, and then consider the result you will get with noise alone-- the ratio will be your SNR that matters as it is based on the final result-- any noise that got filtered out by your measurement system does not apply).  This is analogous to a communication system where we will use the bandwidth of the signal, since the ideal receiver will have filtered out anything outside of that bandwidth but not reject any of the valuable signal. Your ideal measurement system should do the same, which means you will have interest in the spectral occupancy of your signal (and for how long it can be assumed to be stationary, which leads to concern with phase noise since that is typically not stationary... read on below).
You are likely doing an experiment where you are trying to estimate the presence of a specific signal, in the presence of noise, and ultimately will need to determine the best observation time in which to maximize your SNR (and therefore the accuracy of your estimate). Phase noise as your first measurement gave is a great indication of this as it is typically (if not always) NOT stationary over long durations, (and therefore for very long duration the typical power spectral density measurement is not valid). A much better measurement you can make that will allow you to optimize your estimation duration is the Allan Variance (also known as a two-sample Variance) since that will indicate the duration of time in which you can use an assumption of stationarity, and will then give you the averaging time that will maximize your SNR.
I go into more details of the Allan Variance and Allan Deviation for this particular application here: What determines the accuracy of the phase result in a DFT bin?


Choosing Drive Signals for System Identification?

For linear systems, you can completely characterize the transfer function using its frequency response, so a frequency sweep would be one possible choice. However, you would need to ensure that at each test frequency, you allow time for the system's transient response to die out before measuring its steady-state amplitude/phase response.


Computing Hartley transform of zero-padded signal from transform of original signal?

Since the convolution theorem applies to the Hartley transform as well as the DFT, you may be able to use the same technique as can be used for FFT interpolation or spectrum resampling.  Use the Hartley transform of the implied rectangular window as an interpolation kernel, and interpolate to a larger N. Perhaps window the interpolation kernel for a trade-off between accuracy and computational efficiency.

Depending on the size and efficiency of the interpolation needed to meet your requirements, it may be more efficient to do an inverse and forward transform instead.


In order to remove the noise from the signal, you need to subtract the spectral content of the noise from that of the signal ( filtered_signal = unfiltered_signal - noise ). This technique is called spectral subtraction. This will only work well if the noise is rather uniform across the entire length of the signal.

The steps required to do this are

Identify the noise region(s) in your signal and do an FFT on it. A single noise region is enough provided that it's long enough for your FFT.
Do the FFT on the entire signal. This may be done in multiple segments (frames).
In the frequency domain, subtract the result of step 1 from step 2. If you're using multiple segments in step 2, then simply subtract the result from step 1 from every segment. This will effectively filter your signal.
Do an inverse fft, which will bring your signal back into the time domain.

Note that I said that this technique will only work if the noise is uniform and by uniform I mean that its spectral content is rather constant. But you can also have regions in your audio where the noise is not present at all so obviously trying to remove the noise from these noiseless regions should be avoided. If your noise removal process is manual (at least in part), it would be possible (albeit tedious) to label such regions and then have the algorithm skip them.


How can I perform object tracking using noisy images?

If the noise levels added to each pixel are independent, perhaps blurring the images first and then applying a pairwise cross-correlation could do the job.


You want something that has lots of features in the spectral domain and in the time domain. However you have to be careful what other content or meaning you would like to convey or not. Spectra conveys pitch and/or coloration and timbre. Time transients convey rhythms and/or patterns.

For example you could do band limited noise burst with different center frequencies and patterns. That would allow yo to vary localization, spectrum and rhythm without directly conveying pitch. You could also use different percussion samples with different rhythms.


How to do people detection from bird eye view?

It appears that motion is your only cue in this case, assuming that you have a time-series of these images. I would start by trying to track these objects using a Kalman filter with a constant velocity motion model. The fountain and the street sign would be stationary. If your images contain other moving objects, e. g. cars, then you can try using speed to distinguish them from people.

Here is an example of how to track multiple objects using the Computer Vision System Toolbox for MATLAB. The example detects moving foreground objects using the Gaussian mixture models. You may want to replace that part or augment it using the depth information.


Comparing two sound files in embedded-systems (C language)

Matched filtering, or cross-correlation, is a pretty cheap and easy way to determine the similarity between signals. Otherwise you could try implementing the algorithm Shazam uses.


Reference Code for Positive Basis Pursuit Denoising

I assume you're after the following optimization problem:

$$\begin{align*}
\arg \min_{x} \; & {\left\| x \right\|}_{1} \\
\text{subject to} \; & A x = b \\
& x \succeq 0
\end{align*}$$

This is pretty simple problem if we pay attention fo the fast that given $ x \succeq 0 $ then $ {\left\| x \right\|}_{1} = \boldsymbol{1}^{T} x $.

This means the above problem is equivalent of:

$$\begin{align*}
\arg \min_{x} \; & \boldsymbol{1}^{T} x \\
\text{subject to} \; & A x = b \\
& x \succeq 0
\end{align*}$$

Now, this is a Linear Programming problem which can be easily solved using MATLAB or any other LP solver.

For instance, in MATLAB it will be something like that:

numRows = 10;
numCols = 100;

mA = randn([numRows, numCols]);
vB = randn([numRows, 1]);

vE = ones([numCols, 1]); %<! Vector of ones to sum vX

vX = linprog(vE, [], [], mA, vB, zeros([numCols, 1]), inf([numCols, 1]));


the code on the bottom looks like a triple-nested FM oscillator if i understand how the AdvanceOscillator() method works.  likely the SFreeSpinningOscillator class has a method of setting the initial phase of the oscillator.  you may get wildly different tone timbres with exactly the same code and different phases for the 4 oscillators.

if you can get your hands on a DX7 alg chart

you can see that most of their algs are sum of two chains.  alg 2 has a triple nested oscillator chain where your d is the DX7 osc #3.  but i can't tell you what the initial phases of their oscillators are.  there are so many parameters that affect the timbre, sometimes in unpredictable ways.  so doing software to emulate a DX7 will involve some tweeking of knobs


So, yes, since your decimation is a power of two, going through cascaded half-bands is a clever thing to do!

First of all, what's special about a half-band filter? If you implement a filter of 1/2 bandwidth as FIR, then you can choose to design it so that the transition width's center is exactly at ¼ the sampling rate – that allows for symmetry, and that in turn gives you a FIR filter where every other coefficient is zero.

Zero coefficients don't need to be summed. This gives you the steepness of a $K$-tap filter at the cost of $\frac K2$!

Now, let's look at the fact that you're decimating by a factor of 2:

That means that you throw away every other output sample. Polyphase decomposition of your filter allows you to write your single filter as a sum of two alternatively "fed" filters, so that instead of calculating $K$ products and sums per time step, you calculate $\frac K2$. That's another effort reduction by a factor of 2.

So, what you do is:

you pick an existing decimating FIR filter implementation in C or assembly, or write it yourself (it's but a simple for-loop). GNU Radio has one; it's in C++, and uses optimized filter kernels, but it might help illustrate how you split up one filter into $K$ subfilters to reduce the load by a factor of $K$ in a $K$-decimating FIR.
You design a half-band filter. You probably don't need something fantastic; make sure you get as many zeros in there as sensible.
you concatenate $\log_2 64=6$ of these (in fact, that's only nearly efficient; read that article, it's good).

General comment: Using an 8-bit microcontroller to process 16 bit numbers is a bad idea through and through. Atmegas are usually more expensive than similarly peripheral-equipped ARM counterparts, and really, they aren't even remotely adapted to DSP loads. Don't do that to yourself. This task is easy even on a cheap 20 MHz ARM Cortex-M0, and you could have a Cortex-M4F for the price of a larger ATMega. That would have a floating point unit; not that you need that for any of this (you can do all these calculations in fixed point, like you'd have to do it on your AVR), but it does make one's life easier further down the processing chain. I assume you use an AVR "because you already have that, and have experience with that"; but believe me, you're shooting yourself in the knee. More suitable MCUs are cheap, and so are development boards (they literally start at 2$) sufficient to connect your ADC to them.


How do I time-stretch a signal from the frequency domain?

Are you talking about taking an audio recording where someone was speaking too fast, and you want to slow down the speech without changing the pitch? If so, then your original formulation of the problem is incorrect; you can’t do what you want by taking one long FFT of the entire recording and then processIng it.  The usual way this is solved is by using the phase vocoder, where you break the signal into overlapping blocks and transform each block.  This is a deep subject with lots of tricks to avoid artifacts.


The paper actually deals with $\pi/4$-QPSK, not with plain DQPSK. In $\pi/4$-QPSK, the symbols are always rotated by a minimum of $\pm \pi/4$. Rappaport's Wireless Communications textbook has this table:

\begin{align*}
11 \rightarrow & \, \pi/4 \\
01 \rightarrow & \, 3\pi/4 \\
00 \rightarrow & \, -3\pi/4 \\
10 \rightarrow & \, -\pi/4
\end{align*}

Unfortunately, the paper you point to does not specify what mapping it's using; the only clue is equation (2), where we see that the rotations are $n\pi/4$ for $n=1,3,5,7$.

Let's assume that $\text{00} \rightarrow 7\pi/4$. Then, in your case you would have $s_{n-1} = \sqrt{2}/2 + i\sqrt{2}/2$, and $s_n = e^{i7\pi/4}s_{n-1} = 1$. Then, according to the paper,

\begin{align*}
b_1 = & \text{Re}(s_n s_{n-1}^*) = \sqrt{2}/2 \\
b_2 = & \text{Im}(s_n s_{n-1}^*) = \sqrt{2}/2
\end{align*}

I worked out $\text{00} \rightarrow 7\pi/4$ as the only rotation that produces results that make sense. However, it's not often that I see a paper with so many typos and so poorly typeset, and I wouldn't trust it without working out the details by myself first.


Pulse detection against various noise floors

I cannot comment on your question, so I am writing an answer here instead.

Be aware that any time you mention 'threshold', you are implicitly making a statement that you are comfortable with those three probabilities:

Certain probability of false alarms, (PFA).
Certain probability of detection, (PD).
Certain probability of misses, (PM).

Put another way, every time you set a threshold, those three numbers will change accordingly. This is inescapable, in any detection problem. Typically one can set the probability of false-alarm that they are willing to 'live with', and go from there.

For your problem, you mention changing noise variance. One way to circumvent this is to Z-Standardize your data vector before attempting anything. This then centers your data, (de-meaning), and forces the variance to unity. This take care of the problem of you having to re-compute thresholds all the time.

With this pre-processing step out of the way, you can now focus on your Test-Statistic, and your Threshold. For example, a test-statistic can simply be each sample. Or, it can be the average of three samples. Or the average of 10 samples. Either way, you will be comparing your test-statistic, against your threshold. Again, you pick a threshold you "can live with", that is to say, one that gives you false-alarms at a rate that your application deems acceptable.

For example, if a doctor has a false alarm rate of 10% in informing you that you have a cold, that might be acceptable for this application, because it is not the end of the world, and the worst thing that can happen is you drink more orange juice than usual. If on the other hand you are designing a fraud detection system for your bank and you false alarm in detecting fraud 10% of the time, you are going to have some very, very angry customers swearing at your representatives all day for locking their credit cards.

Back to your application: It seems that those bits are the result of a demodulated PAM receiver system, so you probably already have a significant amount of averaging gain. This means that a test-statistic of one sample is good enough to go off of. All you need to do now, is set your threshold (of your standardized data), such that you false alarm as some rate you are comfortable with.

(Remember, that you are dealing with standardized data, so you need only set this threshold once, and then forget about it. As a practical note to assist you in determining the threshold, try to plot the histogram of your data without signal, after you have standardized it. This will be your 'noise PDF'. This can assist you in estimating the optimal threshold ala Neyman_Pearson detector.)


Extracting absolute temperature values from an Infrared video

There is no direct relation between the emitted intensity and temperature, as this depends on the body color and many other factors.

Anyway, you can get an approximation by calibration. For this  you use a reference object and heat it to known temperatures. You will draw a calibration curve (pixel intensities as a fuction of temperature) and possibly map a simple model (polynomial) onto it.


How to determine a quantization table for image compression with DCT?

The standard JPEG quantization table is a actually an 8 x 8 matrix of weights whose elements divide the corresponding elements of the 2D-DCT coefficient matrix of an 8 x 8 image block.

The division effectively reduces the magnitude of the DCT coefficients and most of them are forced to be zero, due to the integer divison rounding.

The standard table is determined to yield a minimum distortion under maximum compression applied to the most typical image (i.e., exhaustive average of large image data set) and typical viewer visual characteristics.

As a result of this quantization, high frequency coefficients are reduced as much as possible.

Note that the JPEG quality parameter provides a means of scaling the standard quantization table weights, but no control over individual elements.

You can also design and use a specific quantization table for your needs, according to the rate-distortion criteria you have in your mind. But this requires you to embed Q table information within the compressed data bitstream.


The FFT output will be bins evenly spaced from DC (bin 0) to 1 bin less than your sampling frequency. Thus your sampling frequency is at N for bins 0 to N-1, where N is the number of samples in your time domain sequence.

Within the first Nyquist zone of $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate, the bins from 0 to N/2-1 are the positive frequency components, while the bins from N/2 to N-1 are the negative components. Specifically bin N-1 is -1, bin N-2 is -2, etc so you could view the domain as a circle from 0 to N-1 and repeating again.

The reason your peaks are a mirror image is because your time domain signal is real, and for real signals the positive and negative frequencies have a conjugate symmetric relationship (same magnitude, opposite phase).

The reason you see a large response in bin 0 is because your time domain signal has a non-zero mean. Bin 0, being the DC term, is the mean of your signal (scaled by a factor of N). You likely don't care about DC, and would therefore benefit by removing the mean of your signal before taking the FFT.

The FFT is an algorithm that implements the DFT efficiently, which is as the name implies a discrete version of the Fourier Transform. I think it is helpful to see that the Fourier Transform is simply a correlation as we sweep frequency. (Where correlation is complex conjugate multiply and integrate, or equivalently discrete as complex conjugate multiply and accumulate). Thus for the DFT, each bin is a correlation of a single rotating phasor (your traditional sine wave is made of two equal and opposite rotating phasors as shown by Euler's Identity: $\frac{1}{T}(e^{j\omega t}+e^{-j\omega t}$).  So for the DFT, given as: $$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j(k\omega_0) n}$$

So we see that for the first bin (k=0), we just have N times the mean of x. If x has a DC value, we would get a strong correlation for that calculation. Also observe when k =1, which corresponds to the next bin, the exponential is spinnning at its lowest possible rate $\omega_0$ corresponding to one cycle over your time domain sequence. So in that calculation we are computing the correlation of the lowest possible frequency that is spinning in the opposite direction ($e^{+j\omega t}$): such a signal in our x(t), after this computation would have the same result as bin 0: we would get N times the mean after having "derotated" the signal that was present there. Thus every bin goes through the same "averaging filter" which is a Sinc function in frequency (not a great filter), so our output at any given bin will be sensitive to some degree of frequencies everywhere in our spectrum (and as we see with the Sinc function there will be discrete frequencies that are at null positions and therefore will be locations where there is no "leakage" into other bins).

This should give you enough insight into some important fundamental points to dig further into the DFT/FFT. In addition here are some other related posts I think will help you:

Frequency measurement from acoustic emission sensor data(VDC RMS)

What proportion of a padded FFT should be actual values

Signal Processing/FFT gives very high magnitudes for low frequencies

You will want to read up a little more about windowing too; the links above and other posts (by search) include more information on that. Ultimately fred harris' paper has been a great resource for me on that topic.


Unsure if my unit step function is correct in relation to a diagram

If the diagram was given to you than your function is incorrect. First of all you don't need the small brackets for the time. If you look at the diagram you can see that you can disassemble it into smaller parts:
Part I: from 0 to 1 you have a curve which starts at (0,0) and ends at (1,3). If you plug in these $(x,y)$ values into the equation $y=kx+n$ you will get $n=0, k=3$ which means that the curve function is equal to $y=3x \Rightarrow f(t)=3t$. Since the interval is from 0 to 1 we get $u(t)-u(t-1)$ and multiplied by the curve function $3t[u(t)-u(t-1)]$.
Part II: from 1 to 3 you have a rectangular pulse with an amplitude of 3 which is $3[u(t-1)-u(t-3)]$.
Part III: from 3 to 4 you have again a curve whose function you retrieve the same way as in part I. You will get $(9-2t)[u(t-3)-u(t-4)]$.
Part IV: from 4 to 6 you have again a rectangular pulse but with an amplitude of 1 so you don't have to multiply the difference of the step functions because it's 1 by default. You will get this $u(t-4)-u(t-6)$.
Part V: Assuming that your signal equals to 0 for $t>7$ you again have a curve whose function you retrieve as in Part I and you get $(2t-11)[u(t-6)-u(t-7)]$.  Summing up all obtained parts your function would be:
$f(t)=3t[u(t)-u(t-1)]+3[u(t-1)-u(t-3)]+(9-2t)[u(t-3)-u(t-4)]+u(t-4)-u(t-6)+(2t-11)[u(t-6)-u(t-7)]$  If you use MATLAB here's the code and the plot:

t=0:0.01:8;
x = 3*t.*(stepfun(t,0)-stepfun(t,1))+3*(stepfun(t,1)-stepfun(t,3))+(9-2*t).*(stepfun(t,3)-stepfun(t,4))+stepfun(t,4)-stepfun(t,6)+(2*t-11).*(stepfun(t,6)-stepfun(t,7));
plot(t,x,'r');
axis([0 8 0 4])
grid on;


You wouldn't normally do that because it isn't worth the cost (hardware.)

FIR filters get computationally very expensive if you need sharp cutoffs.

Digital IIR filters that provide the same cutoff as a long FIR are computationally very cheap.

An analog filter will be equivalent to a digital IIR.

If an analog filter can provide the performance you need, then you are better off implementing a digital IIR filter.  It saves you the DAC/ADC parts and it saves you the mess of an additional analog section with its inherent instabilities and potential to pickup noise.

If the digital IIR that can provide the performance you need is unstable, then it is very likely that the equivalent analog filter will also be very difficult to make work correctly.


Lowest frequency within pitch detection in polyponic sound

The lowest pitch in a chord is likely to have a large number of its possible overtone frequencies match the actual overtones of other notes in a harmonious chord, leading to a good cepstral magnitude, which looks for spectral periodicity (overtone series's).

However a cepstrum might also find more overtone sequence periodicities at pitches lower than the lowest pitch in a chord.


Having Glitches trying to use a Sinc Function Interpolation on a Signal, by chunks of 1024 samples

Ah, that's interpolatef (it's in vector/src/vector_types/time_freq/interpolation.rs).

I'm not fluent in Rust, but I think it calculates the interpolation by:

Taking in a function (a real impulse response), and generating sampled vectors of that, with fractional offsets.
It uses that polyphase filterbank to do the arbitrary resampling

You could hence just "prepend" every chunk with conv_len samples, and then "throw away" the right amount of samples from the beginning of the result.

However, that might be hard: Since this function allows for arbitrary real resampling ratios, conv_length * interpolation_factor might not be an integer.

So, the easiest thing would probably to restrict yourself to rational resampling ratios, here, and always transport a number of samples $n \ge \text{conv_len}$ with $\mathbb N \ni n = \frac m{\text{conv_len}}$ so that you can always throw away an integer amount of samples.


Aliasing due to the Convolution of Gaussian Functions

Short answer: For the convolution of two Gaussian sequences there will be aliasing.

If the sequence produced as the result of linear convolution has infinite domain of support, then there will always be aliasing in the Fourier domain implementation of the linear convolution using the circular convolution implied by the DFT (discrete Fouerier transform) which is used to represent the theoretical DTFT that provides the theorem in your question.

Therefore, if you want an alias free implementation of the theorem in your question, then the resulting signal (result of the true linear convolution) must have a finite length.

That being said, for many practical cases, the amount of distortion due to aliasing will be small, so small to be lost behind any ADC step size or even smaller than represantable by the floating point number system that's used to implement the algorithm. In such cases it can practically be considered to be a sufficient realization of the true lienar convoluion.


Dynamically changing cut-off and sampling frequency of a digital filter

There are a few ways that you could do run-time dynamic filters.

With a program like Matlab or Octave, but I'm guessing that they aren't available to your application.  There might be DSP libraries that have functions for generating filters though.
Calculating from scratch a windowed sinc function.
You could create somewhat crude low-pass filters based off of your starting filter.  For instance, if you wanted to halve the cutoff frequency you would halve the filter sample time.  You would thus create a new filter that had all of the original filter's samples and a sample in between those samples.  You could create the extra samples through splines.  This same approach could be used to create any fractional cutoff frequency change.  A more complicated example would be if you wanted the cutoff frequency to be $\frac{9}{8}f_c$.  If we say that your original filter's samples are at time 0, 1, 2, ..., $n$, then the new filter's samples would be at time 0, $\frac{9}{8}$, $\frac{18}{8}$, ...


Not sure about audio plugins but what you want to do is certaininly doable by using filters (delay line in your case in the subbands. Split the band using using a Perfect Reconstruction Analysis filter bank. Now plugin the delay line and then combine two subbands using synthesis filter bank. Note that you need to be extra careful while plugging in the filter between analysis and synthesis filter bank. You can check some of these papers for ideas on how to do subband filtering.

https://dl.acm.org/citation.cfm?id=2817121

http://ieeexplore.ieee.org/abstract/document/1661283/


How to calculate step and impulse response at multiple sampling rates

This is not the full answer you want, but the problem is that you can't permute sampling and convolution. Let's say that the $D(\cdot)$ operator takes a continuous signal and returns its sampled version. What you're doing is

$$y'(n) = D\left(e^{-at}u(t)\right) * D\left(x(t)\right)$$

What you want to get is

$$D\left(y(t)\right) = D\left(e^{-at}u(t) * x(t)\right) \neq y'(n).$$

It is in general impossible to do what you're trying to do, but I'm sure that good approximations exist. I'm not an expert on what they are or how to use them.


What windowed Sinc might be closest to a Parks–McClellan low-pass filter?

none.

what you're asking for is what window, which multiplies a sinc() function, will get you an optimal impulse response kernel, like you might get from Parks-McClellan.  well, to get that window, you must undo the multiplication.

the problem is that you cannot divide by zero.  any of these "sinc-like" functions that come out of P-McC will pass through 0.  but not necessarily at the same place that the sinc() function goes through zero.


Performance evaluation in terms of channel MSE and BER

From the simulation section of the paper:

  The SER for each Monte Carlo simulation is calculated as follows: only
  when the activity state and the transmitted symbol are both detected
  correctly for a certain user, the corresponding user’s information is
  considered to be correct, and all symbols transmitted by the user
  whose activity state is mistakenly detected are considered to be
  incorrect. To provide a benchmark for performance comparison, we
  consider the oracle LS by assuming the actual active user support set
  known at the BS.

  However, the SER performance gap of the proposed joint CE and MUD
  framework and the oracle LS always remains. This is mainly due to the
  CS-based sparse signal reconstruction algorithm with non-orthogonal
  training matrix suffers significantly, so that the data detection
  error is the main cost of user activity detection.

The non-oracle receiver has to figure out who is transmitting and hence can do additional errors. The oracle knows this information already. That's why you have the difference in the SER.


What characterizies 'causality' for a finite FFT?

I am interested in this question and since no-one has answered yet I will offer my further thoughts (up for debate or confirmation):

I agree with the OP's perspective that the Discrete Fourier Transform (DFT) fundamentally represents when cast to domains that extend to $\pm \infty$ periodic sequences in time (and in frequency) and therefore is fundamentally non-causal due to the time periodicity specifically.

However, to utilize these transforms with context of causality, we can consider the cases that best match an approximation of the Continuous Time Fourier Transform for the same causal representation of our sampled system. Such cases are the ones with no time aliasing, which would occur when $x[n]=0$ for $n\ge N/2$ given $N$ samples with $n \in [0, N-1]$. When we can restrict ourselves to causal sequences specifically (meaning when we know the underlying process is causal) then this restriction can be extended to be $x[n]=0$ for $n\ge N$ given $N$ samples with $n \in [0, N-1]$. This means basically that we know the time duration of the data capture exceeds the expected response time of the underlying continuous time process.

Further in frequency for such an equivalently causal time domain waveform the phase would be going increasingly negative versus frequency indicating the delay that is indeed properly being modeled.

It is the counter case, when for a known causal sequence $x[n]$ has significant non-zero values in the upper-half of the time domain samples right up to the $N-1$ boundary that we cannot be reasonably assured we are not observing the effects of time domain aliasing. We are no longer able to uniquely distinguish between causal and non-causal sequences, or more specifically and practically when we know the sequence is causal; we are no longer able to distinguish from underlying continuous-time responses that went beyond the sample duration of the sequence due to time domain aliasing (so specifically in the OP's case there is no way to know if the solution is representing longer response times than the samples are providing).

This is not very different from what I would typically do when evaluating results in the frequency domain: if I observe strong spectral content right up to the Nyquist boundary I am not assured that my sampling rate is high enough (or filtering is not tight enough) as frequency domain aliasing is likely occurring. While if the spectrum rolls off to sufficiently low levels prior to the Nyquist boundary (and if I can assume the system in question would not have higher spectral content at further offsets) then I am reasonably assured that the digital representation accurately represents the continuous time spectrum. This obviously requires knowledge of the system and what the likely spectral content is beyond the sampling rate, but is very analogous to the time domain challenge presented here.

For further details on matching an underlying CTFT with the DFT scroll down to the "Exact CTFT Result from DFT" section of my answer at this post: Why do we have to rearrange a vector and shift the zero point to the first index, in preparation for an FFT?


Differences between Python and MATLAB filtfilt function

According to the newsgroup https://mail.python.org/pipermail/scipy-user/2014-April/035648.html,

  The difference is in the default padding length.  In matlab's
  filtfilt, it is 3*(max(len(a), len(b)) - 1), and in scipy's filtfilt,
  it is 3*max(len(a), len(b)).

Therefore, I modified the padlen when invoking filtfilt in Python. And the results between MATLAB and Python will be the same.

MATLAB code

input_signal = [1,2,3,4,5,1,2,3,4,5,1,2,3,4,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5]
[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal)

Python with SciPy code

from scipy.signal import butter, filtfilt
input_signal = [1,2,3,4,5,1,2,3,4,5,1,2,3,4,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5]
[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal, padtype = 'odd', padlen=3*(max(len(b),len(a))-1))


Can a wavelet transform give time dependent phase of sinusoids in signal

It is indeed possible (up to a period), with complex "continuous" wavelet transforms, provide your mother wavelet ooscillate sufficiently. You can look at the scalogram modulus and its phase at the appropriate scales, as illustrated from Matlab: Wavelet coherence:

Alternatively, you can convert your signal into an analytic one (using the Hilbert transform) and perform a real wavelet analysis.

Finally, you might dig into other linear (short-term Fourier transform) or bilinear time-frequency (Wgner-Ville, Choi-Williams) or time-scale tools.


How to interpret the effect of different windows in short time fourier transform?

A widnow $w[n]$ truncates and weights (tapers) an input signal $x[n]$, to produce $v[n] = x[n]. w[n]$., for subsequent spectral analysis of $x[n]$. A windows's effect on the input signal's true spectrum $X(e^{j\omega})$ is described by a convolution of $X(e^{j\omega})$ with $W(e^{j\omega})$ (window's Fourier transform);
$$V(e^{j\omega}) = \frac{1}{2\pi} \int_{-\pi}^{\pi} W(e^{j\theta})X(e^{j (\omega -\theta)}) d\theta $$
due to this convolution, two main effects are observed on $X(e^{j\omega})$.
1 - smoothing (smearing) of $X(e^{j\omega})$ due to the main lobe width of $W(e^{j\omega})$,  which results in a loss of spectral resolution.
2- spectral leakage due to peak side lobe of $W(e^{j\omega})$, which results in a loss of weak components shadowed by nearby stronger ones.
Main lobe width of any window type is primarily determined by its length. Increasing the length of any window will therefore decrease its main lobe width (hence increase its spectral resolution capability)
A rectangular window has the narrowest main lobe width and the highest peak side lobe, compared to all other windows. And the remainig window types perform a tradeoff between mainlobe width and the peak side lobes.
Peak side lobe is primarily determined by windows's shape. So by changing the shape (type) of the window you adjust it's leakage amount.


Compression ratio is $$ \frac{N_{uncomp}^b}{N_{comp}^b}$$
where ${N_{comp}^b}$ is the total number bits requried to represent those $3$ pixels $A$,$B$, and $C$ while ${N_{uncomp}^b}$ is the total number of bits required when they are not compressed.

Based on your claims, assuming that, those three pixels had $8$ levels with fixed length coding (FLC) of $3$ bits per pixel then; ${N_{uncomp}^b} = 3 \times 3 = 9$ bits in total. And again based on your variable length code (VLC) of the lossless coding scheme which uses $1$ bit for $A$ and $C$ and $4$ bits for $B$ you have ${N_{comp}^b} = 1 + 4 + 1 = 6$ bits in the compressed case, hence the ratio is $$ \frac{N_{uncomp}^b}{N_{comp}^b} = \frac{9}{6} = \frac{3}{2} $$

which means a $1.5 \times$ reduction in the number of bits required to represent the same source compared to its uncompressed form.

Note that you can also use the reciprocal of that ratio to indicate a similar metric, which would be $2/3$ in this case indicating the size that compressed form will have compared to the original.


Pyramid Match Kernel: How to fit grid around data points?

No, that's not it. Look at the last sentence of your quote:

  ... which may be enforced by scaling the data to some precision and
  truncating to integer values.

So this is exactly what you do. Take all your feature vectors, multiply them by 1000 (or some other factor), and truncate to integer. Then the distance between any two unique feature vectors will have to be >= 1. If you have a new vector y, then you multiply it by the same factor and truncate.


How to pick a sampling frequency for Matlab fft

The sensor is analogue and being input to Arduino for analogue to digital conversion. Arduino has a function Micros() which Returns the number of microseconds since the Arduino board began running the current program, I suggest using this to accurately save the time of analogue read on Arduino. Calculate the difference as shown in the code below and send it through with the sensor reading to Matlab. Then just average the values of the time difference calculated by t2-t1 which will be almost constant, and you will have as accurate as possible the sampling time/Sampling frequency used and so you can very accurately do an FFT analysis using this information. Hope this will be helpful.

//Define Analogue Ports for Sensor input
const int xpin1 = A2; // x-axis of the accelerometer 1

//Define variable t1 for timing
unsigned long t1 = micros();

void loop()
{
  int x1 = analogRead(xpin1);  //read x from 1

  Serial.print(x1);  //print x1 value on serial monitor
  Serial.print("\t");

  unsigned long t2 = micros();
  Serial.print(t2-t1);
  Serial.print("\t");
  t1 = t2;
}


I'm not in the image processing filed, but I would say:
Yes, Gaussian noise means that each of your pixel is assumed to be produced by a stochastic process. Hence you can for instance think of it as a "real" intensity corrupted by an additive normally distributed white noise. Normally distributed means also that if you take $n$ observations of the same pixel you most likely will end up with $n$ wrong measurements. I think this is actually the worst kind of noise since given any pixel, then any other pixel inside a given neighbourhood will be corrupted, hence telling which observation is to trust more is difficult. The way I would go is probably (low pass) filtering.

You can think to the impulsive noise again as an additive noise potentially present in all the pixels, but potentially also localized in only a subset of pixels or observations (either it is or it is not). I personally would not go for averaging or low pass filtering since it will never erase completely the noise.
Assume indeed to have $n$ images of the same scene. Assume that for a given pixel $p$ whose "actual" value is $p=0$ you measure it as $255$ in one image and as $0$ in $n$ images. If you take the mean then you will chose $p = 255/n$ as the actual value. Infact a low pass filter is a process with infinite impulse response and it tends to zero only asymptotically.
In such case for me it looks more effective to take the median, over the $n$ samples. With the median you obtain that you estimate of the actual values of $p$ is exactly $0$.
You can also take the most-frequent value (the mode) and it will perform better than the mean. In general you can also use a median or a mode inside the same image by taking them in a neighbourhood of a given pixel. Infact for example, take a $10\times 10$ sub-image and assume it corresponds to a white background. Assume that 30 over the 100 pixels are black due to the impulsive noise, then you can substitute to each pixel the median over the entire sub-image obtaining a completely white patch. A mean will lead to a 30% grey (76 / 255)


Filters are usually real, not complex, but occasionally they are complex.  What do you do with the complex impulse response?  Convolve it with your data, just like you would with a real impulse response.  This will likely result in a complex output.  What you do with the imaginary part at that point depends on the situation.  Sometimes it makes sense to leave the data complex, sometimes it makes sense to get rid of the imaginary part.

The imaginary part of the filter would particularly affect the output phase, so the fact that the output phase is not what you expected when you only used the real part of the impulse response makes sense.


First of all - I'm not familiar with Python and didn't check your code.
1a) For FFTs it is the usual that the frequency resolution equals the inverse length of your data record.
1b) As far as I understand your code (I may be wrong here), you take the FFT of the step response. The Bode diagram usually shows the FFT of the impulse response. Take care - the result may not be what you expect.
1c) In practice (due to transients, noise,...) the usual approach to measure the FRF is

excite the system with a broadband signal, estimate the spectral density of the input and the output signal, divide the output spectral density by the input spectral density. This approach is well covered in the literature (and called "spectral method" or something similar

exite the system with a sinuosoid, wait until the transient has died out, and read the magnitude and phase from your data. Repeat for every frequency where you want to know the FRF.

There are other approaches to directly measure the FRF, like the Local Polynomial Method (or Local Rational Method), TRIMM, and more - they are covered in (scientific) publications, together with a lot of beautiful maths. The approaches I mentioned are called non-parametric methods (there are also parametric methods, but that may lead to far away for your use case)

The phase jump at Nyquist frequency is to be expected. The FFT spectrum of a rational signal will be symmetric to the DC-frequency - the magnite is axisymmetric, the phase point symmetric. For some FFTs the "negative" and "positive" frequencies are wrapped (this may depend on your FFT implementation). The result is a "jump" the phase.

The input signal should excite the entire frequency range of interest - check the spectral density function of the input signal (i.e. the step).


Finding inverse z transform for two sided ROC?

Since you have two poles, you have to decompose $0.3 <|z|<1$ into two ROCs, $R_A$ and $R_B$ such that: $$R_A \cap R_B \subseteq \{0.3 <|z| < 1\}$$

Your poles are $z=-0.3$, $z=-1$, so you have two choices for each ROC:

$|z| > 1$
$|z| < 1$

and

$|z| > 0.3$
$|z| < 0.3$

The requested ROC is given by $R_A \cap R_B = \{|z| < 1\}\cap\{|z| > 0.3\} = \{0.3 < |z| < 1\}$

In this case (and using your table),

$$y[]=-2.86(−1)^{}[--1]−1.86(−0.3)^{}[]$$


Bluetooth locating with RSSI and only one receiver

If the environment does not contain a lot of 2.4 GHz radio reflecting objects, you can try shielding one side of the BTLE receiver's antenna (it might be tiny) with a large conductive object, the rotating the shield and antenna configuration.   The difference between RSSI's will be greatest when the shield is between the source and receiver antenna compared to the opposite orientation.


Is the fundamental basis function of the DFT not an intregral number of cycles?

First, a small correction:
The basis functions of all Fourier transforms are of the the type $e^{2j\pi k\frac nN}$, which can also be interpreted as $\cos(2\pi k \frac nN) + i\cdot\sin(2\pi k \frac nN)$, but never only real-values sines or cosines; that'd lead to large problems on the algebraic structures around. But that really doesn't matter here:

You're absolutely right, the "completing" sample of the oscillation is always missing – but that's absolutely logical, if you think about it – when you periodically continue the signal you're transforming, you'd get twice the same value (the "start" value) after each other, and that mustn't happen.

Your picture is an excellent illustration of this:

The oscillation is observed for exactly two periods, and each period is 16 samples long; that means that for the $n$th sample, the $(n+16)$th sample must have the same value. For example, the $n=16$th has the same value as the $(16+16)=32$nd sample – which, by principle of 32-periodicity is the same as the $0$th sample.