Least Square Error Estimation: Conditions for $ (A^TA)^{ -1} = A^{-1}(A^T)^{-1} $?

First of all, the minimum norm least square solution is $A^+b$, where $A^+$ is the pseudoinverse. Only when the left inverse $A_L^{-1}$ (or right inverse $A_R^{-1}$) exists, you have $A^+=A_L^{-1}=(A^TA)^{-1}A^T$ (or $A^+=A_R^{-1}=A^T(AA^T)^{-1}$).

Likewise, the projection matrix onto the column space is $P=AA^+$, or $P=AA_R^{-1}=AA^T(AA^T)^{-1}=I$ if $A_R^{-1}$ exists. In other words, if A has full row rank, then the projector $P$ is indeed an identity matrix. If only left inverse exists, then $P=AA_L^{-1}=A(A^TA)^{-1}A^T$ is not an identity matrix.

Finally, $(AB)^{-1}=A^{-1}B^{-1}$ is valid only when $A$ and $B$ are square matrices.


Because that's way more efficient at encoding a lot of typical scenes.
Imagine this: you're filming a complex scene with lots of different objects. The camera shifts ten pixels to the right; with your approach, you'll have to save a relatively large difference for every pixel on average, because the pixels 10 px to the left aren't very similar.
When you first encode the fact that this is just a shift by 10 px, that needs nearly no space at all - and then you'll just need to encode how the scene actually changed.


How to classify audio files based on their ZCR, energy, and enthropy of energy features?

This is a machine learning problem, and it can be classified using any good classification algorithm like 1) Linear Classifiers: Logistic Regression, 2) Naive Bayes Classifier, 3) Nearest Neighbor, 4) Support Vector Machines, 5) Decision Trees, 6) Boosted Trees, 7) Random Forest, 8) Neural Networks. Matlab has in-built classifiers, and you can train your model very easily.

As you have already found out the features you should get some good classification results using any of these classifiers. All of the features that you mentioned are time domain-based. Also, you can incorporate some frequency-based features like 1) Spectral Centroid and Spread (C.G and deviation from CG), 2) Spectral Entropy, 3) Spectral Flux (Measure of Spectral Change over successive frames), 4) Spectral Roll-off (Spectral shape descriptor of an audio signal) which may increase the classification accuracy.

For more audio signal features you can go through this book "Giannakopoulos, Theodoros, and Aggelos Pikrakis. Introduction to audio analysis: a MATLAB® approach. Academic Press, 2014."


Where is the mistake in the convolution?

$3^{rd}$ term in your final answer should have $6 \le n \le 10$, instead of $0 \le n \le 4$. And, the limit of summation will go from $n-4$ to $6$.
$$\sum^{6}_{k = n-4} \alpha^k \ ,  6 \le n \le 10$$
$0 \le n \le 4$ case has already been covered by the $1^{st}$ term of final answer.


Frequency response of a microphone using a sine sweep

Say $\mathbf{h}_1$ and $\mathbf{h}_2$ are the impulse response coefficients of mic1 and mic2, respectively, where $\mathbf{h}_1$ is known. Let $\mathbf{x}$ represent the input signal. Then, we can write the output of mic1 as

\begin{align}
\mathbf{b} = \mathbf{X} \mathbf{h}_1
\end{align}
where $\mathbf{X}$ is a Toeplitz matrix. If you want to find the impulse response coefficients of mic2, then you need to solve the following optimization problem
\begin{align}
\arg \min_{\mathbf{h}_2} ||\mathbf{X} \mathbf{h}_2 - \mathbf{b}||^2_2
\end{align}

The solution to the above problem is a simple least-squares solution. In the above formulation, I tried to express the convolution operation as a matrix multiplication. Find the filter coefficient response and express the transfer function in terms of $H_2(z)$. Once you know the transfer function, it is easy to determine the magnitude and phase response.


Optimum Filter Signal Detection for Non AWGN Channels

AWGN is a specific noise model that is to varying physical circumstances, of varying utility.

Not AWGN is not a specific noise model, but you can look at it as detection not dependent on Gaussian assumptions.  There is no one kind of optimal detector for every kind of non Gaussian noise.

A good book is

https://books.google.com/books/about/Signal_Detection_in_Non_Gaussian_Noise.html?id=tvpSAAAAMAAJ&hl=en

Multiplicative noise is another class as well.


There are 2 properties which are important to know about your noise:

Its Self Correlation
What usually called Auto Correlation. This will tell you the relation between 2 noise samples which are apart (On time) from each other (I'm assuming the signal is stationary, namely its properties doesn't change over time).
Best thing to hope for is White Noise.
Its Distribution
Sometimes, in order to be able to do estimations of the data one should also know the distribution of the noise. The most wide spread model is Gaussian Noise as you wrote. Yet on Image Processing it could be Poisson Noise.

If you have access to "Pure" noise samples analyzing those properties of the noise is really easy.
For the first you just calculate the Auto Correlation Function.
For testing Gaussianity / Normality Test (Reference - Tests of Gaussianity) things gets a little tricky.
As you can see in the Wikipedia Page of Normality there are many tests.
Though poor man's (Yet effective) test would be just looking at the Empiric Histogram and compare it to the model (Curve of Normal PDF based on the data Mean and Variance).


How do I find the convolution of the following signals?

Once you understand the intuition of convolution using the link in the comments, you can work on deriving the equations. These problems are a bit tedious at first and can cause confusion if you try to do too much at once. I'm uploading my hand written notes on how to break the problem down. The key is to think about what happens when there is no overlap, partial overlap when the box enters the triangle from the left, full overlap (occurs only when $t=6$), and partial overlap wen the box begins to exit the triangle.

So now we can check our work using MATLAB! We just need to be careful and understand that conv computes the convolution sum, not integral, but we can still use it to approximate. Here is the code to check your work:

%% Do the problem in MATLAB first
% Original time axis for x(t) and h(t)
t = linspace(-10, 10, 1000);

% Make x(t)
x = zeros( size(t) );
x(  t <= 5 & t >= 3  ) = linspace( 2, 0, length( find( t <= 5 & t >= 3 ) ) );

% Plot the triangle
subplot(4, 1, 1)
plot(t, x)

% Make h(t)
h = zeros( size(t) );
h( t >= 1 & t <= 3 ) = 2;

% Plot the box
subplot(4, 1, 2)
plot(t, h)

% Calculate the convolution and new time axis
y = conv(x, h);
ty = linspace(-20, 20, length(x) + length(h) - 1 );
dy = ty(2) - ty(1);

% Need to multiply by dy to get the same result as doing the integral by hand
% MATLAB conv() does discrete time sum for the convolution but we can
% approximate the integral by multiplying by the time step dy
y = y .* dy;

% Plot the convolution result
subplot(4, 1, 3)
plot(ty, y)
ylim([0 5])

%% Verify using "by-hand" equations
% Make the plot using the derived equations
t = linspace(-10, 10, 100);
y = zeros( size(t) );

tt = t( 4 <= t & t <= 6 );
y( 4 <= t & t <= 6 ) = 10*(tt - 1) - (tt - 1).^2 - 21;

tt = t( 6 < t & t <= 8 );
y( 6 <= t & t <= 8 ) = (tt - 8).^2;

subplot(4, 1, 4)
plot(t, y)
xlim([-20 20])
ylim([0 5])
```


Downsampling in MATLAB - why scaling filter coefficients?

This is done to achieve unity gain, i.e. to have the same energy in the downsampled and original signal.

The sum of filter coefficients sum(h) is the filter gain, by dividing the filter by it, it is set to 1 (unity).
It is multiplied by p to compensate for the missing samples, again to achieve the same engery as in the original signal.


Training own classifiers using trainCascadeObjectDetector for face detection

If you want to train the classifiers with your own database, you will only need 'trainCascadeObjectDetector' function and feed your images into the proper arguments (Positive, negative images). The output classifier will be in your 'outputXMLFilename' as in  traincascadeobjectdetector

trainingImageLabeler is helpful function for classifying positive/negative image. Try this function later when you work with larger database.


Matlab: How to apply constant modulus algorithm in wireless communication

Ok, your code is very difficult to understand, but I think that we can answer your question.

There is a group of processing blocks that every communication systems have. Normally, they appear in the following order in almost every simulation:

Bit Stream -> Modulator -> Channel -> Equalization -> Demodulator -> Received Bit Stream

Note that I suppressed some blocks like pulse shapping, matched filter, FEC, and others. This is the most simple scheme. The transmitter modulate the signal, and sends it through a channel. In the receiver, the samples needs to be equalized in order to perform the detection and demodulation properly.

So, you can conclude that the CMA blind Equalization is somewhere between the channel and the detection.

In your code I only saw the test for the Walsh Code. I am not seeing the complete system implemented (i didn't see any modulator, for example), but I assume that now you can figure out where to include the CMA equalizer.

Let me know if you need more help.


Quadrature Amplitude Modulation: I and Q Same

It doesn't really affect anything at all. If you modulate a signal on the I channel only, then it would look something like:

$$
x_1(t) = m(t) \cos(2\pi f_c t)
$$

If it were on the Q channel only, then it would look like:

$$
x_2(t) = m(t) \sin(2\pi f_c t)
$$

$x_1(t)$ and $x_2(t)$ only differ by a phase shift of 90 degrees on the carrier. Likewise, if you were to split the signal evenly in power across the I and Q channels, it would look like:

$$
x_3(t) = \frac{1}{\sqrt{2}} m(t) \cos(2\pi f_c t) + \frac{1}{\sqrt{2}} m(t) \sin(2\pi f_c t)
$$

Again, in the complex baseband form that is often used for analyzing signals, this corresponds to just a phase shift of 45 degrees from the original signal $x_1(t)$. In most systems, you will have a random phase shift between the transmitter and receiver anyway; this should not have any effect on the performance of your receiver.


According to the scenario you formulated, you only have measurement noise and no input noise. You are correct that you might need a state space model. For example,

$$
\begin{bmatrix}\dot{x} \\ \ddot{x}\end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 0 & -\frac{b}{m}\end{bmatrix} \begin{bmatrix}x \\ \dot{x}\end{bmatrix} + \begin{bmatrix} 0 \\ 1 \end{bmatrix} u + \begin{bmatrix} 0 \\ 1 \end{bmatrix} v,
$$

$$
\begin{bmatrix}y_1 \\ y_2 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix}x \\ \dot{x}\end{bmatrix} + \begin{bmatrix}w_1 \\ w_2\end{bmatrix},
$$

with $u$ always equal to $g$, $v$ the input noise (which assumed to be zero), $y_1$ the measured position, $y_2$ the measured velocity and $w_1$ and $w_2$ zero mean Gaussian noise.

Since we did not take into account input noise, therefore you can also view it as a zero mean Gaussian noise with also zero (standard) deviation. When calculating the covariance matrices you can take the variance (standard deviation squared) on the diagonal, if all random variables are independent (which I assume $w_1$ and $w_2$ are). This would yield,

$$
Q = \mathbb{E}(v^2) = \begin{bmatrix}0 \end{bmatrix},
$$

$$
R = \mathbb{E}\left([w_1\ w_2]^T [w_1\ w_2]\right)= \begin{bmatrix}225 & 0 \\ 0 & 25 \end{bmatrix}.
$$

But depending on whether you want to implement a Kalman filter in continues or discrete time, you might first have to discretize the state space model.


Reconciling Continuous and Discrete Complex Domains

Simple answer: that's how the math works out.

Better answer: In the discrete domain, frequency is periodic with the sample rate. If you sample at 48 kHz, the range of frequencies that you can represent range from -24 kHz to +24 kHz. Let's say so you have a 100 kHz signal. You can sample this at 48 kHz but it would look exactly like a 4 kHz signal. In fact the two are indistinguishable, so it makes no sense to treat them differently.

That's why a circle makes a lot of sense here. Frequency is represented as the angle between the origin at any point on the circle. If you go around the circle an integer number of times, you end up exactly where you started.

So in the continuous domain 4 kHz and 100 kHz are different points. In the discrete domain (when sampled at 48 kHz) they are not.


Efficient way to get power levels from an SDR (RTL-SDR specifically)

As far as I understand your question, you want to estimate the PSD in a large bandwidth.

What you're currently doing is tuning to one part of that bandwidth, receive samples at a low sampling rate, calculate the power in that small bandwidth, retune, and repeat.

So, yes, using a larger sampling rate will obviously reduce the time you need to sense a large bandwidth: A larger sampling rate will give you more spectral information.

The DFT->square method is one classical spectrum estimator. So, yes, instead of tuning $N$ times to receive an $N$th of the overall bandwidth $B$ by sampling at a rate of $\frac BN$, you could just as well sample with a rate of $B$, and calculate $N$-point DFTs in the same time. Since sampling the same time at rate $B$ gives you $N$ times the samples of sampling at $\frac BN$, the information would be the same (Nyquist/Shannon).

You might want to look into other spectral estimators, too – but for now, I don't really see a big disadvantage in simply using the DFT with an appropriate Window – and given the fact that neither the analog nor the digital filters on the RTL dongle probably are overly fantastic, I think you can just stick with a Hann window and be fine.


How come FFT window give half size frequency bands when filter can seems to have an infinite precision?

The DFT is a tool to represent the spectrum of any signal. Usually for an N point DFT only N/2 output coefficients are required because for a real signal the DFT is symmetric about the centre.

No filter has infinite precision, neither infinite bandwidth nor numeric representation. A filter with enough precision can be chosen for an application to seem smooth enough.

The DFT though doesn't have anything to do with a filter. It simply represents the spectrum of the signal at it's input.


Rayleigh fading with frequency selective fading channel

First it is useful to understand what a Rayleigh and Ricean distribution is, and then we can more easily see how that statistic can be used to describe certain fading channels.

A zero-mean complex signal that has an identical Gaussian distribution on both the in-phase and quadrature components will have a magnitude that is Rayleigh distributed. A useful explanation I came across that explained this well is an Agilent (HP) app-note on spectrum analyzer noise.

https://literature.cdn.keysight.com/litweb/pdf/5966-4008E.pdf

And I will paste the relevant figure below, thanks the HP/Agilent/Keysight. The first plot shows typical samples of a complex zero mean Gaussian signal. If we drew the histogram of the distribution on the real and imaginary axis, both would be the familiar bell curve of a Gaussian distribution.

If the complex signal with an identical Gaussian distribution on both the in-phase and quadrature components was not of zero-mean, such as I depict in the figure below, then the amplitude distribution of this process would be Ricean.

So how does this apply to fading channels? The central limit theorem comes into play here. The central limit theorem states that the sum of a large number of independent random variables will have a distribution that tends toward a Gaussian distribution regardless of the distribution of the individual processes.

In the case of fading channels we are seeing the result of many independent signal paths from the transmitter reaching the receiver at different delays. The different arrival times results in each signal being an independent random process. Often in terrestrial land-mobile communications there is no particularly strong direct-path, but a distribution of many reflected paths resulting in a zero mean Gaussian complex signal at the receiver (once the carrier is removed). To the extent that the signal actually has this process (which means to the extent that there are many indirect paths and no particularly single dominant direct path), then the AMPLITUDE or related power level of this received will be well modeled by a Rayleigh distribution.

When there is a dominant signal path present, particularly stronger than the others, with the presence of multiple reflected paths, then the amplitude of the received signal will follow a Ricean distribution. This is more common in satellite communications for example.

The above occurs whether we are dealing with flat fading (the amplitude across the entire signal bandwidth changes as one), or frequency selective fading (the amplitude at particular frequencies within the channel bandwidth is fading independent of the amplitude of other frequencies within the channel bandwidth.

Frequency Selective or Flat Fading

To be frequency selective fading or flat fading is independent of the distribution of the signal strength being Rayleigh. That only has to do with delay spread. (See figure below).

The Delay Spread of the channel is what determines the frequency selectivity of the fading (regardless if it is Rayleigh or not). Typical delay spreads for the macrocellular environment are shown in the figure below:


Understanding the Threshold Process in Harris Corner and Edge Detector

After doing the low and high thresholds to the eigenvalues, there will be 2 types of corners:

weak corners which are close to the low thresholding value.
strong corners which are above the high thresholding value.


Sure.

Roughly speaking the PSD is simply the magnitude squared of the FFT. So any set FFT's that have the same magnitude but different phases will result in the same PSD. For example white noise and a delta impulse have the same PSD.

The following figure shows the impulse responses of a filter that has the same magnitude frequency response.


Additive disturbance in block diagram representation

The pump keeps the tubing pressurized, the desired input flow is tuned via the valve. However the net flow rate that feeds the tank is not this inflow, it is the difference between the inflow and outflow.
The outflow caused by the pump is modeled as a disturbance because, in the current example, the goal is to model the liquid's level inside the tank. The outflow maybe considered a byproduct of the pump's operation.
Although the above diagram describes the plant (the physical process of the hydraulic system), if the control of the plant was studied later a possible design goal could be to maintain a constant and controllable level of the liquid in the tank (this requires feedback control aka a closed loop) and/or reaching the steady state value(desirable level in the tank) fast.
In addition a pump's flow rate tends to decrease as the pressure inside the tubing increases.
Think about a garden's hose, which sinks water into a flower pot at some flow rate and the water drains out of the flower pot through the bottom. Then if the inflow is much larger than the outflow the soil's surface will be flooded not being able to properly absorb the incoming flow of water. On the other hand, if the outflow is close enough to the inflow then the soil will be able to absorb the incoming water flow, not resulting in the flower pot's water overflow.
The block $G_t(s)$ describes the tank's response therefore the inflow-outflow difference is its input. The system $G_v(s)$ is the valve's response to the user input.


Modeling ADC in Active Noise Cancellation

I am also working on ANC, with my understanding may be I could provide some clarity for you.

I could not see your point exactly on ADC modelling, but as far as I know this ADC transfer function will be covered in secondary path modelling.

  "However, there is also the transfer function from the noise source to
  the adaptive filter."

Do you mean modelling of the reference input microphone, as far as I can see only transfer function involved between noise source to adaptive filter is reference input microphone function.

The way you are using secondary path description is also a bit confusing to me.

But here is what I know from this paper page 945

Here P(z) is the unknown primary acoustic path between reference input sensor and error sensor.

Below text directly quoted from the reference provided

  Summing junction represents the acoustical superposition in the space from cancelling loudspeaker to error microphone, where the primary
  noise is combined with the output of the adaptive filter.

  Therefore, it is necessary to compensate for the secondary-path
  transfer function from y(n) to e(n) which includes the digital-to-analog (D/A)
  converter, reconstruction filter, power amplifier, loudspeaker,
  acoustic path from loudspeaker to error microphone, error microphone,
  pre-amplifier, anti-aliasing filter, and analog-to digital (A/D)
  converter.

The compensation as you might have already known is done by offline or online modelling techniques.

If you want to include your primary path modelling as well i.e reference microphone TF, you can do that or if you are considering feedback effect, then again in the feedback compensation this reference microphone transfer function is included.

Please correct me if I am wrong.
Thanks for bringing this up!


If you mean correlation at the origin (no signal shifts), then yes so long as all the overlapping frequencies the two signals are 90 degrees out of phase.

That is, for all positive $\omega$ that are overlapping, $h_1(\omega) = k_\omega h_2(\omega) e^{\pm i \tfrac{\pi}{2}}$, where $h_1$ and $h_2$ are the spectrums and $k_\omega$ is just a real scalar. The negative frequencies will be the conjugate of these.

If you mean correlation with signal shifts (cross-correlation - similar to convolution) then no, because if both signals have the same non-zero frequency then at some shift the sinusoids will not be 90 degrees out of phase, and therefore correlated.

More simply

Imagine a signal consisting of a sine wave and another consisting of a cosine wave of the same frequency. They will have overlapping spectrums (same frequency) but have zero correlation.

Add to the first signal a sinusoid of a certain frequency, and add to the second signal the same sinusoid but 90 degrees phase shifted.

Example

Barbara image with DC component removed, and Barbara image with DC removed and half spectrum multiplied by $i$ and the other half by $-i$ (90 degrees phase shift and conjugate for the negative frequencies)

Their spectrums overlap completely but they have zero correlation (at shift 0).


Why does it make sense to have noise of a different shape than the state vector?

Why does it make sense to model the noise as something which gets multiplied by a matrix $G$ in contrast to simply having a noise vector which has the same shape as the state vector?

The usual reason is because the model has some states that are not directly driven by the process noise. For example, suppose your states are:
$$
\mathbf{x}_k = \left [ \begin{array}{c}
d_k\\
v_k\\
a_k
\end{array} \right]
$$
where $d_k$ is displacement, $v_k$ is velocity, and $a_k$ is acceleration all at time index $k$.

Suppose, as with a car, you can only change the acceleration component of the state (by pressing on the accelerator or the brake [decelerator]):

$$
\mathbf{x}_{k+1} = \mathbf{A} \mathbf{x}_{k} +
\left [ \begin{array}{c}
0\\
0\\
1
\end{array} \right] v_k
$$

In this case,
$$
G = \left [ \begin{array}{c}
0\\
0\\
1
\end{array} \right]
$$


Calculating the Spectral Centroid of a Signal

The function you wrote is basically like calculating the center of mass of a body or for that matter the Mean Value of a Distribution.

If you remember from Probability, the Mean is given by:

$$ E \left[ y \right] = \sum_i {y}_{i} p \left( {y}_{i} \right) $$

Now, let's think of the Frequency Function (The DFT) as a probability function.
In order to do so we'll normalize it to have a sum of 1.

This is done by dividing each element by the sum of elements.
Then just apply the Mean operator to calculate the Mean of this "Distribution":

$$ {C}_{i} = \sum_{k} \underset{{y}_{i}}{\underbrace{k}} \underset{ p \left( {y}_{i} \right ) }{\underbrace{ \frac{{X}_{i} \left( k \right)}{\sum_{k} {X}_{i}
 \left( k \right) } }} $$

The result you'll get is exactly what you have written above.


To illustrate Justme's answer:

Discrete Cosine Transform (DCT) is a lossy

The DCT can't be a lossy algorithm, since there's an inverse operation that restores the original input exactly.

data compression algorithm

Also, it's not a compression algorithm: in- and output have the same size.
So, both your central statements are wrong :(

that is used in many compressed image and video formats,

yes.

including JPEG, MJPEG, DV and MPEG.

What's DV? And: MPEG is a huge family of video compression methods. There's not "the MPEG algorithm", there's dozens of different video compression standards under that name.

In this algorithm, special DCT coefficients are calculated for each 8x8 image block.

That applies to JPEG, and probably a few of the many MPEG codecs. It's not true for all MPEG compressors! (For example, MPEG-H Part II, also called H.265, uses blocks of 64×64, 32×32, or 16×16, 8×8 or 4×4, depending on the image content.)

Then the coefficients are quantized

And that's where the lossiness happens: it's not in the DCT, it's what happens to the output of it!

and the image block is represented as a matrix of these quantized coefficients.

Again, only applies to JPEG. Be clear about that!
The matrix is how it's often visualized. The matrix representation is actually non-existent in the memory or storage formats, usually. On the very contrary, the elements are typically stored in a zigzag diagonal ordering if we imagine the matrix. (that's because that puts values that tend to be correlated closer to each other, making the result better to compress using subsequently applied lossless methods like LZW, Huffmann.)

The algorithm utilizes the fact that human visual system does not distinguish small changes in color or intensity.

No, it uses the fact that human perception often cares about high-frequency changes less than about small changes in low-frequency components. Otherwise, the selective quantization would make no sense.

Hm, you've not written the greatest paragraph. But I think you understand a lot of things correctly. Be more careful in really being aware what exactly does what, and you'll be fine!


Songify seems use Prosodic Modification (Pitch, Time) for monophonic signals:

First you need get the pitch contour and pitch duration from musics you want to follow
(use one pitch track for monophonic signal like YIN,
AMDF, Auto-Correlation)
Split your imput signal in voiced/unvoiced
Apply Pitch Scale Modification in voiced parts to match your extracted pitch contour (you can use PSOLA technique here)
Apply Time stretch to match your extracted pitch duration (you can use PSOLA technique here)


Real-time frequency recognition in python

Here's a recent blog post that may be helpful:

https://jakevdp.github.io/blog/2013/08/28/understanding-the-fft/

You want to install numpy (and probably scipy, too), and read up on the Discrete Fourier Transform (DFT; this is what the FFT computes).  I'm not that up on Python media packages for the playback part. Python's wave module may be enough.

The algorithm could be as simple as computing the FFT of the last, say, 4096 samples, getting the index of the maximum absolute value, and multiplying by the bin width.  Or you can get fancier and use a Hanning window, peak interpolation, overlap processing, and exponential averaging.


Most efficient way to find single dominant frequency (without amplitude) in analog signal

One approach is to use a phase/frequency detector as the circuit determining error and to ensure that the initial loop bandwidth is wide enough to “see” both the weaker signal that would otherwise be a false acquisition as well as the strongest signal desired. The loop will lock onto the strongest signal within the loop bandwidth as long as it is 6 dB above the total power of the other tones (a phase vector diagram illustrates well the need for 6 dB margin as anything less causes abrupt 0/180 degree changes In phase such that a phase lock loop can’t converge to zero error). Once locked, the bandwidth can then be tightened for better SNR performance on the locked signal. This approach also leads to much faster acquisition time in general and is often referred to as “gear-shifting”.
If the signals are much wider than the acquisition bandwidth then more elaborate search and assessment methods are needed prior to setting the center frequency for acquisition. A separate frequency lock loop with wider bandwidth can be employed or other techniques that sweep over the entire bandwidth and access signal strength (such as the FFT or if processing time is available and resources scarce a linear frequency ramp can be employed). Note that the FFT approach can be significantly simplified since the number of bins (samples) needed in the FFT computation only need to be enough to ensure one bin for each step size equal to the loop acquisition bandwidth or less (so if the entire bandwidth of interest is ten times the loop bandwidth in its widest setting for acquisition or a real signal, then a 20 point FFT could be used to determine the location of the strongest bin).


Common Filter Types for Audio Applications

a lowpass with variable cutoff frequency and resonance.

RBJ has a cookbook on how to do this:

http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

Here the cutoff frequency is called w0, and the resonance is called Q:

LPF:        H(s) = 1 / (s^2 + s/Q + 1)

            b0 =  (1 - cos(w0))/2
            b1 =   1 - cos(w0)
            b2 =  (1 - cos(w0))/2
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha

If you want to vary these parameters continuously while passing signal, you are not supposed to use biquads (because updating the coefficients causes glitches?) and you should use state-variable filters instead:

http://www.earlevel.com/main/2003/03/02/the-digital-state-variable-filter/

This is also a 2nd-order lowpass with variable Q.

  Does anyone have insight into filter types commonly used to mimic the sound of traditional analog gear?

Generally any filter transformed to a digital IIR using the bilinear transform is an "analog-like" filter.  Bessel is an exception because the important property of a Bessel filter is its group delay, which is not preserved by the bilinear transform.  Bilinear transform attempts to preserve the frequency response primarily.


Magnitude of the Gradient in Frequency Domain

Remember that as part of computing the gradients of an image, you are really just convolving your image, with (two) spatial filters. This then gives you two 'gradient' images, which you can then combine.

For example, to compute the x-direction-gradient-image $G_x$, you would convolve your image with a filter $f_x$ that might look like this $\begin{bmatrix} 1 & 0 & -1\end{bmatrix}$. Similarly, to attain the y-direction-gradient-image,  $G_y$, you might convolve your image with a filter $f_y$ that looks like this $\begin{bmatrix} 1 \\ 0 \\ -1\end{bmatrix}$.

Then, your final gradient image would be $G = \sqrt{G_x^2 + G_y^2}$.

So far so good. However, you don't need to do the convolution pixel by pixel. Instead, you could simply take the 2D DFT of $f_x$, multiply it pixel by pixel with the 2D DFT of the image, and then take the 2D IDFT back into the spatial domain. This would then also give you $G_x$. (Do the same to find $G_y$).

So the answer is at least a partial yes, in that you can compute the x and y gradient images via the FFT algorithm, (an implementation of the DFT), and then simply treating each of them as an element in a vector, and finding the L2 norm.


Inverse Problem / Deconvolution with Pink Noise

Some Remarks:

2nd Order Knowledge
It seems you only have a knowledge about the Auto Correlation of your data (2nd Order) and not on its distribution. Hence methods you can apply are ones which minimizes only 2nd order functions of the noise (Such as MMSE).
Wiener Filter Minimizes the MMSE
As can be seen in the derivation of the Wiener Deconvolution Filter the only assumption made is independence from the signal. If that holds in your data this is a good solution. Pay attention that Wiener Filter model is that the input data is also stochastic and all data is wide sense stationary.
Richardson Lucy Deconvolution
Pay attention this is a Maximum Likelihood method. Namely it assumes the data is a parameter, not a random data. It also assumes the distribution of the noise and not only its second moment,
Decide the Cost Function
First thing you should decide what would you like to minimize according to your model. Does the MMSE makes sense? Maybe a Least Squares is good enough? May you come with other objective. This is the first step.


Thresholding a difference image with a very high kurtosis

The first thing that jumps to my mind is that you want to test if your data is significantly different from zero. Determining if it is positive or negative should be trivial.

I would set a threshold based on the standard deviation of the results. Anything greater than say $\pm2\sigma$ from the mean (zero) are considered to be significant changes. The suitable value of sigma to choose on how peaked your distribution is and how many and how big changes you expect.

This sort of approach is traditionally used to look for outliers, conventionally with $\pm3\sigma$, which is quite similar to what you want.


Properties of projected 3d points (to 2d), Principal component analysis (PCA)

The new 2D data is the projection of the 3D data onto the tangential plane of the underlying 3D shape.
This should allow you to answer most of your questions by a simple drawing.


Magnitude of intermodulation distortion products?

Certainly the amplitude of a harmonic can exceed the amplitude of the fundamental,
or the amplitude of an intermodulation product can exceed the amplitudes of the individual signals.   But this occurs when the system is so nonlinear that it would
make no sense to use it as an approximation to a linear system.
Indeed, the very
fact that such a nonlinear system is being used (intentionally)
means that it has some properties
that are useful in the application at hand.

Example: If the system is a squarer that produces output $x^2(t)$, then
with $x(t) = \cos(2\pi f_c t)$, the output is entirely the
second harmonic (and a DC term). Since the fundamental is entirely
absent from the output, the harmonic amplitude is infinitely larger
than the fundamental.  If two sinusoidal inputs of frequencies $f_1$ and
$f_2$ are added together before the
squaring operation, then the output has second harmonics at frequencies
$2f_1$ and $2f_2$, intermodulation distortion products at frequencies
$f_1+f_2$ and $|f_1-f_2|$ and no fundamentals at frequencies $f_1$ and
$f_2$.

Moral: Don't use a squarer and pretend it is a linear addition circuit
or linear amplfier. If you are knowingly using a squarer, it is for
some other useful properties (e.g. frequency doubler) of the circuit that
are important to you right now.


How to detect the presence one- and two-dimensional patterns in distributions of spatial points?

Perhaps you could perform blob detection to find the individual dots and then map their location to a 2d space, that preserves their spatial structure, through the use of some measure for example the center of mass of the blob.
Once you have mapped each blob to its coordinates in the 2d space you could perform clustering or fit a regressor since the patterns somewhat resemble piecewise polynomials, you could model them as second or third degree polynomials, and build your algorithm upon RANSAC, it is a probabilistic algorithm that fits a model while discarding the outliers, you could also introduce some extra hyperparameters:

the percent threshold of points that minimize the cost function (you could use the sum of squared errors as a base cost function to build upon) and thus belong in the pattern.
the number of piecewise relations if you opt in the piecewise polynomial approach. Then you could exclude valid points that fit one branch of the piece wise polynomials from next iterations until all such polynomials are found.

Another classic algorithm that can cluster points based on their geodesic structure rather than just the distance would be DBSCAN, choosing a suitable $\epsilon$ for your dataset as the points that belong to the pattern are rather close to each other. A more complex classical machine learning algorithm is ISOMAP, it can also cluster points that based on geodesic properties


Because $Tf{_s}$=1 and $N$ and $n$ are integers.  So the exponent becomes 2$\pi$ times an integer, say $K$. But $e{^{i\theta}} = \cos(\theta)+i\sin(\theta)$. Therefore $e{^{i2K\pi}} = 1$.


Hough Transform in Matlab without the built-in function

Check this out. It is an easy to follow implementation. You can compare your implementation against this. Here is a no loop version.

Another advice is to add comments in the code above so that it could be understood better.


What are advantages of having higher sampling rate of a signal?

Sampling at a higher frequency will give you more effective number of bits (ENOB), up to the limits of the spurious free dynamic range of the Analog to Digital Converter (ADC) you are using (as well as other factors such as the analog input bandwidth of the ADC). However there are some important aspects to understand when doing this that I will detail further.

This is due to the general nature of quantization noise, which under conditions of sampling a signal that is uncorrelated to the sampling clock is well approximated as a white (in frequency) uniform (in magnitude) noise distribution. Further, the Signal to Noise Ratio (SNR) of a full scale real sine-wave will be well approximated as:

$$SNR = 6.02 \text{ dB/bit} + 1.76 \text{dB}$$

For example, a perfect 12 bit ADC samping a full scale sine wave will have an SNR of $6.02\times 12+1.76 = 74$ dB.

By using a full scale sine wave, we establish a consistent reference line from which we can determine the total noise power due to quantization. Within reason, that noise power remains the same even as the sine wave amplitude is reduced, or when we use signals that are composites of multiple sine waves (meaning via the Fourier Series Expansion, any general signal).

This classic formula is derived from the uniform distribution of the quantization noise, as for any uniform distribution the variance is $\frac{A^2}{12}$, where A is the width of the distribution. This relationship and how we arrive at the formula above is detailed in the figure below, comparing the histogram and variance for a full-scale sine wave ($\sigma_s^2$), to the histogram and variance for the quantization noise ($\sigma_N^2$), where $\Delta$ is a quantization level and b is the number of bits. Therefore the sinewave has a peak to peak amplitude of $2^b\Delta$. You will see that taking the square root of the equation shown below for the variance of the sine wave  $\frac{(2^b\Delta)^2}{8}$ is the familiar $\frac{V_p}{\sqrt{2}}$ as the standard deviation of a sine wave at peak amplitude $V_p$.  Thus we have the variance of the signal divided by the variance of the noise as the SNR.

Further as mentioned earlier, this noise level due to quantization is well approximated as a white noise process when the sampling rate is uncorrelated to the input (which occurs with incommensurate sampling with a sufficient number of bits and the input signal is fast enough that it is spanning multiple quantization levels from sample to sample, and incommensurate sampling means sampling with a clock that is not an integer multiple relationship in frequency with the input). As a white noise process in our digital sampled spectrum, the quantization noise power will be spread evenly from a frequency of 0 (DC) to half the sampling rate ($f_s/2$) for a real signal, or $-f_s/2$ to $+f_s/2$ for a complex signal.  In a perfect ADC, the total variance due to quantization remains the same independent of the sampling rate (it is proportional to the magnitude of the quantization level, which is independent of sampling rate). To see this clearly, consider the standard deviation of a sine wave which we reminded ourselves earlier is $\frac{V_p}{\sqrt{2}}$; no matter how fast we sample it as long as we sample it sufficiently to meet Nyquist's criteria, the same standard deviation will result. Notice that it has nothing to do with the sampling rate itself. Similarly the standard deviation and variance of the quantization noise is independent of frequency, but as long as each sample of quantization noise is independent and uncorrelated from each previous sample, then the noise is a white noise process meaning that it is spread evenly across our digital frequency range. If we raise the sampling rate, the noise density goes down. If we subsequently filter since our bandwidth of interest is lower, the total noise will go down. Specifically if you filter away half the spectrum, the noise will go down by 2 (3 dB). Filter 1/4 of the spectrum and the noise goes down by 6 dB which is equivalent to gaining 1 more bit of precision! Thus the formula for SNR that accounts for oversampling is given as:

Actual ADC's in practice will have limitations including non-linearities, analog input bandwidth, aperture uncertainly etc that will limit how much we can oversample, and how many effective bits can be achieved. The analog input bandwidth will limit the maximum input frequency we can effectively sample. The non-linearities will lead to "spurs" which are correlated frequency tones that will not be spread out and therefore will not benefit from the same noise processing gain we saw earlier with the white quantization noise model. These spurs are quantified on ADC datasheets as the spurious-free dynamic range (SFDR).  In practice I refer to the SFDR and usually take advantage of oversampling until the predicted quantization noise is on level with the SFDR, at which point if the strongest spur happens to be in band, there will be no further increase in SNR. To detail further I would need to refer to the specific design in more detail.

All noise contributions are captured nicely in the effective number of bits (ENOB) specification also given on ADC data sheets. Basically the actual total ADC noise expected is quantified by reversing the SNR equation that I first gave to come up with the equivalent number of bits a perfect ADC would provide. It will always be less than the actual number of bits due to these degradation sources. Importantly, it will also go down as the sampling rate goes up so there will be a diminishing point of return from oversampling.

For example, consider an actual ADC which has a specified ENOB of 11.3 bits and SFDR of 83 dB at 100 MSPS sampling rate. 11.3 ENOB is an SNR of 69.8 dB (70 dB) for a full scale sine wave. The actual signal sampled will likely be at a lower input level so as not to clip, but by knowing the absolute power level of a full scale sinewave, we now know the absolute power level of the total ADC noise. If for example the full scale sine wave that results in the maximum SFDR and ENOB is +9 dBm (also note that this level with best performance is typically 1-3 dB lower than the actual full scale where a sine wave would start to clip!), then the total ADC noise power will be +9dBm-70 dB = -61 dBm.  Since the SFDR is 83 dB, then we can easily expect to gain up to that limit by oversampling (but not more if the spur is in our final band of interest). In order to achieve this 22 dB gain, the oversampling ratio N would need to be at least $N= 10^{\frac{83-61}{10}} = 158.5$  Therefore if our actual real signal bandwidth of interest was 50MHz/158.5 = 315.5 KHz, we could sample at 100 MHz and gain 22 dB or 3.7 additional bits from the oversampling, for a total ENOB of 11.3+ 3.7 = 15 bits.

As a final note, know that Sigma Delta ADC architectures use feedback and noise shaping to achieve a much better increase in number of bits from oversampling than what I described here of what can be achieved with traditional ADC's. We saw an increase of 3dB/octave (every time we doubled the frequency we gained 3 dB in SNR). A simple first order Sigma Delta ADC has a gain of 9dB/octave, while a 3rd order Sigma Delta has a gain of 21 dB/octave!  (Fifth order Sigma Delta's are not uncommmon!).

Also see related responses at

How do you simultaneously undersample and oversample?

Oversampling while maintaining noise PSD

How to choose FFT depth for ADC performance analysis (SINAD, ENOB)

How increasing the Signal to Quantization noise increases the resolution of ADC


Will the resulting spectrogram from this algorithm be 1D or 2D?

Spectrograms are two-dimensional, usually calculated from a one dimensional input vector.

As an example say that you want a spectrogram with k-pts of frequency domain span and M-pts of time-domain span.   With NO overlap an input vector of length N = k x M will fill the spectrogram.  Every k samples a k-pt DFT is performed and makes one of the M time-axis vectors on the spectrogram.

For a running spectrogram with an overlap of k/2, then N/2 input samples will be needed to create a k x M spectrogram, and the time scale will be cut in half.

Usually the logarithm only comes in to scale the amplitude of the output of the DFT.  i.e., each of the k complex samples in the DFT output becomes a scalar value computed from the log of the complex sample magnitude.


real time signal analysis of output from an audio device

I'm writing a test program that can detect pops/cracks/noise (unknown frequencies and duration) in the audio output of an audio device.

Pops/cracks are impulsive in nature and therefore quite broadband. noise is very broad but I would assume here that is anything additive to the sinusoid.

You can catch all three, given reasonable assumptions, using the same way that Total Harmonic Distortion (THD) is measured.

To do that, setup a very narrow band-stop filter, also known as a Notch filter, and apply it to the recorded signal directly. The cut-off frequency of the Notch filter should be exactly on the frequency of the sinusoid you are recording.

The signal at the output of the notch filter will contain everything except the sinusoid signal.

Which means that you would be able to capture cracks and pops, because they would spill across the spectrum as well as any additive noise to the signal.

To quantify the extent to which your signal is "noisy", you can then apply an envelope detector, which in your case is basically a moving average filter followed by a threshold.

The idea there is that you integrate the "residual signal" (what is left after you have removed the sinusoid) over a period of time depending on the sort of temporal resolution you want to achieve and then you threshold that quantity. Basically, if the total sum of the signal over 5 seconds (for instance) becomes larger than some threshold (similar to the $\epsilon$ you seem to already have used), then that time interval is marked as "noisy".

For a notch filter you can design one using plain simple pole-zero placement or use a Twin-T.

If the Q-factor of those filters is not satisfactory, in other words, if the Notch filters you obtain with these methods are not selective enough, then you can simply subtract the sinusoid at the input from the recorded signal. Of course, this has to be done "carefully" because any phase difference inserted in the middle would result in an offset of your $\epsilon$ and to cut a long story short, prior to subtracting the sinusoid, you would first have to make sure that its phase is aligned, in which case we are into the Phase Locked Loop land. The PLL will basically be generating a local sinusoid that is locked on to the phase and frequency of the incoming sinusoid and you can use that to subtract it from your recorded signal. This would make it probably the sharpest of notch filters.

You could of course sum all the bins except the bin that contains the sinusoid you are trying to measure but the challenging aspect there is going to be the resolution of the DFT.

Hope this helps.

EDIT:

Just as a note, I thought I would just say here that what you are doing there should be done in a multi-threaded application, so that the recording thread proceeds uninterrupted. When you write:

  Record 10 seconds of audio samples.
  Apply window, run FFT on the recorded samples.
  Record the next 10 seconds of samples.

then, it is good to keep in mind that a single thread application would miss part of "reality" between steps 1 and 3 when it was processing the frame (but not recording). You can minimise this "gap" error by processing smaller chunks of audio but within limits. To truly do it in as close to real-time as possible way, you would need to have one thread recording sound and stuffing a queue with frames and another thread on the other side of the queue popping frames for processing.


The rfft function returns complex values ordered in specific way. In order to retrieve the amplitude of your DFT you must take the absolute value of it. Easiest way to do that is to call y_a = np.abs(y).

I also noticed in your case that you are using fftfreq function to obtain corresponding frequency vector. In case of rfft you must use the rfftfreq. This will eliminate this nasty wrap-around on your plot. Also you should not forget about amplitude normalization by half of your window energy.

Here is the complete example:

from scipy.fftpack import rfftfreq, rfft
import numpy as np
import matplotlib.pyplot as plt

def spectrum(sig, t):
    f = rfftfreq(sig.size, d=t[1]-t[0])
    y = rfft(sig)

    return f, y

if __name__ == "__main__":
    N = 1000
    t = np.linspace(0, 4*np.pi, N)
    x = np.sin(10*np.pi*t)+np.random.randn(N)

    f, y = spectrum(x, t)

    plt.subplot(211)
    plt.plot(t,x)
    plt.grid(True)
    plt.subplot(212)
    y_a = np.abs(y)/N*2 # Scale accordingly
    plt.plot(f, y_a)
    plt.grid(True)
    plt.show()


Certain varieties of generalized means have been applied in signal processing, such as the Lehmer mean. This paper discusses several generalized means, mainly from a theoretical point of view.

Note, however, that all applications of generalized means deal with signals that are inherently non-negative. This is also the case in signal processing, where there are many applications with non-negative signals, such as frequency domain denoising, image processing, fuzzy logic, etc. The application of generalized means in image processing is discussed in this thesis.

I think that rather than looking for new definitions of generalized means for bipolar signals, it is probably more useful to define the signals properly such that generalized means can be applied.

A completely different alternative to generalized means might be general non-linear filters, such as Volterra filters, which don't have the restriction that the input signal must be non-negative.


How to identify causality, stability and ROC from the pole-zero plot?

For systems with more than one pole (with different radii), there are $3$ types of ROCs:

inside the circle with a radius corresponding to the smallest pole radius; the corresponding time-domain sequence is left-sided
outside the circle with a radius corresponding to the largest pole radius; the corresponding time-domain sequence is right-sided
in a ring defined by the radii $r_1$ and $r_2$ of any two poles with $r_1<r_2$, with no other pole at a radius in the interval $(r_1,r_2)$; the corresponding time-domain sequence is two-sided

In your example you have case $3$ above, because that's the only ROC that includes the unit circle, and, consequently, corresponds to a stable system. Both other possible ROCs (inside the circle with radius $0.5$ and outside the circle with radius $2$) do not include the unit circle, and, consequently, do not correspond to impulse responses of stable systems. Since a ROC equal to a ring between two poles corresponds to a two-sided sequence, the system can't be causal. A causal system has a ROC outside a circle (case $2$ in the list above).


Starting with your last question: to extrapolate you do just use your fitted coefficients and extend them to a greater $x$ range.

However, I think this approach is inherently flawed for data such as yours over any significant time period, especially when using a polynomial model. The problem is that any polynomial has a limited number of turning points. Once you pass these the data will tend to head off to $\pm \infty$ which is unlikely to model anything well.

Frequency domain analysis may be a possibility, but I expect you would run into similar problems. Another approach could be to use some sort of moving average or auto-regressive model.

Banks and similar financial institutions invest millions into modelling this sort of data so it definitely isn't easy.


Yes, all answers given by you are fine.

Causal: If ROC is outside the outermost pole

Left sided: If ROC is inside the innermost pole

Stable: If ROC contains the unit circle (marginally stable if it touches unit circle)


Is there any alternative characterization of sparsity of a signal in compressed sensing

Indeed, there are ways in which sparsity, or information content, may be estimated
at the acquisition device. The details, practicality, and actual usefulness of doing so
is debatable and heavily dependent upon the context in which it is applied. In the case
of imaging, one could determine areas of an image which are more or less compressibile
in a predetermined basis. For example, see "Saliency-Based Compressive Sampling for Image Signals" by Yu et al. In this case, the additional complexity requirements placed on the acquisition device provide marginal gains.

With regards to your questions about making determinations as to the usefulness of Compressed Sensing on a given signal at the time of acquisition: If the signal in question adheres to
any kind of model known a priori, Compressed Sensing is possible. Accurate recovery is
simply dependent upon ratio between the number of measurements taken and the degree to which
the sampled signal adheres to your model. If it is a bad model, you won't get past the phase
transition. If it is a good model, then you will be able to calculate an accurate
reconstruction of the original signal. Additionally, Compressed Sensing measurements are,
in general, future proofed. If you have a given number of measurements for a signal which
are insufficient in number of accurately recover the original signal using the model you
have today, then it is still possible to devise a better model tomorrow for which these
measurements are sufficient for accurate recovery.

Additional Note (edit):
The acquisition approach mentioned in your question sounded quite near to adaptive Compressed Sensing, so I thought the following might be of interest for the readers of this question. Recent results by Arias-Castro, Candes, and Davenport have shown that adaptive measurement strategies cannot, in theory, offer any significant gains over non-adaptive (i.e. blind) Compressed Sensing. I refer readers to their work, "On the Fundamental Limits of Adaptive Sensing" which should be appearing in the ITIT soon.


Trying to compare two speech samples that have been recorded by smartphones in pockets

The distance emitter-receptor and the presence of obstacles in the trajectory(clothes for example) will afect the recorded signal for sure. So some frequencies can be attenuated or just disappear.


Comparison Between Guided Filter (Edge Preserving Filter) and Gaussian Filter

Create  a synthetic image with abrupt change form black to white (Better yest from Dark Gray to Bright Gray).
Add Gaussian Noise to it.
Filter it with Guided Filter and Gaussian Filter.
Look at the results.
Draw a 1D Line which crosses it and look at the Original, Noisy, Filters Guided, Filtered Gaussian.

Here are the results:

Noisy Image

Gaussian Filtered Image

Guided Filtered Image

As one could see, the Gaussian Filter smooth the image across high contrast edges while the Guided Filter keep the edges.

The code is located at GitHub - DSP StackExchange Q29041.


You can't really have a covariance of a matrix.  What you can have, is a covariance matrix of a set of vectors.  So, if you think of the rows of your matrix A as two vectors in 3D: [2 3 4] and [5 5 6], then the covariance matrix of this set of two vectors is C = A' * A (A transpose times A).  Note that if you shuffle the rows of A in a different order, C does not change.

Now that you have C, you understood it correctly.  You can also do interesting things with the eigenvalues and eigenvectors of C.

Oops... Totally forgot.  Before you do anything, you have to compute the mean vector (column-wise mean of A) and subtract it from every row.


Getting started with wavelet-based real-time anomaly detection for temporal signals using Python

I suggest that you look at the Wikipedia article:

https://en.wikipedia.org/wiki/Change_detection

and first focus on the statistical aspects of the problem.  The references cited reflect most of the foundational material of the topic.  The article is of high quality. You should also think about the systems aspect of your problem like, the architecture of the system in the sense of issues like having a set of essentially special purpose detectors, or something that allows some interactions.  Your problem statement says that only a single anomaly is present at a time, and this should be reflected in the organization of system.  Questions such as Do you need automatic gain control (AGC) or do you need a constant false alarm (CFAR) behavior do have  indirect effects.

The major departure from the usual change detection problem, is $a(t)$. Typically, it is in $a(t)$ that one wishes to find changes.  I don't think that this is necessarily a problem but $a(t)$ will propagate through your detectors.

Change detectors work on changes from something typical and $a(t)$ seems like a good place to look for "typical". This is probably where your likelihood ratio test or score metric is going to reside.

I'm guessing that your essential problem is that $a(t)$ is from a remote sensor, and your anomalies are associated with a transmission channel.  If so, transmission glitches are not typically subtle low SNR events that necessitate near optimal techniques needed to find them. They are obvious and are often managed by simple heuristics. Much that follows will be irrelevant if this is the case.

Finally, Wavelets. The answer is, you are going to determine this for yourself, by trial and error and thinking through what makes sense. One heuristic is that time-bandwith matching has a small loss relative to something like a matched filter.

The choice of basis functions is specific to your application.  People use Wavelets. People use STFT, People use Wigner.  People use Cyclostationarity. People use running histograms.

It all depends on how parsimonious a basis set represents your signals of interest, but you can't neglect $a(t)$. Those basis function have a direct bearing on characterizing "typical", which is why I strongly recommend you think about the statistical organization of your system, and then work on selecting your basis functions.


How to measure distance between two 3d objects using cameras

Since nobody has answered this for you yet, I'll give it a shot.

You have to understand how to do ray tracing math, which means vectors.  You also have to be able to calibrate your camera so you know for each pixel, what is the direction of the the ray coming from the focal point.

It makes a big difference if the objects are in motion or are still.  If they are in motion you will need to deal with blurry edges.

Since the screen is two dimensional, you will only be able to make two dimensional direction measurements.  However, with something like a ball with a known diameter, you can calculate how far it is from the focal point using trigonometry.

Having two cameras on the scene can give you much more accurate depth measurements.  Again, this requires knowing the rays and doing trigonometry.

If your surface is flat, you can define a mapping between your pixel locations and a coordinate system on the surface.

Many years ago, I wrote a program using a webcam, where when I identified a sheet of notebook paper on the floor by clicking the four corners with my mouse, the program calculated the location of the camera relative to the paper within a tenth of an inch.

Once you have the objects located in 3D space, finding the distance between them is merely an application of the distance formula, aka Pythagorean Theorem.

Hope this helps.

Ced

=============================================

Followup (too long for an add on comment above):

@M.M., I've never used OpenCV so I searched on "lens distortion matrix openCV" and found this:

Camera calibration With OpenCV

This approach is different than what I did, though somewhat similar.

Suppose $(m,n)$ is the pixel location and $(x,y,1)$ is the direction vector where the z-axis is straight out from the lens.  I used calibration points, and a standard best fit approach to find the coefficients of:

$$ x = A_x m^2 + B_x n^2 + C_x m n + D_x m + E_x n + F_x $$
$$ y = A_y m^2 + B_y n^2 + C_y m n + D_y m + E_y n + F_y $$

This required finding the inverse of a 6x6 matrix to solve.  I wrote the program in VisualBasic3 (I think) on an i386, to give you an idea of when.

Ced

=============================================

More detail:

Spheres are easy because they look like a solid circle from any angle.

Find the center of the circle in pixel coordinates.  A centroid calculation is a brute force way to do this, but there are more efficient ways.

Translate the pixel coordinates of the center to a ray emanating from your lens.

Then process several points around the edge.

For each point:

1) Translate to a ray

2) Calculate the angle to the center ray.

3) Figure the distance to the center of the sphere from the lens

Average the distances (they should ideally be identical) to get a good read.

Calculate the location in 3D coordinates of the center of the sphere from the center ray direction and the calculated distance.

Calculate the distances between the spheres using the distance formula.

Angle calculation formula for two vectors $\vec a$ and $\vec b$

$$ \theta = cos^{-1} \left( \frac{ \vec a \cdot \vec b }{ \| \vec a \| \| \vec b \|  } \right) $$

Distance to sphere where $r$ is the radius

$$ d = \frac{ r }{ \sin \theta } $$

The center of the sphere:

$$ \vec c_k  = d \frac{ \vec a_k }{ \| \vec a_k \| } = ( x_k, y_k, z_k ) $$

The distance between spheres:

$$ d_{k,l} = \sqrt{ ( x_k - x_l )^2 + ( y_k - y_l )^2 + ( y_k - y_l )^2 } $$

Your surface will require a different approach, but since you haven't provided any details, I can't help you there.

Ced


Octave is a free mostly MATLAB compatible application. You shouldn’t rule it out. Also, your signals would look more alike if you removed the low frequency trend obvious in both plots. You may have a valid ambiguity function although I think you may be calculating an ambiguity on signal plus noise not the proper calculation of just signal.

There is more than one ambiguity function, this is the wide band single transmitter single receiver version. Most people think Doppler is just a frequency shift but it is actually a time dilation, although the frequency shift model is very often sufficient, particularly when looking at steady lines like molecular emission lines.

In acoustics the motion of the propagation media can also cause Doppler.

This code should be ok as long as the interpolation is accurate, which assumes that the signal has sample more than about 5 points per cycle.

 % wide band ambiguity function
    clc
    clear all
    close all
    %
    tb=-1;
    tend=1;
    dialation=.8;
    %
    %signal=[zeros(1,100) randn(1,200) zeros(1,100)];
    %signal=[zeros(1,100) exp(1j*2*pi*(4*[0:199]/200)) zeros(1,100)];
    signal=[zeros(1,100) exp(1j*2*pi*(4*[0:199]+.01*[0:199].^2)/200) zeros(1,100)];
    conj_rev_signal=conj(signal(end:-1:1));
    time=linspace(tb,tend,length(signal));
    delay_time=linspace(2*tb,2*tend,2*length(signal)-1);
    tau=linspace(dialation,1/dialation,512);
    out=zeros(length(tau),2*length(signal)-1);
    for i=1:length(tau)
    out(i,:)=conv(signal,interp1(time,conj_rev_signal,time*tau(i),'spline'));
    end
    figure(1)
    imagesc(delay_time,tau,abs(out))


will sound waves going to point A be affected by that wind

Yes. There are two major effects of wind on sound.

Noise: Wind is moving air and moving air is sound. It's mostly DC air movement so it's not audible per se, but it when the airflow hits objects it often creates wide band turbulences which can be quite loud. Particularly problematic is wind hitting the receivers i.e. a microphone or your ears. It can generate a lot of pick up noise. Wind screens are one way to deal with this.
Waveform shaping: Sound waves are moving with resect the fluid media (i.e. air). If the media moves as well, than the overall movement of the sound wave will be the sum of both motions and this deforms the wave shape. If the wind is relative constant you still get roughly spherical waves, but they are not concentric anymore, i.e. the center is moving with the wind. So your glass shatter would be significantly louder downwind then upwind.


modified Bessel function for Speech enhancement Using MMSE in matlab

The Bessel function just turns up as a result of the MMSE estimation problem given the assumed statistical model (cf. Appendix A of the paper). This modified Bessel function is implemented in Matlab as besseli.


What's the difference between nonstationary and time varying?

My understanding :
Non-Stationarity is mentioned in statistical sense where channel coefficients are non-stationary random processes. You can look for Wide-Sense Stationary Processes and Strict-Sense Stationary Processes and their definition and properties.
Time-Varying is with regards to time-domain channel impulse response characteristics. Time-variance indicates that depending upon when any specific input is provided to the channel, the output will change because the channel is varying with time.
A particular channel model can be time-varying but wide-sense stationary. This would mean, depending upon when any specific input is provided to the channel, the output changes but the statistical properties of the output like correlation function, mean etc will follow wide-sense stationarity conditions.


Is a negative Es/No for an SCPC RF signal possible, and what would such a value signify?

Assuming the number is reported in dB, a negative value of $E_s/N_0$ means is that the energy per symbol $E_s$ is less than the noise spectral density $N_0$.

Even then, a negative $E_s/N_0$ seems hard to obtain if the signal's power is larger than the noise power (as measured on the spectrum analyzer). Assuming that the maximum symbol rate is used ($E_s=2B$, where $B$ is the bandwidth), and that the noise power $\sigma^2=N_0B$, then $$\frac{E_s}{N_0} = \frac{P}{2\sigma^2},$$ where $P$ is the signal power.

So, under these assumptions, a negative $E_s/N_0$ (in dB) means that the signal's power is less than twice the noise power.


It looks like you could use a low pass filter or another noise reduction filter. Then maybe an edge detection filter.


Build discrete IIR from a transfer function which is not in terms of $s$

There exist genetic and evolutionary algorithms that attempt to guess and optimize a set of pole/zero placements that approximate an arbitrary frequency response.

There is also a algorithm called Frequency Domain Least Squares : http://www2.units.it/ramponi/teaching/DSP/materiale/ES_5_2.m

The poles and zeros found can then be used to create an IIR recurance function.


No-lag causal bandpass filter kernel formula

What you need is not a "no-lag" filter, as such a thing doesn't exist if the filter is still causal (how should it filter, without knowledge of the future and not being aloud to wait).

What you need is a filter with a constant group delay. Any linear phase filter has that. Any time-symmetric filter has linear phase. And by the simple method of shifting a time-symmetric filter by half its length, you make it causal. And have half its length in delay, at any frequency.

So, any filter kernel that is symmetric works. There's far too many bandpass filters that have linear phase, and which one you'd use completely depends on your application. But it's usually easier to find a tool to design a linear-phase filter than to find one that gives a different kind of filter, so you should be fine.


I've had this post up with no takers so I will just put what I have found so far.

The algorithm that has resulted in the least noise for my application that I have tried is:

Let $X =[x1,x2,....,xn]$ be a sample of a function over discrete time intervals where x is a real number and represents the output of a particular pixel at $\delta t*i$ where $0<i<n$ .

Find a linear function $f$ that best approximates the data in the form (I tried polynomials of varying degrees but a linear function performed the best):

$f(x) = m *x + b$

Evaluate the function at ${t }= {{n*\delta t} \over{2}}$

So $f({{n*\delta t} \over{2}}) \approx Median $

My motivation was to find the median but it may be a somewhat bold assertion.  I wanted the median value of each pixel but the chances of picking a pixel due to noise was too great.  Plus because of the physical constraints I know that the data is either decreasing or increasing over the intervals with respect to the phase of data collection.

I do this for both hot and cold images to form the image using the notation above $I=\bar{H} - \bar{C}$.


Seeking segmentation methods for this light microscope image

People were approaching similar problems with methods deriving from Watersheds thresholding. Here, is a sample of Nuclei Seperation, and here is an OpenCV sample. You might try it out.


How to chose a kernel function for a FIR low pass filter in audio processing?

Yes the sinc kernel $$h[n] = \frac{\sin(\omega_c n)}{\pi n}$$ is the theoretical result for the ideal lowpass brickwall filter whose gain is one and cutoff frequency is $\omega_c$ radian per samples.

Your filter impulse response will be samples of that function. Truncate it symmetrically to some finite length, and then apply a window on those samples. Windowing will help you achieve a much smoother response.


IIR low pass filter: floating point and fixed point

setting aside the DF1 vs. DF2 issue (and we know there are several other forms besides those two, all with advantages and disadvantages), there is an answer, sans phase response, to the OP's original question here.


How to convert a filtered time-domain signal into frequency one?

Taking the FFT will give you a frequency representation of the data:

y=fft(sig3.data);


Carrier Tracking Loop: Tracking noise vs Phase Noise Trade-off

The question is specific to optimizing loop bandwidth for the decision directed carrier tracking loop with one sample per symbol. In other posts such as this one PLL for Phase Demodulation and Carrier Tracking I have detailed the consideration for not making the loop bandwidth too low due to increasing contributions from LO phase noise as well as less ability to track system dynamics (rate of Doppler change for example). I have also suggested as the loop BW gets to large we start introducing other noise sources and lose signal by tracking out the carrier such that an optimization in loop BW is involved. The OP is making the good point that it would be impossible to track out the carrier if we are only sampling 1 sample/symbol. So the question is what others factors would limit the maximum carrier tracking loop bandwidth, or would the go be to provide the maximum bandwidth that can be achieved in a loop with an update rate $R$ of a sample/symbol (i.e. < $R/2$).

Below I show more detailed simulations of this type of decision directed carrier tracking loop. Specifically we see the closed loop transfer of noise from the input to just prior to the decision as a high pass filter (as the loop tracks out the slower changes in phase as part of tracking the carrier). We note that with a narrow loop bandwidth the gain levels out close to 0 dB at the high end of the Nyquist bandwidth, while in the wide loop bandwidth the gain of the higher frequencies is +6.8 dB, which leads to noise enhancement.

The expected input noise will typically rise significantly as dominant phase noise as we approach the carrier (low frequency offsets) as is typical for local oscillator phase noise. However at some offset from the carrier the decreasing phase noise will intersect the noise floor which will then dominate from other white noise sources such as amplified thermal noise and quantization noise in fixed point implementations. Thus we see how a trade can exist with the loop bandwidth setting where the savings from further attenuation of lower frequency noise is offset by the noise enhancement from the wider loop bandwidth setting. Note that the decision directed phase detector is equally sensitive to AM and PM components on the signal, and the AM noise components of the input signal will be translated to PM components on the NCO in the loop and will thus not cancel out the AM noise but add to it as an uncorrelated noise source. Ultimately if an optimization is desired, and to confirm if this particular noise contribution will even affect the receiver,the simulation and analysis would need to be done with the expected noise profile for the specific system being analyzed, and importantly while operating at senstivity. So there is not a one answer fits all but hopefully this makes it clear how the SNR could be degraded by the carrier tracking loop if the bandwidth is set too high and what the considerations are for analyzing this.

Simulation Result with Phase and AWG Noise

Here is a complete simulation of the Decision Directed Carrier Tracking Loop with a 2nd Order PI filter (as a typical implementation) showing how the loop bandwidth can be optimized for pre-decision SNR performance, specifically showing conditions where the noise can degrade as the loop bandwidth continues to get wider, even if it is not at the point of tracking out the modulation.

Given the following Carrier Tracking Loop (Applicable to BPSK, QPSK, M-PSK and QAM modulations):

The following is the Loop Model and equation for the open loop gain:

$P$: Proportional gain constant

$I$: Integral gain constant (per sample with time in seconds)

$IT$: Integral gain consant (normalized with time units in samples)

$A$: complex magnitude of samples at discriminator input

$K_d$: Discriminator gain includes parasitic 1 sample delay in loop

$H_{LF}(z)$: PI Loop Filter

$K_{NCO}$: NCO gain, control word given in frequency of Hz (this way we can monitor frequency directly)

$$H_{LF}(z) = \frac{Pz + IT-P}{z-1}$$

$$K_{NCO} = 2\pi\frac{z}{z-1}$$

$$K_d = \frac{A^2}{z}$$

$$G_{OL}(z) = 2\pi PA^2\frac{z-\bigg(\frac{P - IT}{P}  \bigg)}{(z-1)^2}$$

With this Loop Model I created three test cases, for all cases $A=1$ and $T=1$ (normalized magnitude and time is in samples):
$$\begin{bmatrix}  \text{Loop BW} & -3 \text{ dB BW}       & \text{Settles In}  & I     & P    & G_{CL}(z=\pi)   \\
\text{---------}&\text{---------------}&\text{-------------}&\text{-------}&\text{------}&\text{----------------}\\
\text{Narrow}  & 0.011 \text{ cyc/samp} & 200 \text{ samples} & 1e-4  & 7e-3 & 1.045 \\
\text{Medium}  & 0.056 \text{ cyc/samp} & 20 \text{ samples} & 0.01  & 0.07 & 1.257 \\
\text{Wide}  & 0.202 \text{ cyc/samp} & 5 \text{ samples} & 0.072  & 0.21 & 2.206 \\
 \end{bmatrix}$$

I then created a noise sample of $2^{15}$ complex samples with the following phase noise target values for the PSD:

Freq $1=2\pi$,   Phase Noise (dBc):

1e-2 cycles/sample, -15 dBc

1e-1 cycles/sample, -45 dBc

0.5 cylces/sample, -59 dBc

To this phase noise data I added AGWN with a total power of -40 dBc. This created experimental data of complex samples with both AM and PM noise components, with the phase noise contributing AM only with increasing density toward the lower frequencies, and AWGN contributing AM and PM components equally.

The loop performance was characterized by determining the closed loop transfer function from the input of the phase rotator to the output of the phase rotator which is given by:

$$G_{CL}(z) = \frac{1}{1+G_{OL}(z)} $$

It is clear that signal component, here normalized to 1 would not be affected in magnitude by the phase rotator, so we can assess SNR from the noise directly after it has passed through the above transfer function.

The Decision Directed Phase Detector responds equally to small scale AM and PM changes (A small change in amplitude can not be distinguished from a small change in phase), so any AM components would get translated to (uncorrelated) PM noise at the output of the phase rotator as the loop tries to correct phase offsets that don't exist. Therefore the resulting total noise at the output taken at the phase rotator output would be:

$$S_o(\omega) = S_{am}(\omega) + |G_{CL}(\omega)|^2S_{\phi}(\omega) + |G_{CL}(\omega)|^2S_{am}(\omega) $$

Where $S_o(\omega)$ is the power spectral density of the noise alone. Which means all the AM noise original on the signal would still exist without modification since the loop does not correct for any AM, the PM noise would be filtered by the loop, but the AM noise would also be translated to PM by the loop and this component of translated AM noise would also be filtered by the loop.

So accounting for all that, using this experimental data, here were the results after passing the data set through the narrow, medium and wide bandwidth loop filter responses:

$$\begin{bmatrix}  \text{Loop BW} & \text{Total Integrated Noise}  \\
\text{---------}&\text{-----------------------------}\\
\text{Narrow}  & -64.6 \text{ dBc} \\
\text{Medium}  & -69.8 \text{ dBc} \\
\text{Wide}    & -67.1 \text{ dBc}\\
 \end{bmatrix}$$

I assumed the optimum loop bandwidth for minimizing noise would be around the corner frequency where the phase noise intersects with the white noise floor but haven't confirmed this, nor did I do a detail sweep beyond these three test cases. I wanted to show an example that demonstrates how wider loop bandwidths can cause the noise level to start increasing as bandwidth gets wider and for this loop it is because of noise enhancement for wider loop bandwidths. We can see how this must be the case given the two closed loop zeros at $z=1$, and that the poles move closer to $z=-1$ as loop bandwidth increases.

Also to note, the resulting overall noise levels in this test case are inconsequential for demodulating most communication waveforms, so it would make no difference to degrade the tracking noise by 3 dB so that we could operate at the wider loop bandwidth, with all its other benefits. The point here was to show that it can be optimized but also that it very much depends on the overall system and there is not a one answer for all other than to evaluate this on a candidate system with similar detail.

The above is the concluding data and remarks, but further plots are shown below with regards to the test data, loop step responses and closed loop transfer functions.

Experimental Data - Phase Noise Only (time domain)

Power Spectral Density of Phase Noise Data

Power Spectral Density of Phase Noise with AWGN Data

Closed Loop Step and Frequency Response From Input to Phase Rotator Output


Well I think I found a solution to this, mostly by trial and error.

So what I did and its giving quite good results is :
1) If we have M(No. of frames) X N(13;No of MFCC Coefficients) matrix;

ie X=X[1 2 .....N(13) ; 1 2.....M(No of frames) ]

2) Find the DTW of [X(1;1 2 3.....M)] with the DATABASE [REF(1;1 2 3...M)]
   Here we are taking K=1(Not 2...12)
   Similarly find DTW of [X(2;1 2 3...M)]  and DATABASE [REF(1;1 2 3...M)]

3) On generalizing, Find DTW of [X(k(kth MFCC coef); 1 2 3...m)];

4)So at the end we will get : K(13) values ;

5) Find the mean of these 13(K) values and let that be your Reference for      deciding;

*Disclaimer: I am not certain if this is what all the papers imply, I just found it by trial and error and it gives accurate results


Method for modifying frequency-domain signal to produce predetermined time-domain signal

Zero-padding is one domain is equivalent to Sinc (or Dirichlet) convolution in the other, which is O(n^2), so that won't save you.  For less computation, you might want to try different filter design methods (compositing existing FIR filters, IIRs, etc.) or finding a way to recompute your filter less often.


Deconvolution Using Response to an Heavy Side

If we can assume no noise (Or the SNR is very high) you can get the response by applying the inverse filter in frequency domain.

Lets say $ y [n] $ are the signal samples.
Given $ x [n] $ the samples of the ideal signal you can apply on both of them the DFT to get $ Y [k] $ and $ X [k] $.

The response is given by the Inverse DFT of the division $ \frac{Y[k]}{X[k]} $.

If there's some noise you need to regulate the result.
Easy choice would be the Wiener Filter or more specifically Wiener Deconvolution.


Decimation vs Mean in downsampling operation

What you are calling decimation is in fact downsampling: keeping one of every $k$ samples from the original signal.

Proper decimation involves a digital low-pass filter in order to eliminate all components above $\pi/k$ (or $f_s/(2k)$ in continuous time) since these components cannot be represented at the new lower sampling frequency (and will result in aliasing).

After low-pass filtering, you downsample your signal keeping one out of every $k$ samples.

The second method you propose is actually a low-pass filter followed by downsampling. But in this case, the low-pass filter you use is a moving-average filter, which is very simple, but it is not a very good low-pass filter. Gain in the pass band is very variable, and attenuation in the stop band is poor.


Chosing polyphase interpolation coefficients

At Peter K's suggestion I am submitting my solution as the answer. In the poly phase filter the 9 coefficients were all symmetric about the center but gave poor ripple performance. When arranged in the format of a poly-phase filter the 3 coefficients in the 3 branches did not add up to the same value for all three branches. By making a slight edit to the 1 coefficient at +/-3 positons from the center all the branches added up to the same value and the results are as shown.

To get better performance then that shown I need to increase the # of coefficients to N=15 (wanting to keep N a multiple of L=3 and odd).

Brian


Why Fourier series if Fourier transform can be calculated for both periodic and aperiodic?

You can't take the Fourier transform of a periodic signal, the integral diverges for all multiples of the period.

This can be handled by the theory of distributions, but the Fourier series is a better fit.


Small Scale Fading time varying or position varying

Small scale fading results from a superposition of a large number of propagation paths between TX and RX, which can be constructive or desctructive depending on their phase offset. As the phase offset varies rapidly with the propagation distance, the phase offsets between different paths are likely to vary rapidly too, leading to a pseudo-random behaviour of the received amplitude. We refer to this as small scale fading.

Technically, this behavior requires some form of change to the propagation conditions. So again, technically, if you did an experiment in a closed chamber with noone inside and no movements, there would not be any fading but the amplitude would remain constant.

However, in reality, there is always some movement. Most of the propagation paths are due to reflections and reflecting objects can move on their own (think of wind, people, cars, etc.). So in practice, even if TX and RX are stationary, you should expect small-scale fading to be present. A way to reduce it is using directive antennas that put stronger emphasis on the line of sight component, if that is an option.


Your derivation doesn't give the correct result because the limit $\lim_{t\to\infty}e^{-j\omega t}$ doesn't exist. You seem to have confused that limit with $\lim_{t\to\infty}e^{-\omega t}$, which does exist and which equals zero for $\omega>0$.
There are many ways to derive the Fourier transform of the unit step function $u(t)$ (see e.g., this answer), but I'll show you one which doesn't involve more than just engineering math. At least if you accept some intuitively pleasing facts such as that the derivative of the unit step equals the Dirac delta impulse:
$$u'(t)=\delta(t)\tag{1}$$
It's important to note that $(1)$ holds for any function with a discontinuity at $t=0$ that is otherwise constant, satisfying $u(0^+)-u(0^-)=1$. Hence, $(1)$ only defines $u(t)$ up to a constant. Taking the Fourier transform of $(1)$ gives
$$j\omega U(\omega)=1\tag{2}$$
It would be tempting to conclude from $(2)$ that $U(\omega)=1/j\omega$ holds, however, this is not correct. The reason is the additive constant that got lost due the derivative in $(1)$.
Hence, the correct conclusion from $(1)$ is
$$U(\omega)=c\delta(\omega)+\frac{1}{j\omega}\tag{3}$$
where $c$ is the constant that can be added to any function satisfying $(1)$ without violating Eq. $(1)$. We just have to find the specific value of $c$ for which $(3)$ holds, with $U(\omega)$ being the Fourier transform of the unit step.
It's easy to see for which function $c=0$ holds. Since $1/j\omega$ is purely imaginary, we know that the corresponding time domain function is odd, it has a jump of height $1$ at $t=0$, and it is otherwise constant. The only function satisfying all of these requirements is a scaled signum function, i.e.,
$$\frac12\textrm{sgn}(t)\Longleftrightarrow \frac{1}{j\omega}\tag{4}$$
Since the unit step function satisfies
$$u(t)=\frac12+\frac12\textrm{sgn}(t)\tag{5}$$
we finally arrive at the Fourier transform of $u(t)$ by combining $(4)$ and $(5)$:
$$\bbox[#f8f1ea, 0.6em, border: 0.15em solid #fd8105]{
U(\omega)=\pi\delta(\omega)+\frac{1}{j\omega}}\tag{6}$$


Engine for Google / Browsers to detect manipulated images?

Overview:
Detecting fake images is part of the research field of media forensics. This is an ongoing research niche.
There already exists a fairly mature set of computational tools that can identify various types of image manipulations. They can also indicate the region in the image that has been altered. Concurrently, the improvements in photo-editing tools, like Adobe Photoshop, had made the detection of manipulated images more challenging.

Would it be possible for Google and web browsers to scan images and photos, to detect fakes?

Yes, it is indeed possible. Integration into a popular web browser is probably pending, but the technology already exists. You can try it for yourself here
There is also an open source repository for those interested to contributed or modify here

For example, a web browser could check for abnormal color gradients and signals that indicate if a photo is fake?

A popular example of fake photo is the missile launch by Iran:

The actual lunch had only 3 missiles and 1 missile-launcher failed:

The manipulated image used copy-move forgery from the original image to synthesize the fake 4th missile. A clone-tool detection algorithm can detect this:

Google could fingerprint photos into ASCII code, then compare to elements copied from older photos. If that was done in conjuction with people appearing in news articles, accidents, or porn, forgetting all the animal and technology photos, it wouldn't be a programming impossibility?

There are examples of easy to detect fakes when the manipulator has used existing images. Typically this happens when the manipulator uses stock photography. So the answer is Yes, you are right! However, this is a small subset of fake images and more skilled manipulator have better ways of covering their tracks. Additionally, it does not work when the fake image is completely new.

I was imagining a google images page where all the fake photos were shown.

This question raises the point of provenance of images. The image in raw format from the sensor (camera) is the original. Every operation subsequently done to the image is by definition a manipulation. This includes saving the image in lossy format like jpeg. (Skilled image manipulators will cover their tracks by doing a double jpeg compression!) However, some manipulations are malicious, whereas others are not. There can indeed be a Google page with all fake photos (ie. photos with lower provenance than the 'original' photo). This might even become a requirement in photo-journalism in the near future!
There were multiple interesting sub-questions raised here. I hope I have at least partially answered most of them.


Time Frequency Signal Analysis in Mathematica (MMA)

If you are using an older version of Mathematica (pre v.8) and are interested in wavelets - yes, you need an add-on to perform wavelet analysis. More about it here. If you are using v.8 or above then everything wavelet-related is built-in. If you are interested in Fourier analysis - it is built in (since v.1). To showcase the capabilities of the time-frequency analysis functionality:

Suppose we have the signal $2 \exp \left(-\frac{t^2}{10}\right) \cos (5 t)$

We can visualise its Gabor transform

Wigner transform (notice the artifacts)

and the (somewhat) corrected Gabor-Wigner transform

A second test signal

and its spectrogram calculated using partitions of length 256, offset 1 and BlackmanHarrisWindow

its periodogram (with the same options as above)

Now we move onto wavelets
Scalogram after a continuous wavelet transform with a Morlet wavelet

and the respective scalogram after performing a discrete wavelet transform

Just because I find it pretty - a scalogram of a simulated noise in 3D


I've been reading about this and there are multiple ways to do it, using different size N.  My Matlab is rusty, so here they are in Python (N is length of input signal x, k is arange(N) = $[0, 1, 2, ..., N-1]$):

Type 2 DCT using 4N FFT and no shifts

Signal [a, b, c, d] becomes

[0, a, 0, b, 0,  c,  0,  d,  0,  d,  0,  c, 0, b, 0, a].

Then take the FFT to get the spectrum

[A, B, C, D, 0, -D, -C, -B, -A, -B, -C, -D, 0, D, C, B]

then throw away everything but the first [A, B, C, D], and you're done:

u = zeros(4 * N)
u[1:2*N:2] = x
u[2*N+1::2] = x[::-1]

U = fft(u)[:N]
return U.real

Type 2 DCT using 2N FFT mirrored (Makhoul)

[a, b, c, d] becomes [a, b, c, d, d, c, b, a].  Take the FFT of that to get
[A, B, C, D, 0, D*, C*, B*], then throw away everything but [A, B, C, D] and multiply it by $e^{-j \pi k \over 2 N}$ (half-sample shift) to get the DCT:

y = empty(2*N)
y[:N] = x
y[N:] = x[::-1]

Y = fft(y)[:N]

Y *= exp(-1j*pi*k/(2*N))
return Y.real

Type 2 DCT using 2N FFT padded (Makhoul)

[a, b, c, d] becomes [a, b, c, d, 0, 0, 0, 0].  Take the FFT of that to get
[A, B, C, D, E, D*, C*, B*], then throw away everything but [A, B, C, D] and multiply it by $2e^{-j \pi k \over 2 N}$ to get the DCT:

y = zeros(2*N)
y[:N] = x

Y = fft(y)[:N]

Y *= 2 * exp(-1j*pi*k/(2*N))
return Y.real

Type 2 DCT using N FFT (Makhoul)

Signal [a, b, c, d, e, f] becomes [a, c, e, f, d, b], then take the FFT to get [A, B, C, D, C*, B*].  Multiply by $2 e^{-j \pi k \over 2N}$ and then take the real part:

v = empty_like(x)
v[:(N-1)//2+1] = x[::2]

if N % 2: # odd length
    v[(N-1)//2+1:] = x[-2::-2]
else: # even length
    v[(N-1)//2+1:] = x[::-2]

V = fft(v)

V *= 2 * exp(-1j*pi*k/(2*N))
return V.real

On my machine, these are all roughly the same speed, since generating exp(-1j*pi*k/(2*N)) takes longer than the FFT.  :D

In [99]: timeit dct2_4nfft(a)
10 loops, best of 3: 23.6 ms per loop

In [100]: timeit dct2_2nfft_1(a)
10 loops, best of 3: 20.1 ms per loop

In [101]: timeit dct2_2nfft_2(a)
10 loops, best of 3: 20.8 ms per loop

In [102]: timeit dct2_nfft(a)
100 loops, best of 3: 16.4 ms per loop

In [103]: timeit scipy.fftpack.dct(a, 2)
100 loops, best of 3: 3 ms per loop


Does the DTFT of $\frac{u[n-1]}{n}$ exist?

Note that the sequence

$$x[n]=\frac{u[n-1]}{n}\tag{1}$$

is in $\ell^2(\mathbb Z)$ because

$$\sum_{n\in\mathbb{Z}}|x[n]|^2=\sum_{n=1}^{\infty}\frac{1}{n^2}=\frac{\pi^2}{6}<\infty\tag{2}$$

but it is not in $\ell^1(\mathbb Z)$ since

$$\sum_{n\in\mathbb{Z}}|x[n]|=\sum_{n=1}^{\infty}\frac{1}{n}=\infty\tag{3}$$

We know that for sequences in $\ell^1(\mathbb{Z})$ the DTFT $X(e^{j\omega})$ exists and it is a continuous function of $\omega$. However, it is not necessary for a sequence to be in $\ell^1(\mathbb{Z})$ for the DTFT to exist. We can extend the definition of the DTFT to sequences in $\ell^2(\mathbb{Z})$, which implies that the DTFT converges in a mean square sense (and not uniformly), and that $X(e^{j\omega})$ may have discontinuities.

In that extended sense, the DTFT of $(1)$ exists, just as it exists for the ideal discrete-time differentiator or the discrete-time Hilbert transformer, the impulse responses of which also only decay as $1/n$.

The DTFT of $x[n]$ is

$$X(e^{j\omega})=\sum_{n=1}^{\infty}\frac{e^{-jn\omega}}{n}\tag{4}$$

The series $(4)$ is the complex Mercator series

$$\sum_{n=1}^{\infty}\frac{z^{n}}{n}=-\log(1-z),\qquad |z|\le 1,\quad z\neq 1\tag{5}$$

with $z=e^{-j\omega}$. It is the Taylor series for $-\log(1-z)$ and it converges for all $|z|\le 1$, $z\neq 1$. From $(5)$ we can immediately write down the $\mathcal{Z}$-transform of $x[n]$:

$$X(z)=\sum_{n=1}^{\infty}\frac{z^{-n}}{n}=-\log(1-z^{-1}),\qquad |z|\ge 1,\quad z\neq 1\tag{6}$$

Since $(6)$ also converges for $|z|=1$ except for $z=1$, we can obtain the DTFT directly from $(6)$ by setting $z=e^{j\omega}$. For $e^{j\omega}=1$, i.e., $\omega=0$, the DTFT has a singularity.

  $$X(e^{j\omega})=-\log(1-e^{-j\omega})\tag{7}$$

Equation $(7)$ agrees with WolframAlpha's result.
Also take a look at this related answer.

Let's visually check if it is plausible that the series $(4)$ really converges to $(7)$. I've computed the sum over the first $100$ terms of $(4)$, and I've plotted the real and imaginary parts of that partial sum together with the real and imaginary parts of $(7)$:

From the figure, it appears that the series $(4)$ indeed converges to $(7)$ in a mean square sense. This is of course no proof but a simple sanity check.


Understanding of Random Process/Random Variable

Let me explain it in another way. Consider you have 6 different function of time.
You only throw your dice once and regarding the outcome you choose on of six functions and the chose one is one outcome of your random process.

Now expand this concept to all possible functions of time instead of only 6 functions and instead of dice throwing consider another form of random event which could append a probability to all (infinite) possible functions. The set of all these possible functions is your ensemble.


What does it mean to deconvolve the impulse response

A linear time-invariant (LTI) system is completely characterized by its impulse response $h(t)$. The output signal $y(t)$ given an arbitrary input signal $x(t)$ is given by the convolution of this input signal with the system's impulse response:

$$y(t)=(x\star h)(t)\tag{1}$$

where $\star$ denotes convolution. In the frequency domain, Eq. (1) becomes a simple multiplication:

$$Y(f)=X(f)H(f)$$

where $Y(f)$, $X(f)$, and $H(f)$ are the Fourier transforms of the signals $y(t)$, $x(t)$, and $h(t)$, respectively.

Deconvolution now refers to the process of "undoing" the effect of an LTI system. I.e. you're given the signal $y(t)$ and you want to recover the original signal $x(t)$. E.g., you might want to remove the reverberation added to a signal by a given room. In this case you need to undo the effect of the convolution of the room's impulse response with the original signal.

Naive deconvolution can in principle be done by applying an inverse filter to the signal $y(t)$. The inverse filter ideally has a frequency response $G(f)=1/H(f)$. Note that such a filter does not necessarily exist, e.g. if $H(f)$ has zeros at certain frequencies. Another disadvantage of this naive deconvolution is the effect it has on measurement noise. In frequency regions where $H(f)$ has a small magnitude, the magnitude of $G(f)$ becomes large, and the noise is highly amplified. This is why applying an inverse filter is usually not a good solution. If the system's impulse response is known and if the noise can be modeled as (wide-sense) stationary, Wiener deconvolution can be used. In that case, the original signal is estimated using a mean square error criterion, which takes into account the impulse response as well as the noise. If the impulse response is unknown (or can't be estimated), there are methods available referred to as blind deconvolution.


Smoothing can be done in many ways, but in very basic and general terms it means that you even out a signal, by mixing its elements with their neighbors. You smear/blur the signal a bit in order to get rid of noise. For example, a very simple smoothing technique would be, to recalculate every signal element f(t) to as 0.8 of the original value, plus 0.1 of each of its neighbors:

f'(t) = 0.1*f(t-1) + 0.8*f(t) + 0.1*f(t+1)

Note how the multiplication factors, or weights, add up to one. So if the signal is fairly constant, smoothing doesn't change it much. But if the signal contained a sudden jerky change, then the contribution from its neighbors will help to clear up that noise a bit.

The weights you use in this recalculation function can be called a kernel. A one-dimensional Gaussian function or any other basic kernel should do in your case.

Nice example of one particular kind of smoothing:

Above: unsmoothed signal
Below: smoothed signal

Examples of a few kernels:


These are interesting results.  Your comment "The X axis represents the distance travelled in metres (from 3.5m away to 0.5m away from the sound source). " seems to contradict "With this explanation, '0' represents 0.5m from sound source, while '3.0' represents 3.5m from sound source."  The latter seems to be the one that fits your data better.

What is more pronounced is the fluctuations you are getting.  Taking the speed of sound at 343 meters per second:

$$ \frac{343~\text{meters}/\text{second}}{500~\text{cycles}/\text{second}} = 0.686~\text{meters/cycle} $$

With standing waves, the nodes are at half the wavelength, which corresponds closely to the size of your fluctuations in the graph.  So as your robot is moving through the room, it looks like you are measuring the standing wave nodes.