Is there a way to determine the quality of an impulse-response relationship?

I don't yet see the connection of your input with "impulses" in that if I understand correctly you are looking at the settled state from prior changes rather than dynamic responses. What you are asking for I believe is a statistical correlation between "Stream Flow" and "Stream Height"; given a stream flow (in the steady state), what is the "settled" stream height that would be sustained by that flow. I would argue that to be a correlation problem to the degree you are looking for a linear relationship between the two. We have other posts that detail how to compute the correlation coefficient, if you need further info on that you could look here:

Noise detection


Unexpected imaginary part in the fft of a zero-padded cosine

A signal that's discrete in frequency is periodic in time. By using an FFT you inherently make this assumption and the FFT calculates the spectrum of a signal that's infinitely repeated in time.
The FFT is ONLY real valued if the time domain signal has even symmetry, i.e. $x(-t) = x(t)$. For a discrete periodic signal of period N this becomes $x[n] = x[N-n]$. That holds true for a pure cosine (at least for some frequencies), but not for the zero padded version.


Given an experimental FFT plot, how to find the best match from a set of generated FFT plots?

Regarding the problem with high correlation values, this happens when both signals are not zero-mean.

Assume two signals $x[k]$ and $y[k]$ of length $N$ each. Each signal can be decomposed into the sum of a constant and a mean-free signal as follows:

$$x[k] = \overline{x} + x_0[k]$$

with $\overline{x}$ being the mean

$$\overline{x} = \frac{1}{N} \sum_{k=0}^{N-1}x[k]$$

and $x_0[k] = x[k] - \overline{x}$. Do the same to get

$$y[k] = \overline{y} + y_0[k]$$

Correlating both sequences at zero lag yields (ignoring scaling here, but this does not change the issue)

$$c = \sum_{k=0}^{N-1}x[k] \cdot y[k] = \sum_{k=0}^{N-1}\overline{x} \cdot \overline{y} + \overline{x} \cdot y_0[k] + \overline{y} \cdot x_0[k] + x_0[k] \cdot y_0[k]$$

Since $x_0[k]$ and $y_0[k]$ are zero-mean, the cross-terms cancel out so the final result is

$$c = N \cdot \overline{x} \cdot \overline{y} + \sum_{k=0}^{N}x_0[k] \cdot y_0[k]$$

So if both correlated signals are not mean-free, the correlation will always contain the product of both means times the amount of samples. If you don't want this and the mean is not relevant, you should subtract it prior to correlating.

Regarding the Matlab question. You can pack all reference signals next to each other in a matrix and then correlate with all at once at zero-lag. If you have $M$ reference signals of length $N$ each, you could do (assuming 'signal' as a column vector and the reference signals saved in columns of 'ref' as well)

c = sum(repmat(signal, 1, M) .* ref)

to obtain a vector row vector c containing the correlations of your experimental signal with all the reference signals ref. No need to calculate the whole correlation or convolutions which is indeed a huge overhead. If you do not want to use repmat, you can also have a look at bsxfun which can achieve the same result and is generally interesting when you want to do element-wise operations with a vector on each row or column of a matrix.


It's called anti-aliasing and happens when you low-pass filter your signal (image) and then decimate it (scale it down). Hence, the way to get this is take a large picture with sharp edges, and image->scale image, using anything but "None" as method. Alternatively, blur filter (==low pass) or selective blur will also do.


Is there a transformation filter to decode light signal through Glass Bricks?

In theory this is possible, in practice probably very hard, but I wouldn't say impossible. See a recent article by the Silberberg group for example, that demonstrates imaging through 'thin' turbid layers using scattered incoherent light. Also, search for works pioneered by A. P. Mosk on that matter (for example this one and this one). The main challenges are algorithm speed, or, by the time you'll find the correct transmission matrix, it'll change due to thermal fluctuations.


Sinusoidals are one of the easier signals to analyze given how you can manipulate the expression by algebraic tricks.

for this function particularly, I would convert sine to its euler form first:

Then use the Fourier Transform definition:


Compressive Sensing: What Class of Signals Are Exactly Model Sparse?

In the case of general Compressed Sensing, we assume that we have no prior information about the location of the sparse support of a given signal. If we knew the support exactly, that is, if an oracle told us where the non-zero coefficients are located, we would only need $K$ measurements to recover our signal accurately. Since we do not generally have such a genie, we require more measurements to both determine both the support and the signal values on that support, hence $M$ on the order of $O(K \operatorname{log}N)$.

The authors in the Model-Based Compressed Sensing technique make the observation that if an a priori assumption can be made about the location of signal support, fewer measurements are required (since we know more about the support, statistically, we need to deduce less information from our measurements).

When a signal obeys a tree-structured, or hierarchical, model for its support, we have an easy way of predicting the support. This knowledge allows us to lower the number of required measurements for accurate reconstruction. This kind of model is very common since many classes of interesting natural signals can be represented compactly using Wavelets. The multiresolution nature of wavelet analysis groups significant (i.e. non-zero) coefficients into a tree-structure: significant coefficients at one resolution scale are highly correlated with spatially co-located coefficients in neighboring resolution scales.

(from: http://www.dominikdeak.com/coding/wavelet-compression/figure_17.png)

In the example you give about a mixture of Gaussians, it is not immediately clear that Model-Based CS can be used, since, as you say, the parameters (the peaks) of these Gaussians are randomly distributed within the signal. This kind of signal can be sparsely-represented (pretty much ideally) with some redundant dictionary of Gaussians of different parameters.
The sparse support when represented within this dictionary would not necessarily obey
any kind of model, the measurements alone would be needed to determine the support (Gaussian parameters).

On the topic of Compressed Sensing recovery of Gaussian Mixture Models, I ran a quick google search and came up with the following on arXiv:

Yu and Sapiro, "Statistical Compressed Sensing of Gaussian Mixture Models."


Does "keying on" a sine wave at a zero-crossing reduce its bandwidth?

It is instructive to have a look at the expression for the spectrum of a truncated sine wave to see how the choice of the window length affects the decay of the spectrum. Let's define a windowed sinusoid $s(t)$:

$$s(t)=\sin(2\pi f_0t)\cdot\text{rect}(t,T)=\begin{cases}\sin(2\pi f_0t),&-T/2<t<T/2\\0,& \text{otherwise}\end{cases}\tag{1}$$

The Fourier transform of $s(t)$ is the convolution of the transforms of the sinusoid and the rectangular window:

$$\begin{align}S(f)&=\frac{1}{2j}\left[\delta(f-f_0)-\delta(f+f_0)\right]\star T\text{sinc}(Tf)\\&=\frac{T}{2j}\left[\text{sinc}(T(f-f_0))-\text{sinc}(T(f+f_0))\right]\tag{2}\end{align}$$

With $\text{sinc}(x)=\sin(\pi x)/(\pi x)$ we get from $(2)$

$$\begin{align}S(t)&=\frac{T}{2 j}\left[\frac{\sin(\pi T(f-f_0))}{\pi T(f-f_0)}-\frac{\sin(\pi T(f+f_0))}{\pi T(f+f_0)}\right]\\&=\frac{1}{2\pi j}\frac{(f+f_0)\sin(\pi T(f-f_0))-(f-f_0)\sin(\pi T(f+f_0))}{f^2-f_0^2}\\=&\frac{1}{2\pi j(f^2-f_0^2)}\Big\{f\big[\sin(\pi T(f-f_0))-\sin(\pi T(f+f_0))\big]+f_0\big[\sin(\pi T(f-f_0))+\sin(\pi T(f+f_0))\big]\Big\}\tag{3}\end{align}$$

The spectrum $(3)$ decays as $\propto 1/f$ because we have $f^2$ in the denominator and $f$ in the first term of the numerator. If the window length $T$ is chosen such that an integer number of periods of the sinusoid are inside the window, then the limits of the window will be exactly at zero crossings of the sinusoid. This condition is equivalent to $Tf_0=n,\;n\in\mathbb{Z}$. With this condition the first term in the numerator of $(3)$ vanishes because the arguments of the two sinusoids differ by $2\pi n$, and, consequently, both sinusoids are identical. As a consequence, in that case the $f$ in the numerator disappears and the spectrum $S(f)$ will decay as $\propto 1/f^2$ instead of $\propto 1/f$.

In sum, if the windowing does not cause any discontinuities (i.e., the window limits occur at zero crossings of the sinusoid), the spectrum of the windowed sinusoid will decay as $\propto 1/f^2$. This is generally the case for continuous functions. Such functions can of course decay even faster, but only if they also have continuous derivatives, which is not the case for a windowed sinusoid. If the windowed sinusoid is discontinuous, i.e., the window limits do not occur at zero crossings, then the corresponding spectrum only decays as $\propto 1/f$.


This is very simple. You need an inter-disciplinary domain where all (and more) of your interests are satisfied, while being very publishable, employable and fun. Easy!! Jump into the bleeding edge of Speech Recognition.

There is TONS of current research is Automatic Speech Recognition (ASR). The domain is very multi-disciplinary (applied math, linguistics, computer science, optimization, electrical engineering) with applications to just about every part of life where a speech interface can work (healthcare, mobile assistants, in-car, entertainment/gaming, education, etc)

For example, state of the art systems like Baidu's Deep Speech 2 use convolutional neural networks (CNNs) to process speech spectrograms (DSP!). CNNs have revolutionized image processing since 2012. The CTC objective that is optimized during training is a very interesting loss function for supervised training of sequences.

Finally, have a look at these 2016 Stanford University student projects, notably those on speech recognition. This should give you some ideas of what scope or size of project makes sense.

Conditioning WaveNet on Learned Formant Characterizations for Speech
Audio Enhancement
Music-Speech Discrimination
Semantic Analysis of Political Speeches
Language Recognition and using Speech Recordings
Reinforcement Learning for Feature Selection in Affective Speech Classification

Good luck!!


Can't understand multi intensity valued mask

...masks are more refined images only emphasizing required portions and setting remaining to 0.

A mask is usually a binary bitmap that uses 1 at some location $(x,y)$ to indicate that it "selects" or "uses" or "includes" that pixel in the region of interest it defines.
The values that you see in your mask are not so much related to the mask but to its labeling.

...in the given mask, how come there are different intensity values for nuclei?

Each intensity value "tags" a different region. In other words, if your mask is in some $M$, then by evaluating $M == u$ where $u$ is a positive integer (say for example 5), then you would "select" all (connected) pixels that belong to the region of interest that has been labeled as "5".

Is it due to the fact that some nuclei are strongly represented and some are weakly represented?

No.
Also, "no" to #2 and "yes but not exactly" to #3: To get the behaviour described in #3 you would have to first "select" the mask (which casts it to a boolean) and then use it in logic expressions.
So, to get your bounding boxes all you have to do is iterate through the labels of the mask, extract it (with something like "where mask==label, set to 1 else 0" and find the min, max of the $x,y$ components of the pixel location.
Hope this helps.


Matlab plot of QPSK system doesn't agree perfectly with theoretical BER curves

The simple explanation is that there is an error in your simulation. Here's one that works in MATLAB:

% number of symbols in simulation
Nsyms = 1e6;
% energy per symbol
Es = 1;
% energy per bit (2 bits/symbol for QPSK)
Eb = Es / 2;
% Eb/No values to simulate at, in dB
EbNo_dB = linspace(0, 10, 11);

% Eb/No values in linear scale
EbNo_lin = 10.^(EbNo_dB / 10);
% keep track of bit errors for each Eb/No point
bit_err = zeros(size(EbNo_lin));
for i=1:length(EbNo_lin)
    % generate source symbols
    syms = (1 - 2 * (randn(Nsyms,1) > 0)) + j * (1 - 2 * (randn(Nsyms, 1) > 0));
    % add noise
    syms_noisy = sqrt(Es/2) * syms + sqrt(Eb/(2*EbNo_lin(i))) * (randn(size(syms)) + j * randn(size(syms)));
    % recover symbols from each component (real and imaginary)
    syms_rec_r = sign(real(syms_noisy));
    syms_rec_i = sign(imag(syms_noisy));
    % count bit errors
    bit_err(i) = sum((syms_rec_r ~= real(syms)) + (syms_rec_i ~= imag(syms)));
end
% convert to bit error rate
bit_err = bit_err / (2 * Nsyms);

% calculate theoretical bit error rate, functionally equivalent to:
% bit_err_theo = qfunc(sqrt(2*EbNo_lin));
bit_err_theo = 0.5*erfc(sqrt(2*EbNo_lin)/sqrt(2));
figure;
semilogy(EbNo_dB, bit_err, 'bx', EbNo_dB, bit_err_theo, 'r', 'MarkerSize', 10, 'LineWidth', 2);
xlabel('Eb/No (dB)');
ylabel('Bit error rate');
title('QPSK bit error rate');
legend('Simulation','Theory');
grid on;

Note that the theoretical expression for bit error rate for BPSK/QPSK modulation is:

$$
P_b = Q\left(\sqrt{\frac{2E_b}{N_0}}\right)
$$

keeping in mind that $E_b$ is the energy per information bit. The somewhat subtle distinction between $E_b$ and $E_s$, the energy per symbol, is something that often trips up people new to the subject. This difference also explains why QPSK and BPSK have the same bit error rate when expressed as a function of $\frac{E_b}{N_0}$; you don't get any bit-error performance benefit by moving to QPSK, although you can achieve a given bit rate with less occupied bandwidth.


Real signal from I,Q, sample rate and center frequency

I gather you are trying to reconstruction a real (non-imaginary) version of your I & Q signal. The following scheme should work:

Up sample your IQ from 10 MHz to 24 MHz. This involves zero padding and applying a low-pass filter to remove the extra images that appear due to the zero padding. See below for more details.
Mix the signal with $\exp(j2\pi7\times 10^6t)$. This shifts the signal to the 7 MHz.
Take the Real component of the result i.e. throw away the imaginary component.

For the upsampling process you have a couple of options. You can upsample by 12 and then decimate by 5. You would up sample by 12, apply the low-pass filter, and then down sample by 5. The filter would need a pass band from 0 - 4MHz (exact upper frequency depends on the nature of your IQ signal). The stop band would need to be at 5 MHz. You might be able to get away with a higher stopband frequency - e.g. if you have no frequency content between 4 MHz and 5 MHz, you could use a stopband starting at 6 MHz. The filter design is slightly dependent on the frequency analysis of your IQ signal and you would have to figure out your requirements for the filter specifications e.g. pass-band ripple and stopband attenuation.

For practical purposes you probably want to break the up sample by 12 into two separate steps e.g. Up sample by 4 and then up sample by 3 - each step would have its own filter applied.

An alternative method is to up sample by a factor of 3 or 4, apply the low-pass filter and then use a spline filter to interpolate the signal to required 24MHz sampling rate.


Estimation of the amplitude of a frequency component in a signal

Generally, the frequency estimation problem decouples from the amplitude + phase estimation problem.

As I said in the comments, you could just do:

$$
\hat{A} = \left|\sum_{n=0}^{N-1} x_n e^{-j 2\pi\hat{f} n} \right|
$$
where $\hat{f}$ is your frequency estimate.

Regarding your issues:

because the number of points isn't high enough to hit the right resolution

I'm not sure what you mean by "resolution" here? Resolving the amplitude?

second is that I need a really fast way

The suggested technique will only take $N$ complex multiplications and $N-1$ complex additions, plus whatever you need to calculate the absolute value (two complex multiplications and one addition).

That's pretty fast, and linear in $N$.  And certainly faster than an FFT (which would be $N\log_2 N$).


Finding the fundamental period by using four elementary operations (+,-,x,/)

Well eventually operations such as the Fourier Transform and autocorrelation are defined with those four operations.

Are comparisons and conditional decision-making allowed?  Are you allowed to pick a pair of numbers from an array based on one of those numbers being a maximum?

My other answer is a place to begin.


Very interesting problem.

Your pulse locations can be found by using a sliding window RMS measurement.  Once found, I would recommend centering a DFT frame on the pulse and applying a VonHann window (broadens narrow peaks and narrows broad peaks), then take the FFT.  Looking at the data, you should get a pretty distinct and clear profile.

Your backgound noise may contain some small peaks too, so you may want to analyze a few "quiet periods" to determine what they are so you can discount them in your evaluation.

Precise frequency determination shouldn't be necessary (or even possible) so a simple estimator, like Jacobsen's or MacLeod's, should suffice.

Like wise, a rough amplitude estimate should also be adequate for each peak.  (See my answer here for phase and amplitude estimation: FFT Phase interpretation of input signal with non-integer number of cycles in FFT window)  Don't use the window if you are going to calculate these.

If you know the frequencies ahead of time (or want to do it twice), selecting a DFT bin that has your frequencies near bins will improve these estimation methods considerably.


The "x" axis you refer to is the frequency axis. To answer your second question first: The frequency axis is typically given as $\omega$ and waveforms as function of $\omega$ can be represented at Real and Imaginary components given by the cosines and sines, or more concisely as Magnitude and Phase as given by the complex exponential since

$$Ae^{j\omega t} = Acos(\omega t) + j Asin (\omega t) = I(t) + j Q(t)$$

Either is a representation of a complex function having real and imaginary components. For representing the frequency spectrum of a time domain function we first map our time domain function to the frequency domain with the Fourier Transform which correlates the time domain function of interest to these basis functions (either cosines and sines or much simpler the complex exponential, either with magnitude = 1). Once in the frequency domain, the result will be complex, so two plots are typically used, since it takes two real numbers to represent one complex number. This can be plots of magnitude vs frequency, and phase vs frequency, or alternatively plots of real component (I) vs Frequency and another of the imaginary component (Q) vs frequency.

When the functions are observed in the time domain (meaning how the function is changing versus time instead of versus frequency as used above) we can represent them as you suggest in terms of sine and cosine components in quadrature. This specifically is a phasor representation; for more detail on that see my explanation on the end where I use the phasor representation to help explain negative frequencies.   (There is no reason we can't use phasors to show functions in the frequency domain as well; it would be a 2 dimensional plot showing the magnitude and phase versus frequency instead of versus time- how it is plotted really comes down to what is the best way to convey the information given the function at hand).

The basis function for the complex exponential spectrum is generally given as $e^{j\omega t}$ as given by the Fourier Transform (which is a correlation to this basis function at any given $\omega$, since correlation in general is multiply and integrate). Generally the Fourier Transform is given as:

$$FT\{g(t)\} = G(\omega) = \int_{-\infty}^{+\infty}g(t)e^{-j\omega t}dt$$

Note since $\omega$ is a continuous function in this context (we are solving for G as a function of $\omega$), there are an infinite number of actual basis functions as the FT is an infinite-dimensional space.

DIRECT SOLUTION, k is the independent variable

If we applied that in this case, with $\omega$ and $k$ as independent variables, we get the following result:

$$FT\bigg\{\frac{1}{2}e^{jk\omega t} + \frac{1}{2}e^{-jk\omega t}\bigg\} = G(\omega)$$

$$= \frac{1}{2}\int_{-\infty}^{+\infty}e^{jk\omega t}e^{-j\omega t}dt + \frac{1}{2}\int_{-\infty}^{+\infty}e^{-jk\omega t}e^{-j\omega t}dt$$

The first term under the integral converges to an impulse (infinite height, area = 1, so 1/2 after being multiplied) when $k = +1$, as does the second term when $k = -1$. For all other k the result is zero.

The amplitudes represent the areas of the impulses given and shows that the Fourier Transform can be given as a function of k alone in this case (since there is no dependence on $\omega$):

$$G(k) = \frac{1}{2}\delta(k+1) + \frac{1}{2}\delta(k-1)$$

We could show this as a surface plot of the magnitude for all values of $k$ and $\omega$ as shown below, with the dashed line representing infinity, which is where an impulse would occur versus k for any fixed value of $\omega$:

Thus for all $\omega$ we would obtain the following magnitude spectrum as a function of $k$ alone:

product of kw as the independent variable

If we decide to make the frequency axis be the product of $k \omega$ (since the frequency of the cosine function is given by that product), we can simply multiply the result given above by $\omega$ to result in the following magnitude spectrum:

We can similarly choose to make $\omega$ be the independent variable and would result in impulses at $\omega = \pm \frac{\omega}{k}$. All of this is in agreement that it doesn't matter what $\omega$ is, in all cases the result is an impulse in the "k domain" when k = $\pm 1$.

Further Background on the FT of the Cosine Function

The following may help provide further intuition for those familiar with the Fourier Transform of a cosine function.

Generally we describe the Fourier Transform of $cos(\omega t)$ as two impulses in the frequency domain $\omega$, but this applies to a specific value of $\omega$, so the time domain function of interest is given as $cos(\omega_o t)$ where $\omega_o$ represents a constant value. The Fourier Transform of $cos(\omega_o t)$ converges to two impulses in the frequency domain, but the Fourier Transform of $cos(\omega t)$ does not converge- it is infinite for all frequencies $\omega$!

What the amplitudes represent is clearest in the case where the frequency term $k\omega$ is held constant, and we define another variable such as $\Omega$ to represent all possible values of $k\omega$, as in that case the Fourier Transform (FT) of a cosine function is simply two impulses in frequency, one at a positive frequency and the other at a negative frequency.

To show this first in general, I will use capital omega $\Omega$ to avoid confusion with the $\omega$ in your formula to represent in general any angular frequency ($\Omega = 2\pi f$), and $\Omega_o$ to represent a specific frequency.

$$FT\{\cos(\Omega_o t)\} = G(\Omega) =  \frac{1}{2}\delta(\Omega-\Omega_o) + \frac{1}{2}\delta(\Omega+\Omega_o)$$

And we can see how this relates directly to the general relationship from Euler's identity expressing the real sinusoidal cosine in terms of two complex exponential frequencies:

$$\cos(\Omega t)= \frac{1}{2}e^{j\Omega_o t} + \frac{1}{2}e^{-j\Omega_o t}$$

And with that we see how when we use the FT to correlate to $e^{j\Omega t}$ the result is two impulses, as when we set the frequency to either of these values $\pm \Omega_o$ the integral given by the FT above goes to infinity since $e^{j\Omega_o}e^{-j\Omega_o}=1$, but would integrate to zero anywhere else.

As a plot vs frequency this would appear as follows showing two impulses that exist at the positive and negative frequency $\Omega_o$. This is a magnitude plot, and typically to be complete a phase plot would also be shown (showing the phase versus $\omega$), but in this case it is trivial since the phase is 0 at both impulses:

Or with $k\omega$ variables, showing the magnitude spectrum that would result with both $k$ and $\omega$ as constants:

This relationship between the frequency and time waveforms as I have shown them applies for when $\Omega_o$ is a constant, otherwise we will be dealing with arbitrary modulation functions and would need to know the modulation waveform specifically (how the independent variable is changing with time) to determine the frequency spectrum.

Negative Frequencies?

To further understand the meaning of a negative frequency, which is applicable to angular frequencies by definition, note that a negative frequency represents a single exponential frequency given as $e^{-j\omega_o t}$, with $\omega_o$ in this case any positive number. It may be helpful to visualize this with reference to the relationship I gave in the first paragraph, showing $e^{j\omega t}$ in terms of I and Q components.  To give negative frequencies a physical interpretation, we can plot the function of $e^{j\omega t}$ as a function of its real and imaginary components versus time resulting in a phasor representation, specifically the plot will map out a unit circle, and rotate around the origin at rate $\omega$ as depicted below. Here we show the angle increasing in a positive direction with time, and thus this represents a positive frequency (A negative frequency would be depicted by instead having the phasor rotate in a clockwise direction):

If this is still not clear, it may be helpful to know that the form $Ae^{j\theta}$ is identical to $A\angle \theta$

If you take the Fourier Transform of a specific exponential frequency with frequency term $-\omega_o$ given as $e^{-j\omega_o t}$, the result is a single impulse at that frequency: $\delta(\omega+\omega_o)$.  (While as we showed above the cosine function has two exponential frequencies; a positive and a negative).

The Fourier transform when presented on a graph with a positive and negative axis represents frequencies in the exponential form, so each impulse shown in the frequency domain is a $Ae^{j\omega t}$ in the time domain. When (and only when) the plot gives positive and negative frequencies where the negative frequencies are the complex conjugate of the positive frequencies given, then in that case we can also represent the same spectrum with just a positive frequency axis; as that is the only way the time domain signal can be real (using a basis function of cosine as you described); if the negative and positive frequencies are not related by a complex conjugate, then the time domain signal must be complex and a positive and negative frequency axis is required to represent the spectrum.

This is clear from the plot below showing both $Ae^{j\omega t}$, and $Ae^{-j\omega t}$ as complex phasors rotating around the origin with time (magnitude A and angle linearly proportional with time as  $\omega t$ , which is by definition a constant frequency since $f = d\phi/ dt$).  The dotted line shows the result of summing those two phasors at any instantaneous time t, showing as long as the two phasors are related by a complex conjugate (equal magnitude opposite phase) the sum will be on the real axis and thus will be a real function with no imaginary components.

Each spinning phasor shown is $e^{j\omega t}$ and $e^{-j\omega t}$ and their summation which stays on the real axis is a (real) cosine function.

Dilation Property of the Fourier Transform

Not specific to your question but worth mentioning since that is occuring here is the dilation property of the Fourier Transform (scaling the independent variable by a constant term k) which has the following general relationship:

$$F\{g(kt)\} = \frac{1}{|k|}G\bigg(\frac{f}{k}\bigg)$$

We could apply the above property directly by transforming the known Fourier Transform of $cos(\omega_o t)$ as another approach to arriving at a solution, although this case would be trivial to do directly as a Fourier Transform on the $cos(k\omega_o t)$ function itself.


Sampling will be at about 1 GHz, output should be about 10 kHz

So, you're decimating by a factor of $10^5=2^5\cdot5^5$.

This very much says that you'd normally go ahead and successively decimate.

In an FPGA, halfband filters are especially cheap and efficient to implement; better yet: they're extremely easy to design, and you don't have to change the shape the halfband filter for each step (though you certainly can, to some gain).

I'd probably (you'll have to do your own effort estimates; Sylvain Muneaut of the Osmocom project had a handy cascaded FIR decimation effort calculator, can't find it...) do the following, from a pure gut feeling

1 GHz -> Halfband
500 MHz -> same Halfband
250 MHz -> ¼-band
62.5 MHz -> simple $\frac15$-band filter
12.5 MHz -> simple $\frac15$-band filter
2.5 MHz -> good $\frac1{25}$-band filter
100 kHz -> excellent $\frac1{10}$-band filter

That way, the easiest to calculate filters (i.e. the ones with the least non-zero coefficients) run at the highest rate, and and the closer you get to the bandwidth you want to "protect", the flatter your passbands and the steeper your transition bands can get.

For complexity reasons, you could consider transferring the 62.5 MS/s or 12.5 MHz out of your FPGA into a PC and processing it there on a multicore CPU – CPU power is much cheaper than FPGA power these days, for memory intense tasks like long FIRs, at least.

Generally, the question is whether you're not doing something questionable when sampling a 10 kHz signal at 1 GHz; but that's due to the physical reality of the system you're observing, and it's hard to help without knowing that system.


There exist perfect reconstruction filter banks that do what you suggest.

One example of the choice of $H_i(z)$ for $i= 0 \ldots N-1$ is the FFT filters (followed by the IFFT filters $F_i(z)$).


Pattern recognition in time series 4x3000 vector

Assuming your time series are the same length, take your data and produce spectrograms, $S_k$, for each row of the $4 \times 3000$ data matrix $D$ that you have. Since these individual time series should be similar, this means you can try to extract features by stacking these spectrograms next to each other horizontally into a large matrix. So if you have $S_k$ as the spectrogram for row $k$, you can create the overall spectrogram as $S = \left[ S_1, \cdots, S_4\right]$.

From here, find the modulus of the elements (since they will be complex), and define the result as $\bar{S} = |S|$. Then I recommend performing a Non-negative Matrix factorization on $\bar{S}$ such that you get the low-rank approximation $\bar{S} \approx FW$ with some small specified number of features. $F$ in this case will be your features (they will be column vectors in $F$) and $W$ will be your new data in feature-space. You can look at $W$ to see where a given feature from $F$ appears. Given you've generated $F$, all you need to do to do to find $W$ for individual/new time series results is compute the spectrogram for that time series, call it $\hat{S}$, and then do $\hat{W} = F^{+} |\hat{S}|$ where $F^{+}$ is the pseudo inverse of $F$. Now you can look at $\hat{W}$ to get an idea of where certain features are being used and also build models using this lower dimensional representation.

Note that the number of features you use will have to be played with. I would look at the modulus of a sample spectrogram to get a feel for how much the frequencies are changing over time and consider roughly how many dominant features you think could be used to approximate it and then set the number of features as that to start. As an example, I have done this stuff above to identify different sounds in an audio clip and I only had to reduce it to the number of dominant sounds (there were 3 of them) so that hugely minimized my dataset size. Images from this example are shown below:

One thing you might consider as well is taking the separate data streams and trying to fuse them using some sort of filtering algorithm or using Hidden Markov Models. Then just perform the above spectrogram computations on the filtered result.


Adaptive Piecewise Constant Approximation (APCA) with wavelets/DWT

Their wavelet coefficients $$W_c = [4.5, 0.5, 1, -1, 1, 1, 0, -1]$$ or $$W_c = [4.5, 0.5, 1/\sqrt{2},-1/\sqrt{2},1/2,1/2,0,-1/2]$$
look perfectly fine to me. In the first case, they use $a = (x+y)/2$ and  $d = x-a$ repeatedly. What might differ is the overall scaling, to account for energy preservation. In Matlab for instance, I get the same result as their by dividing the $3$-level coefficients by $\sqrt{8} = \sqrt{2}^3$:

 wavedec([7, 5, 5, 3, 3, 3, 4, 6],3,'haar')/sqrt(8)

Their algorithm is quite well detailed in the paper, so it should be implementable quite easily. If you want additional help for Matlab like-code (which can be reinterpreted for R or Python), check the tutorial A practical Time-Series Tutorial with MATLAB, starting page 66.


There are direct and indirect effects of distance to the Bit Error Rate (BER) between two transceivers.

The direct effect that distance has on the BER is through the attenuation of the signal. As the transmitter get's further away from the receiver, the power of the received signal is diminished which makes the Signal to Noise Ratio lower which makes signal reception more difficult and the likelihood for a misinterpreted symbol higher.

This is the simplest case with the transceivers operating in free space where theoretically there are no other objects or sources of "noise".

The indirect effect of distance on the BER in a realistic environment (i.e. one containing surfaces that the signal can bounce off, such as buildings) is the introduction of multipath propagation. This means that multiple copies of the same signal arrive at the receiver at slightly different times having followed slightly different paths. Depending on the frequency of the signals, these slight delays may lead to destructive interference. In other words, two versions of the same signal cancel out each other. This condition leads, again, to a decrease in the SNR and signal distortions which leads to a higher probability for errors in symbol reception. A similar (indirect) effect is introduced due to the Doppler shift of the frequencies if the velocity of movement of the transmitter is comparable to the wavelength of the used frequencies.

There are multiple ways by which you can include these effects to a simulation.

When simulating the full transceiver stack, i.e. Source, encoder, modulator, channel, demodulator, decoder, Target then Free Space Loss can be applied at the output of the modulator (as part of the "channel" building block) as an attenuation of the produced signal given a distance L.

Multipath propagation can also be applied at the "channel" stage, either by convolving the modulator's output with a suitable impulse response or by generating an "envelope" signal from a Rayleigh distribution.

For more information, please see this, this (this) and this link.


Blending artifacts in OpenCV image stitching

You need a edge detection algotithm witch will filter the critical edges and then create out of them an binary image. That image will give you an region to blur and the rest should tay the same. This should avoid the water stains or ghost efects.


Frequency change detection in step-wise frequency chirp

I suggest monitoring the phase versus time directly instead of frequency. Frequency is the derivative of phase so the slope of the phase would indicate the frequency. Detrend the phase slope for the starting frequency and then the point in time where the phase starts to ramp up should be easier to detect. The window in which to detect this change will be balanced with the SNR of the signal itself.

Since frequency is the derivative of phase, it will be more sensitive to high frequency noise (the derivative is a high pass function), which is why monitoring phase directly could be a more reliable approach.


Given what you explained, you can count the number of times the sequence is
001, 100 or 010. Then, $\hat{p}$ is defined as that count divided by the number of sequences, $N$.
If we define the null as the hypothesis that the proportion of good sequences  is 1.0,  then the test is
$H_0: p = 1.0$ versus $H_1: p \neq 1.0$
As far as testing the null, ( I'm assuming N is large here so that the normal
approximation to the binomial holds )  you can calculate the test statistic
as $Z =  \frac{ \hat{p} - p_{0}}{\sqrt{\frac{\bar{p}(1-\bar{p})}{N}}}$
and use $p_{0} = 1 $ and $\bar{p} = 0.5$. Then, reject the null if $|Z| > Z_{\frac{\alpha}{2}}$ where $\alpha$ is the critical value of the test.
The CI  would be $\hat{p} ~ \pm Z_{\frac{\alpha}{2}} \times
{\sqrt{\frac{\bar{p}(1-\bar{p})}{N}}}$.
If a digit in one of the three digit sequences is is missing or there's an extra, and you know this, then you can just count that sequence as a bad sequence.


Please consider that the wavelet transform coefficients represent the local details of the signal. In (incomplete) DWT or a Gaussian pyramid, there is also a trend signal after the transform that contains the global behavior.

The first condition on this local-global separation is that constant or near constant parts of the signal get mapped to zero in the local details. The scalar product of a (wavelet) function with a constant function gives its DC contents. So DC of wavelet equal 0 means that the computed coefficients are high at regions of interest like edges.

The same goes for higher moments. If one wants to regard slow linear or quadratic "gradients" in the signal as not essential for the contents, or as part of the trend signal, then the first resp. second moment of the wavelet must vanish.

Or to consider it from a slightly different angle: Each scale of the wavelet transform corresponds to the application of a high-pass filter. One defining property of a high-pass filter is that the amplitude response around the frequency zero is small, ideally zero. Vanishing moments of the wavelet function correspond to a flat graph of the Fourier transform of the wavelet at frequency zero.


Continuous-time mathematical formula for deconvolution filters

The FT of the convolutional inverse is the point-wise inverse in the FD. So the convolutional inverse of a gaussian might have numerical problems. The FT of a gaussian is a gaussian. Its point-wise inverse thus can be large in many places.

For a box function, the situation doesn't look all that good either. Its FT is a variant on $\operatorname{sinc}$, which has many zeroes and loci of low magnitude.

You can certainly determine the convolutional inverse of a gaussian, but the inverse of a box function is not even a function. Convolving with a box will cancel any function that is periodic with a period multiple of the edge length.


how to create signal with constant $f_1$ frequency after T?

Find the instantaneous phase at the end of the $T$ interval and keep adding it to the phase that makes the $\sin()$ "go round".

A plain simple oscillator at some frequency $f$ has its phase increasing at a constant rate.

Here is an example (in GNU Octave but easily transferable to other platforms as well):

Fs = 44100; % Sampling frequency in Hz
T = 1; % Total duration of the signal in seconds
f = 120; % Frequency of the oscillator
t = 0:(1./Fs):(T - (1./Fs)); % Time vector
p = 2.0 .* pi .* t; %Phase vector
y = sin(f.*p);

Now y contains our 120 Hz sinusoid sampled at 44.1kHz.

Notice here that phase (p) is a "straight line" that grows at a constant rate of $\frac{2 \pi f}{Fs}$.

When the oscillator is "chirping", the rate of change of the phase is variable. For example:

Fs = 44100; % Sampling frequency in Hz
T = 1; % Total duration of the signal in seconds
f0 = 1; % Start chirp at
f1 = 120; % end chirp at (Better keep f1>f0)
t = 0:(1./Fs):(T - (1./Fs)); % Time vector
p = 2.0 .* pi .* t; %Phase vector
y = sin(f0.*p + 2.*pi.*(((f1-f0)/(2.*T)).*t.^2));

Now, this looks like:

And the rate of change of the phase of the $\sin$ is:

Now, from the first example, the phase "step" is f.*p(2)-f.*p(1) == 0.017097 and from the second example, after doing:

s = diff(f0.*p + 2.*pi.*(((f1-f0)/(2.*T)).*t.^2)); % Get the first derivative of phase to find the **rate of change**.

The "last" rate of change is s(end) == 0.017097.

So, to have the oscillator stay at the chirp's concluding frequency, keep accumulating the phase at the constant rate you are after from the last known value of the chirp's phase (to avoid phase jumps). For example:

r = f0.*p + 2.*pi.*(((f1-f0)/(2.*T)).*t.^2);  % The phase as above.
z = r(end) + cumsum(ones(1,T.*Fs).*0.017097); % A running sum of duration T.*Fs (samples).

Now, if you feed the combined result of r, z into the oscillator y with something like y = sin([r,z]); it would look like:

...with a rate of change of:

Hope this helps.


Which algorithm could be used to compare two waveforms from an oscilloscope in XY mode?

You can't get the exact value because you don't know how log the signals spent at each stage.
Consider the case where you have two square with period 1, and shifted by $\delta$ compared to each other.
The XY plot will always be a rectangle with the width corresponding to the amplitude of the X signal and height corresponding to the amplitude of the Y signal. In other words, you have always the same input for your function to compute the correlation, while the correct output varies.
You still can get an idea of how similar two signals are, but it is not the signal cross correlation


What is the best approach to classify an image containing either text or image?

The segment with text (I assume with a background) should contains text clearly visible from the background, then the text segment should have a lot of peaks in the transitions from text to background. The segment icon normally shouldn't have many peaks because in general there are few variations in an image. So if you make a stft of a part of a segment the text should have more harmonics than the icon.
However if the icon is similar to the text or if it has many peaks this idea can't works. You can also find the color density, a single color text on a single color background has the color density concentrated in two points. This is only an idea and it must to be improved.


Hi: I think there are two different vewpoints to your question so I'll mention both.

First Viewpoint: ( Time Domain )

What you have written for the covariance is not wrong but you need to understand what you're calculating. Those numbers you are calculating are unscaled covariances of an MA(n) process  assuming A) that the variance of the noise term in the  MA(n) model, $\sigma^2$ = $1$ and B) the coefficients of the MA(n) model ( i.e: the $\theta$ ) are all 1.   So, the covariance of the shifted MA(n)series does stay the same because the 5 numbers ( A) the variance and B) the 2 autocovariances at lags 1 and lag 2 and lag -1 and lag -2 which are happily symmetric ) are staying the same. You get the other zeros because, for some reason, when you shifted the sequence in the time domain, some extra zeros arose. I'm not sure how that happened but it shouldn't. atleast not in the time domain.

Second viewppoint ( Frequency domain )

I'm just learning DSP so can't say a lot here but, from a time series frequency domain perspective,  if $y$ has  covariance generating $g_{y}(z)$, then that $X(z)$ filter can be viewed as an MA filter of order(n) so that if the MA(n) filter is applied to $y(t)$, resulting in say, $X(t)$, then the covariance generating function of $X(t) = B(z)B(z^{-1}) g_{y}(z)$ where $g_{y}(z)$ is the autocovariance generating function of $y$.

All of this quite nicely covered in Hamilton's "Time Series Analysis" if you have that book. I highly recommend it for time-series but it's mostly time domain material with frequency domain coverage somewhat scattered. Still, book is really great.


normalize fft produced by fftw_plan_dft_r2c_1d

Reading the manual, even though only half the data is returned, the full length FFT is calculated:

So you should normalize by whatever the full length is normalized to: fft_size.


Images in focus have sharper edges, so applying an edge detector and measuring the energy of the output puts you on the right track.

A simple technique is to compute the sum of the laplacian over the images, possibly center-weighted if you deal with everyday photographs (people tend to put the object of interest in the center, so it's better to have the center in focus than the borders)

This criterion is for example what is used in many "Stack focusing" applications for macro photography.

There are other approaches used in de-hazing or tone-mapping... Search for "Local contrast measures".


Fitting new images from a SVD/PCA calculation

The "projection" that is referred to is a vector projection. To calculate the projection of vector $\mathbf{a}$ onto vector $\mathbf{b}$, you use the inner product of the two vectors:

$$
\mathbf{a_{proj}} = \langle \mathbf{a}, \mathbf{b} \rangle \mathbf{b}
$$

$\mathbf{a_{proj}}$ in this case is the vector component of $\mathbf{a}$ that lies in the same direction of $\mathbf{b}$. In Euclidean space, the inner product operator is defined as their dot product:

$$
\langle \mathbf{a}, \mathbf{b} \rangle = \mathbf{a} \cdot \mathbf{b} = \sum_{i=1}^{n}a_ib_i
$$

where $n$ is the number of components in the vectors $\mathbf{a}$ and $\mathbf{b}$ and $a_i$ and $b_i$ are the $i$-th component of vectors $\mathbf{a}$ and $\mathbf{b}$, respectively. Intuitively, by calculating the inner product of the two vectors, you find "how much of" vector $\mathbf{a}$ goes in the direction of vector $\mathbf{b}$. Note that this is a signed quantity, so a negative value would mean that the angle between the two vectors is greater than 90 degrees, as illustrated by an alternative definition for the projection operator:

$$
\mathbf{a_{proj}} = |\mathbf{a}| \cos(\theta) \mathbf{b}
$$

where $\theta$ is the angle between the two vectors.

So, given a vector $\mathbf{a}$ and a bunch of basis vectors $\mathbf{b_i}$, one can find "how much of $\mathbf{a}$" goes in each of the directions of each of the basis vectors. Typically, those basis vectors will all be mutually orthogonal. In your case, the SVD is an orthogonal decomposition, so this condition should be satisfied. So, to accomplish what you describe, you would take the matrix of eigenvectors $\mathbf{U}$ and calculate the inner product of the candidate vector $\mathbf{y}$ with each of the matrix's columns:

$$
p_i = \mathbf{y} \cdot \mathbf{u_i}
$$

The scalar value $p_i$ that you get from each inner product represents how well the vector $\mathbf{y}$ "lined up" with the $i$-th eigenvector. Since the eigenvectors are orthonormal, you could then reconstruct the original vector $\mathbf{y}$ as follows:

$$
\mathbf{y} = \sum_{i=1}^n p_i \mathbf{u_i}
$$

You asked whether this representation is unique; I'm not sure exactly what you mean, but it is not unique in the sense that a given vector $\mathbf{y}$ could be decomposed by projection onto any number of orthonormal bases. The eigenvectors contained in the matrix $\mathbf{U}$ are one such example, but you could use any number of others. For instance, calculating the discrete Fourier transform of $\mathbf{y}$ can be viewed as projecting it onto an orthonormal basis of complex exponential vectors of varying frequency.


Unable to understand SNR equation in Matlab

To express power $P$ in $\mathrm{dB}$ one performs this operation:

$$P_{\mathrm{dB}}=10\log_{10}(P)$$

The inverse operation would then be

$$\frac{P_{\mathrm{dB}}}{10}=\log_{10}(P) \implies 10^{\frac{P_{\mathrm{dB}}}{10}}=P$$

And that's exactly what you are doing in that MATLAB equation (i.e. expressing the SNR in linear scale from the data in decibels).

Note that the .^ means that the operation is being carried out element-wise.


How to get relative rotation matrix from two orientation values in android?

Your error might be due to the indices used. Nevertheless, it is always better to state the proper procedure to do it, so that you could find your error.

Let us assume that a 3D scene point $P$ is related to the normalized coordinates of the camera $i$ with the extrinsic parameters $P'_i=R_iP+T_i$ and similarly to camera $j$ with $P'_j=R_jP+T_j$. To find the relative pose we could simply use these two equations to go from the camera coordinates of view $i$ to the real scene and then to the second view $j$. We could then write:

$P=R_i^{-1}(P_i-T_i)$ (from camera $i$ to world)

$P'_j = R_j (R_i^{-1}(P_i-T_i)) +T_j$.  (from world to camera $j$)

If we re-arrange:

$P'_j = R_{rel}P_i +T_{rel}$, where $R_{rel}=R_jR_i^{-1}$ and $T_{rel}=-R_{rel}T_i+T_j$.

where the subscript rel refers to "relative" and is what we seek to find. Keep in mind that, as rotations are orthogonal, $R_i^{-1}=R_i^T$. You can always compute relative orientations in similar manners.

Throughout the computation we used the convention that the transformation $J=[R|T]$ is a mapping from scene to the camera (the inverse is hence the vice versa). The final pose transforms a point in camera $i$ to camera $j$.


What are some techniques to recover an original image from three different compressed versions?

First of all, there are many solutions. Something better under one measure may be worse under another measure. So first you need to think about what is your quality measure(s) to evaluate a restored image.

Commonly, people use the "mean square error (MSE)" (or its log version, known as the peak signal noise ratio). Assume the ground-truth image is $X$, and the recovered image is $\hat{X}$. Then MSE is defined as
$${\rm{MSE}}(X,\hat{X}) = \sum_{i,j}(X[i,j]-\hat{X}[i,j])^2/N$$
where $N$ is the number of pixels in $X$, and $X[i,j]$ denotes a pixels located at $i$th row and $j$th column in $X$.

Without loss of generality, say we have $n$ compressed counterparts are $C_1,C_2,\cdots,C_n$, and we reconstruct $\hat{X}$ as
$$\hat{X}=f(C_1,C_2,\cdots,C_n)$$
where $f(\cdot)$ is a method to recover the original image from its compressed counterparts. And your objective is to find an appropriate function $f(\cdot)$ and minimizes the MSE.

To simplify the problem, people often restrict the function $f(\cdot)$ to some fixed type but with parameters to be learned from training data. For example, one simplest function can be the linear combination i.e.
$$f(C_1,C_2,\cdots,C_n) = \sum_{k}w_kC_k$$
where $w_k$ is the weight of using the $k$th compressed image in the recovered image. This is a classic regression problem and can be solved easily (go and check the normal equation, which gives the closed form solution to find all $w_k$s). However, you should note the following things:

1). if you donot have any training data, the above method does not
work
2). if your image is very large or your have many training images,
computing the closed form solution is not a good idea, because of too
many computations. Instead, you shall seek for methods like the
stochastic gradient descent or similar ones.
3). if I were you, I will go and read more about linear regression
stuff. You will see interesting topics like how to prevent overfit by
adding regularization, and how to include the Bayesian idea in this
regression problem.
4). my guess says the linear regression method is not the best way,
but it is a very good start point. Later your may try other
$f(\cdot)$, like kernel regressions, or adding more terms, e.g. the
total variation term.


What is difference between frequency and pitch?in nutshell and simple words?

Frequency is mathematically defined as the number of cycles per second. So it is a more strict word mathematically. It is represented numerically by the unit called Hertz.

$f=1/T$, where $T$ represents the one-period length of a waveform. This makes frequency quantifiable.

Pitch on the other hand, is a perceptual characteristic of a sound frequency, so it's a word used to describe how humans perceive a tone as different or not. That's why it is mostly used in the field of music more than mathematics. It is closely associated to the frequency of a sound, but not all frequencies can by described as different pitches (e.g. 440 Hz and 441 Hz would practically sound the same pitch, but a 440 Hz and 466 Hz would sound as two different pitches).

That's how musical scales were invented---to create a scale that classifies the fundamental frequencies and label them as distinct pitches (or notes). A good musical scale represents the foundation of music theory for different styles of music.


How to take the linear convolution of these two signals?


Please help me understand this FM carrier frequency inequality

Checking the book you linked to, this equation refers to FM modulation over a (weakly) non-linear channel, modeled as a memoryless third-order non-linearity (Eq. 82). The third-order non-linearity produces signals at twice and at three times the carrier frequency $f_c$ (Eq. 84). Note that also the frequency deviations of these additional signal components at $2f_c$ and at $3f_c$ are doubled and tripled, respectively.

If $\Delta f$ is the frequency deviation of the original FM signal, the signal component at $2f_c$ has a frequency deviation of $2\Delta f$, and according to Carson's bandwidth rule, that signal component extends down to approximately

$$f_1=2f_c-2\Delta f-W\tag{1}$$

where $W$ is the bandwidth (highest frequency) of the message signal. The original FM signal centered at $f_c$ extends up to approximately

$$f_2=f_c+\Delta f+W\tag{2}$$

In order to be able to separate the original FM signal from the signal component at $2f_c$ created by the non-linearity, we require $f_2<f_1$, i.e.,

$$f_c+\Delta f+W<2f_c-2\Delta f-W\tag{3}$$

or, equivalently,

$$f_c>3\Delta f+2W\tag{4}$$

which is the same as Eq. $(85)$ in the book.


Intro Question to Signal Processing (Low-Pass Filter)

Based on your description it sounds like the idea is to use the FFT to figure out where your signal energy is in terms of frequencies.  Once you know that you know what the passband of your filter should be, and from that you can decide on a reasonable cutoff frequency.  You want to cutoff as much of the noise energy as possible which means you want to have the cutoff frequency as close to the passband frequency as possible, but you need some space for the transition band to make the filter implementable with a reasonable number of taps.

By the way, I often like to look at the signal energy in the log scale because it helps to get a better understanding of the weaker portions of the signal.  In Matlab you do it like so:

plot(20*log10(abs(fft(signal))))


Auto/cross correlation of a cosine in the context of an adaptive beamformer

The $90$ degree shift introduces the sine component to the cross-correlation vector $P$ at lag 1 ($\phi_{dx}(1)$).

The block should be considered to be a 'high-level' description, and does not mean that a shift of one sample is equivalent to  a $90$ degree phase shift, or anything like that. Rather it should be thought of as "shift by as many samples as necessary for a phase shift of $90$ degrees".


What distance metric can I use for comparing images?

The following is not intended to be an answer, but is a statistic that will help us choose an appropriate image comparison technique based on the characteristics of the images you are analyzing.

The first step is to plot a "delta histogram" as follows:

for (x,y) in [0, width] x [0, height] begin
    delta = abs( SecondImage(x, y) - FirstImage(x, y) )
    hist[delta] += 1
end

Given a plot of this histogram, we will know a bit more about the "magnitude" of the changes you are looking for, and will make better recommendations.

(Alternatively, post some sample images. Remember that if the sample images aren't representative of the image differences you are interested in, we might make inferior recommendations.)

You can also test Structural Similarity (SSIM) on your image set and post your results here. Remember that SSIM is designed to mimic human's ability to recognize obstructiveness of image degradation, so it would detect pixelation but maybe not blurring.

If your images are not photographic images (or, are scientific images that are not ordinary subjects of photography), then please also post examples of their 2D autocorrelation, suitably cropped and scaled.

Face recognition is too big a topic to be discussed in a single question. Blurring arises in multiple context in face recognition - it can be a data quality issue, or it can be done intentionally as an intermediate step in data processing.

In face recognition we want to detect the identity of faces, therefore we have to ignore image differences that are not caused by identity differences. The basic category of differences that should be ignored in face recognition are: pose, illumination and facial expression.

A general approach to ignore irrelevant differences is called normalization, which attempts to apply various operations and transforms on the input image to obtain a "canonical" or "preprocessed" image, which in turn can be used for identification.

A second approach is to extract features from images that are highly-invariant from the irrelevant factors.

The quality of a face image is subject to the capturing device and the environment where it was captured. When a face image is captured without cooperation of the subject (such as from a security camera), poor image quality is an unavoidable consequence and had to be remedied by software so as not to hamper identification.

In cooperative capture, a computerized measure of image quality is good: the operator can be notified of quality problems and the image can be re-taken.

Blurring can also be an example of malicious tampering of biometrics in order to evade detection (along with occlusion and disguise). If the image is encoded digitally, a digital checksum and cryptographic signature is sufficient to solve the problem completely. If the blurred image is submitted in physical print by an impersonator, a computerized measure of facial image quality can be used to reject such submissions.

The lack of 2D-localizable features or interest points in a certain part of facial image can be a sign of intentional blurring.

However, the broad category of digital image tampering (by a skilled user of image-editing software), can only be dealt with Digital image forensics which compares pixel statistics against known camera models.


Periodic signals in Continuous and discrete time

A periodic continuous-time signal satisfies $x(t)=x(t+T_0)$ for all $t$. The period $T_0$ doesn't need to be a rational number. A periodic discrete-time signal satisfies $x[n]=x[n+N]$ for all integers $n$. The period $N$ is an integer.
If you sample a periodic continuous-time signal, you don't necessarily get a periodic sequence. E.g., sampling the periodic signal $x(t)=\sin(\omega_0t)$ results in the sequence $x_d[n]=x(nT)=\sin(\omega_0Tn)$, which is only periodic if $\omega_0T=2\pi k/l$ with integers $k$ and $l$.
In general, sampling a $T_0$-periodic continuous-time function with sampling interval $T$ results in a periodic sequence only if
$$kT_0=lT,\qquad k,l\in\mathbb{Z}$$
is satisfied.


Does a power spectral density that follow a power law imply that the time hisory follows a power law?

should imply that also the time history has a power law

No. The PSD is an incomplete description of the signal since it has lost the phase information. There are infinitely many time sequences that have the same PSD. For example a unit impulse, white noise and a linear sweep have the same PSD (roughly speaking) but are completely different time-domain wave forms.
In your case $u(t) \cdot t^{\alpha}$ is one signal that fits your PSD but there are many others that will fit too but have a totally different shape.
This being said, your observation of a conflict has some merit. Your first equation simply states that $g(t)$ is mean free which in turns implies that the PSD at $f=0$ must be zero, but $0^{-5/3}=1$. I'm guessing that your power law applies only to a limited frequency range which doesn't include $f=0$ and hence there is no conflict.


How to extract white objects in an image(attached)? The program should then display each N as a separate figure.

It's all about detecting intensities of pixels. It'll be better if you use MATLAB. Such binary images will have just 2 intensity levels 0(black) and 1(white).

No one is going to do it for you unless you are ready to pay. Buy and read a copy of Digital Image Processing by Gonzalez


Convolution involving turning each pixel value to the most represented pixel value of the neighbourhood

This is a non-linear filter, the operation being performed is not a convolution and cannot be represented by a filter kernel. It is sometimes called a mode filter, by analogy to the median filter, since the value taken by a pixel is the mode of the distribution of the neighboring values. OpenCV does median filtering (cvSmooth with CV_MEDIAN parameter), I don't think it does mode filtering though.


There is no way to regain any lost data at all. However, there are many phenomena that might distort the image in a recoverable way, not actually losing any data. You seem to try to undo a "point spread function", that is eg. some consistent blur. This would be not the only way an image can be distorted, however it is by far the most common.

In that case you would assume that the original image is convolved or filtered by an (initially unknown) point spread function. There may be the inverse operation to that. You could call this "equalization" of the frequency space of the image, then.

Not every inverse operation would be represented by an "edge sharpening" operation that is known a priori. You could maybe try to look at the image a virtually point-shaped light source gives, preferrably in front of a dark background. A laser pointer or a bright LED might serve as an example.


Explanation of Spectral Line Splitting when AR process is overmodeled

One way of doing power spectrum estimation of random processes actually requires finding a set of ARMA model parameters that best approximates the given observed random process $x[n]$.

In your case, you selected an all-pole signal model with $p$ poles as the roots of the denominator polynomial $A_p(z)$ with coefficients $a_p[k]$, and a $b(0)$ related with power scaling of the process. Each pole represents a peak in the spectrum for $0 \leq \omega \leq \pi$.

So, PSD estimation is reduced to finding (estimating) all-pole model parameters $a_p[k]$ and $b[0]$, from the given data by first estimating the auto-correlation sequence values $r_x[k]$ for $k=0,1,2,...,p$  , and then solving a set of linear equations (Prony or Yule-Walker etc) to find those $a_p[k]$ and $\epsilon_p$.

Note that there are a number of approaches in estimating the ACS $r_x[k]$ of $x[n]$, such as Auto-Correlation, Covariance, Burg's, Modified Covariance, etc...

Certain of these techniques (such as the autocorrelation and Burg's method) produce an artifact known as spectral splitting when they are forced to fit a high order model to an actually low order data; aka a model-data mismatch problem. There is not much you can do to prevent it. I don't know an analytic derivation of the relation between the chosen model order $p$ and the number of splitted peaks you will observe. That's also dependent on the observed noise instance being processed. So that number of peaks is also random.

However, it's stated that modified covariance method of ACS estimation can overcome the spectral splitting problem.

Furthermore, if you continue reading, you will see in the section named Selecting the Model Order (p. 445, Eq.8.23) that by trying to make a good guess of model order you may reduce the chances of generating that line splitting problem.


Capture legitimate huge increase/drop in fuel level from noisy measurements

Presumably, fuel theft and re-fueling require the vehicle to be stopped and its engine turned off (?). Every time that you start the vehicle up, re-set the system to be tracking the fuel level in the tank at start-up. Comparing the "fuel level at (previous) shut down" to the "fuel level at start up" could tell you what is going on.

Sloshing can persist of course, even during brief stops. Adding a low pass filter / integrator to the sensor would improve instantaneous readings.

It would be interesting to obtain the Fourier Transform of the fuel level sensor signal while the vehicle is moving, provided that it is sampled quickly enough. Depending on your specifications for fuel level accuracy, you might find that a simple low pass filter / integrator might be equally effective in "dampening" the sloshing effect on the signal. This is how common "anti-slosh" modules work anyway.

Hope this helps.


Kronecker delta instead of dirac delta as correlation function of white noise

Dirac delta function has a continuous argument, but Kronecker delta function has a discrete argument.
Your example is a discrete signal so Kronecker delta is used.


Plotting the magnitude response of a filter?

If you want to code your own application that plot the magnitude response, you first need to extract the poles and zeros from your transfer function in the $Z$ domain. The process that follows can either be analytic or graphical. I will try to cover both, starting with the analytic approach, then graphical.

Extracting the poles and zeros

Taking your time domain equation. :
$$y[n] = 0.0976⋅x[n]+0.1952⋅x[n−1]+0.0976⋅x[n−2]+0.9429⋅y[n−1]−0.3334⋅y[n−2]$$
Transfer function in the Z domain is
$$H(z)=\frac{Y}{X}=\frac{0.0976+0.1952z^{-1}+0.0976z^{-2}}{1-0.9429z^{-1}+0.3334z^{-2}}=\frac{0.0976(1+2z^{-1}+1z^{-2})}{1-0.9429z^{-1}+0.3334z^{-2}}$$

Solving the numerator and denominator for $z=0$ will yield some poles,zeros and gains. I will not detail this step as this subject is widely covered. In your case, you will find:
$$zeros = \{-1, -1\}$$
$$poles = \{(0.4746 + 0.3289j), (0.4746 - 0.3289j)\}$$
$$K_{num}=0.0976$$
$$K_{denom}=1$$
The analytic approach

Once you have the poles and zeros, you can rewrite your transfer function in this different form :

$$H(z)=\frac{K_{num}}{K_{denom}}\cdot\frac{(z-zero_{1})(z-zero_{2})...(z-zero_{n})}{(z-pole_{1})(z-pole_{2})...(z-pole_{n})}$$

Or in a condensed form :
$$H(z) = \frac{K_{num}}{K_{denom}}\cdot\frac{\prod_{n=0}^{N_{zeros}}{(z-q_{n})} }{\prod_{n=0}^{N_{poles}}{(z-p_{n})}}$$

The magnitude response of your filter is basically the magnitude of your transfer function when $z=e^{j\omega}$. We can define  $|H(z)|\biggr\rvert_{z=e^{j\omega}}$

You then get :
$$|H(z)|\biggr\rvert_{z=e^{j\omega}} = \biggr\rvert \frac{K_{num}}{K_{denom}}\cdot \frac{\prod_{n=0}^{N_{zeros}}{(e^{j\omega}-q_{n})} }{\prod_{n=0}^{N_{poles}}{(e^{j\omega}-p_{n})}}\biggr\rvert$$

Translate that into a code, and you get something like : (matlab example)

h = ones(1,n);
knum=0.0976
kdenom=1
%Zeros
for k=1:length(q)
    v = exp(1i*w)-q(k);
    h = h .* v;
end
%Poles
for k=1:length(p)
    v = exp(1i*w)-p(k);
    h = h ./ v;
end
h = h * knum/kdeom
plot(abs(h))

The graphical approach

What we will see here, is exactly what we just saw in the analytic approach, but we will try to visualize it a little bit. Let's plot you poles and zeros into the $Z$ plane:

The unit circle, or $z=e^{j\omega}$, contains all the frequencies from $\omega=0$ to $\omega=Nyquist=\frac{2\pi f_{s}}{2}$

In order to know the frequency response of your filter at a specific value of $\omega$. Draw a line from each poles/zeros to the corresponding point on the unit circle.

Take the products of the line length originating from a zero and divide by the product of the line length originating from a poles. You'll get the magnitude response of your filter.

In other word:

$$ |H(z)|\biggr\rvert_{z=e^{j\omega}} = \frac{d_{3}d_{4}}{d_{1}d_{2}}\cdot\frac{k_{num}}{k_{denom}}$$

Take few seconds to understand what we just did there, and you will see that it is exactly what the analytical approach suggests.

Finally

I wrote some matlab code to plot the frequency response of a filter not so long ago and I posted this question to get help getting this done. It might also help you.

Good luck


Different results applying DFT (formula only) and FFT on Matlab?

The Fast Fourier Transform (FFT) is a class of algorithms that efficiently implement the Discrete Fourier Transform (DFT). However, what you computed by hand is the Continuous-Time Fourier Transform (CTFT), which is substantially different from the DFT. The DFT is applied to finite length discrete-time sequences, whereas the CTFT is applied to continuous-time functions of possibly infinite duration:

$$\begin{align}X[k]&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/N}\quad&\text{(DFT)}\\
X(j\omega)&=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\quad&\text{(CTFT)}\end{align}$$

However, you can use the DFT to approximately compute the CTFT, but there will always be at least one of the following two errors (and usually both):

Truncation error: if the signal is too long (or has infinite duration as in your example), it must be truncated to fit inside the DFT window.
Aliasing error: the continuous-time signal must be sampled, and sampling will usually introduce aliasing unless the signal is band-limited and the sampling theorem is satisfied.

Since the signal in your example is neither band-limited nor time-limited you will get both errors mentioned above.

More information on approximating the CTFT by the DFT can be found in this answer.

The following Octave script shows how the approximation works for the example in your question.

% choose parameters (sampling frequency and time interval)
Fs = 50;        % sampling frequency
T1 = 2;         % lower and upper
T2 = 20;        % integration limits

Ts = 1/Fs;
DT = T2 - T1;
N = round(DT/Ts);   % DFT length

% sampled time-domain function on [T1,T2]
t = (0:N-1)*DT/(N-1) + T1;
xd = (t-2).*exp(-t);

% DFT approximation of CTFT
Xd = Ts * fft(xd);
fd = (0:N-1)*Fs/N;

% phase correction for offset T1
Xd = Xd .* exp(-1i*T1*2*pi*fd);

% exact expression for CTFT
N2 = round(N/2);
fc = fd(1:N2);
Hc = exp(-2*(1+1i*2*pi*fc))./(1+1i*2*pi*fc).^2;

% plot results
subplot(2,1,1), plot(fc,20*log10(abs(Xd(1:N2))),fc,20*log10(abs(Hc)))
    title('magnitudes (dB)'), xlabel('f'), legend('DFT','CTFT'),
    axis([0,Fs/2,-110,-10]), grid on
subplot(2,2,3), plot(fc,20*log10(abs(Xd(1:N2))) - 20*log10(abs(Hc)))
    title('magnitude error (dB)'), xlabel('f'), grid on
subplot(2,2,4), plot(fc,unwrap(angle(Xd(1:N2))) - unwrap(angle(Hc)))
    title('phase error (rad)'), xlabel('f'), grid on

The plot below shows the magnitudes of the CTFT and its approximation by the DFT, as well as the magnitude and phase approximation errors.


What is zero-phase filtering and forward-backward filtering?

Zero-phase filtering is a non-causal procedure, so it cannot be done in real time, only offline (or pseudo real-time, i.e., with a sufficient delay). A zero-phase filter needs to have a purely real-valued frequency response, and, consequently, it must have an impulse response that is even with respect to the time index $n=0$, i.e., it is non-causal.

Zero phase filtering with IIR filters is achieved with forward-backward filtering, as implemented in Matlab's filtfilt function. The resulting total frequency response is the squared magnitude of the original IIR filter's frequency response. Since the squared magnitude is real-valued, the resulting filter is a zero-phase filter. Of course there's a delay because you have to feed the time-reversed output of the first filter pass back to the input of the filter, after which you need to time-reverse the output. More details about forward-backward filtering can be found in this answer.


Can we normalize data to unit variance without distorting the phase?

"Normalizing the variance": Since variance is a real number, normalizing it means multiplying the signal with a real number.

That doesn't change the phase, ever.


Gesture recognition using only simple techniques

You can try looking at Hu invariant moments. They can be constructed from basic moments, and are rotation, scale, reflection, and translation invariant.

Calculate them for a set of training contours first, and then apply them to the test contour.

There are implementations in Matlab and OpenCV, as far as I remember.


Let's start the right way (as Tim's comment says): let's start with the signal model of what you're doing.
Let's assume that the measurements at time $k$ are just positions $x_k$ and $y_k$ and that these are also part of the state, which also includes $\dot{x}_k$, $\dot{y}_k$, $\ddot{x}_k$, and $\ddot{y}_k$:
$$
\mathbf{x}_k = \left[   x_k \ \dot{x}_k \  \ddot{x}_k \ y_k \  \dot{y}_k \
 \ddot{y}_k   \right]^T
$$
Then the output equation is
$$
\mathbf{z}_k = \mathbf{H} \mathbf{x}_k + \epsilon_k = \left [ x_k\ y_k \right ]^T + \epsilon_k
$$
with
$$
\mathbf{H}  = \left [  \begin{array}{cccccc} 1 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 0 & 0
\end{array} \right ]
$$
and $\epsilon_k$ is the measurement noise and it is distributed jointly normally with zero mean and covariance $\mathbf{R}$ (called the observation_covariance in pykalman). The observation_matrices is what $\mathbf{H}$ is called in pykalman.
The state update equation will then just be
$$
\mathbf{x}_{k+1} = \mathbf{A} \mathbf{x}_k + \mathbf{B} \eta_k
$$
where $\mathbf{A}$ is the state transition matrix and $\eta_k$ is the process noise.
For a standard displacement, velocity, acceleration model driven by noisy jerk inputs, this means:
$$
\begin{align}
x_{k+1} &=& x_k + \Delta T \dot{x}_k + \frac{1}{2} \Delta T^2 \\
\dot{x}_{k+1} &=& \dot{x}_k + \Delta T \ddot{x}_k \\
\ddot{x}_{k+1} &=& \ddot{x}_k + \eta_k
\end{align}
$$
with similar expressions for $y$.
That means
$$
\mathbf{A} = \left[
\begin{array}{cccccc}
1 & \Delta T & \frac{1}{2} \Delta T^2 & 0 & 0 & 0\\
0 & 1 & \Delta T& 0 & 0 & 0\\
0 & 0 & 1& 0 & 0 & 0\\
0 & 0 & 0 & 1 & \Delta T & \frac{1}{2} \Delta T^2\\
0 & 0 & 0 & 0 & 1 & \Delta T\\
0 & 0 & 0 & 0 & 0 & 1
\end{array}
\right]
$$
and
$$
\mathbf{B} = \left [
\begin{array}{cc}
0 & 0 \\
0 & 0 \\
1 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 1 \\
\end{array}
\right ]
$$
Then the process noise covariance is $\mathbf{B}\mathbf{B}^T \Sigma_\eta$ where $\Sigma_\eta$ is the covariance of $\eta_k$. This is called the transition_covariance in pykalman.
If I do that and take some guesses for the process and measurement noise scales, then I get the result below.

As you can see, the early part of the plot gets perturbed by missing data.  I'm not quite sure the best way to deal with that.
New Code below
import numpy as np
import pandas as pd
from pykalman import KalmanFilter
from numpy import ma
import matplotlib.pyplot as plt

# Reading data
observations = pd.read_csv('Q66687.csv')
rearwheel_y = observations['rearwheel_y'].to_list()
rearwheel_x = observations['rearwheel_x'].to_list()

# Locate and clip relevant snippet of motion, removing unecessary data at the beginning and end
rearwheel_y = rearwheel_y[385: 600]
rearwheel_x = rearwheel_x[385: 600]

#rearwheel_y = np.sin(20*time) + 10

# The numpy.ma module provides a convenient way to address the issue of dropout, with masked arrays.
# When an element of the mask is False, the corresponding element of the associated array
# is valid and is said to be unmasked.
# When an element of the mask is True, the corresponding element of the associated array
# is said to be masked (invalid).
rearwheel_x = ma.masked_values(rearwheel_x, 0)
rearwheel_y = ma.masked_values(rearwheel_y, 0)
measurements = np.array([rearwheel_x, rearwheel_y])

# initial_state_mean
Xbar = np.array([rearwheel_x[0],0,0,rearwheel_y[0],0,0]).transpose();
# initial_state_covariance
P0 = 100*np.identity(6)
n_timesteps = len(measurements[0])
n_dim_state = 6
filtered_state_means = np.zeros((n_timesteps, n_dim_state))
filtered_state_covariances = np.zeros((n_timesteps, n_dim_state, n_dim_state))

# Kalman-Filter initialization
DeltaT = 0.01668335
A = np.array([[1, DeltaT, 0.5*DeltaT**2, 0, 0, 0],
              [0, 1, DeltaT, 0, 0, 0],
              [ 0, 0, 1, 0, 0, 0],
              [0, 0, 0, 1, DeltaT, 0.5*DeltaT**2],
              [0, 0, 0, 0, 1, DeltaT],
              [0, 0, 0, 0, 0, 1]])
B = np.array([[0, 0],
             [0, 0],
             [1, 0],
             [0, 0],
             [0, 0],
             [0, 1],])
SigmaEta = np.array([[1,0], [0,1]]) # Assume eta covariance is unit...
#Q = B @ B.transpose() * 0.0001
Q = np.identity(6) * 0.00001
R = np.array([[1,0],[0,1]])*0.00005
H = np.array([[1,0,0,0,0,0],[0,0,0,1,0,0]])
kf = KalmanFilter(transition_matrices = A,
                  observation_matrices = H,
                  transition_covariance = Q,
                  observation_covariance = R,
                  initial_state_mean = Xbar,
                  initial_state_covariance = P0)
# iterative estimation for each new measurement
for t in range(n_timesteps):
    if t == 0:
        filtered_state_means[t] = Xbar
        filtered_state_covariances[t] = P0
    elif t != 0:
        filtered_state_means[t], filtered_state_covariances[t] = (
        kf.filter_update(
            filtered_state_means[t-1],
            filtered_state_covariances[t-1],
            observation = measurements[:,t])
        )

# plot of the resulting trajectory
plt.figure()
plt.plot(rearwheel_y, 'k.', label = 'Observations')
plt.plot(filtered_state_means[:, 3], "y-", label="Filtered Positions", markersize=1)
plt.ylim(180, 200)
plt.grid()
handles, labels = plt.gca().get_legend_handles_labels()
by_label = dict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys())

Older version below
OK, so if you change your Kalman filter to something like this:
kf = KalmanFilter(transition_matrices = np.array([[0.1,0.9], [0, 1]]),
                  observation_matrices = np.array([[1., 0]]),
                  transition_covariance = np.array([[1, 0],[0,1]]),
                  observation_covariance = np.array([[100]]),
                  initial_state_mean = np.array([[0,0]]),
                  initial_state_covariance = np.array([[0.13, 0], [0, 0.13]]))

you get:

You can change the observation_covariance up or down to get smoother or more jagged tracks. The smoother tracks will have more delay.

Original post
To get anything more sensible out of it, you'll have to choose a better model.  The model you chose before your recent edit had that the measurement changed as a simple random walk.  If that's all you know about how the measurement changes, then what you have is as good as it gets.
If you know more, for example if the measurement has some inertia, then you will need to choose different parameters, especially for your transition matrix.
Without understanding a bit more about what your measurements are, I won't be able to help much. Have a look at this answer and see if that helps.  That answer looks at including position and velocity into the state, so that the transition_matrices parameter is a $2\times 2$ matrix.
I couldn't load your data sensibly, so I just generated an offset sinusoid and changed the KF parameters.  This can yield something like the plot below.  However, I'd need to know more about your signal model to get it to work sensibly for your situation.

Updated Code Below
import numpy as np
import pandas as pd
from pykalman import KalmanFilter
from numpy import ma
import matplotlib.pyplot as plt

# Add time
time = np.linspace(-np.pi, np.pi, 201)

# Locate and clip relevant snippet of motion, removing unecessary data at the beginning and end
#rearwheel_y = rearwheel_y[385: 600]

rearwheel_y = np.sin(20*time) + 10

# The numpy.ma module provides a convenient way to address the issue of dropout, with masked arrays.
# When an element of the mask is False, the corresponding element of the associated array is valid and is said to be unmasked.
# When an element of the mask is True, the corresponding element of the associated array is said to be masked (invalid).
# rearwheel_y = ma.masked_values(rearwheel_y, 0)

# time step
dt = 1
# initial_state_mean
X0 = 0 #rearwheel_y[0]
# initial_state_covariance
P0 = 0
n_timesteps = len(rearwheel_y)
n_dim_state = 2
filtered_state_means = np.zeros((n_timesteps, n_dim_state))
filtered_state_covariances = np.zeros((n_timesteps, n_dim_state, n_dim_state))
# Kalman-Filter initialization
kf = KalmanFilter(transition_matrices = np.array([[0.01,0.99], [0, 1]]),
                  observation_matrices = np.array([[1., 0]]),
                  transition_covariance = np.array([[1, 0],[0,1]]),
                  observation_covariance = np.array([[1000]]),
                  initial_state_mean = np.array([[0,0]]),
                  initial_state_covariance = np.array([[0.3, 0], [0, 0.3]]))
# iterative estimation for each new measurement
for t in range(n_timesteps):
    if t == 0:
        filtered_state_means[t] = X0
        filtered_state_covariances[t] = P0
    elif t != 0:
        filtered_state_means[t], filtered_state_covariances[t] = (
        kf.filter_update(
            filtered_state_means[t-1],
            filtered_state_covariances[t-1],
            observation = rearwheel_y[t])
        )

# plot of the resulting trajectory
plt.figure()
plt.plot(rearwheel_y, 'k.', label = 'Observations')
plt.plot(filtered_state_means[:, 0], "y-", label="Filtered Positions", markersize=1)
plt.ylim(8, 12)
plt.grid()
handles, labels = plt.gca().get_legend_handles_labels()
by_label = dict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys())


Your approach will only denoise anything, if your signal is actual periodic. That's typically the case if you use a sine wave or periodic random noise or sweep as an excitation signal for a measurement. It can also work, if you vibration source is periodic (e.g. an engine running an constant rpm).
In this case $D$ should be chosen to represent EXACTLY one period. $N$ is simply determined by the length of your data set. The longer the better.
This works best if the excitation clock and acquisition clock are phase locked, i.e. derived from the same clock source. If that's not the case, it may be required to do some sort of clock recovery or clock tracking first.


Is the backwards transform just an IFFT (for forward transform)?

Yes.
Compare the "what FFTW actually computes", here


Finding set of orthogonal basis functions for composite signals

Your $S_1,\ldots, S_5$ already are your orthogonal bases, exactly the way your first equation writes them down.
To check that: write down the definition of your inner product, and insert, for example, $S_a$ and $S_b$. You'll see that $\left<S_a,S_b\right>=0 \iff a\ne b$.
This follows quite directly from the linearity that a (real) inner product needs to have:
\begin{align}
\Big<\sum_i A_{a,i}\sin(2\pi f_{a,i}t -\phi_{a,i}),
\sum_k A_{b,k}\sin(2\pi f_{b,k}t -\phi_{b,k})
\Big>&\\
= \sum_i
\Big<
A_{a,i}\sin(2\pi f_{a,i}t -\phi_{a,i}),
\sum_k A_{b,k}\sin(2\pi f_{b,k}t -\phi_{b,k})
\Big>&\\
=\sum_i\sum_k
\Big<
A_{a,i}\sin(2\pi f_{a,i}t -\phi_{a,i}),
 A_{b,k}\sin(2\pi f_{b,k}t -\phi_{b,k})
\Big>
\end{align}
and sines of different frequencies are orthogonal.


Convolution for non-signal-processing background

I disagree with the assertion in the answer by user7090 that the calculation in question is not a convolution; it is indeed a convolution.  Let us consider the case $D=2$. Ignoring $Z$ for now,
the five terms in that sum
$\displaystyle \sum_{k = i-2}^{i+2} x(t_k)e^{-|i-k|/\Delta}$ on the left in the OP's question are
$$\begin{align}
x(t_{i-2})e^{-|i - (i-2))|/\Delta} &= x(t_{i-2})e^{-2/\Delta}\\
x(t_{i-1})e^{-|i - (i-1))|/\Delta} &= x(t_{i-1})e^{-1/\Delta}\\
x(t_{i-0})e^{-|i - (i-0))|/\Delta} &= x(t_i)e^{-0/\Delta} = x(t_i)\\
x(t_{i+1})e^{-|i - (i+1))|/\Delta} &= x(t_{i+1})e^{-1/\Delta}\\
x(t_{i+2})e^{-|i - (i-2))|/\Delta} &= x(t_{i+2})e^{-2/\Delta}
\end{align}$$
and so the sum is very much a convolution of the  sequence $x$
whose $i$-th term is $x(t_i)$
with a length-$5$ sequence $h$ whose $n$-th term is $e^{-|n|/\Delta}$
for $|n| \leq 2$ in this case (and a length-$(2D+1)$ sequence
 with $|n| \leq D$ more generally).

What you, the OP, need to understand is the following (and don't
get bogged down in words like convolution).  You are doing a
smoothing of the data that you have collected.  Each $x(t_i)$ of the
original data is being replaced by a $\tilde{x}(t_i)$ which is formed
as a weighted sum of $x(t_i)$ and its immediate neighbors.  The
weights assigned to $x(t_i)$ is the largest, the weight assigned to
its nearest neighbors $x(t_{i+1})$ and $x(t_{i-1})$ are smaller (but equal),
the weights assigned to slightly more distant neighbors $x(t_{i+2})$ and $x(t_{i-2})$ is yet smaller, and so on, till the we get to $x(t_{i+D})$ and $x(t_{i-D})$ which get the smallest weights of all. Anything more
than $D$ away gets a weight of $0$.  The weights in question are
normalized in the sense that the total weight is $1$, but let
us look at the raw weights which are alleged to be those $e^{-|i-k|/\Delta}$.
What you really have (before normalization of the weights) is
that sum (let us use $D = 3$ for simplicity) which is

$x(t_i) + \left.\left.e^{-1/\Delta}\right[x(t_{i+1}) + x(t_{i-1)}\right]
+ \left.\left.e^{-2/\Delta}\right[x(t_{i+2}) + x(t_{i-2})\right]
+ \left.\left.e^{-3/\Delta}\right[x(t_{i+3}) + x(t_{i-3})\right]$

$= e^{-3/\Delta}x(t_{i-3}) + e^{-2/\Delta}x(t_{i-2}) +
e^{-1/\Delta}x(t_{i-1}) + x(t_i) + e^{-1/\Delta}x(t_{i+1}) + e^{-2/\Delta}x(t_{i+2})
+ e^{-3/\Delta}x(t_{i+3})\quad (1)$
where on the right side we have the elements occurring in the
natural order, flanking the central element $x(t_i)$ on both sides.
Now, the sum of the raw weights exhibited in $(1)$ is
$$e^{-3/\Delta} + e^{-2/\Delta} +
e^{-1/\Delta} + 1 + e^{-1/\Delta} + e^{-2/\Delta}
+ e^{-3/\Delta} \tag{2}$$
which is just the element $Z$ in the question, and so the
normalized weights are just the ones exhibited in $(2)$
divided by $Z$. A better way of writing $Z$ might be
$$Z = \sum_{k=-D}^D e^{-|k|/\Delta} = 1 + 2\sum_{k=1}^De^{-k/\Delta}
\tag{3}$$
and perhaps a better way of expressing the first sum might be
$$\tilde{x}(t_i) = \frac{1}{Z}\sum_{k=-D}^D x(t_{i-k})e^{-|k|/\Delta}$$
which expresses more clearly that we are forming a weighted sum
of $x(t_i)$ and its nearest neighbors and that the weights assigned
to neighbors at distance $k$ decrease exponentially as a function
of $k$ (here $k \geq 0$ of course since we are not distinguishing
between the directions of separation).


Are all images 2d or also we have 3d images?

The term image derive from  Latin  imago. I fully agree with wikipedia saying that it  denotes an artifact

  that depicts visual perception [...] that resembles a subject—usually
  a physical object—and thus provides a depiction of it.

As such, it is a representation of some reality, combining spatial "ordinal" cooordinates (with some natural order, in 2D and 3D for us humans) and  "cardinal" coordinates (luminance/chrominance values, spectra, categories).

The dimension attributed to images to signal/processing people is generally the one of the "ordinal" cooordinates (2D or 3D). To other practitioners who plot 2D images as 3D plots (using elevation for scalar fields in the "cardinal" dimension, non-image-processing people often call 2D image "tridimensional" (bottom-right "Arnold Schwarzenegger" in Terminator), because they look at this artifact in our natural 3D world:

.

Here, $x$ and $y$ are the ordinal dimensions, $f(x,y)$ makes the third. This visualization of images in the resulting topological/topographic map is exploited for instance in the field called "mathematical morphology",  where surfaces and regions are interpreted as geographical objects: valleys, mountains, watersheds, erosions, etc.

If you add (a 1D) time (a very ordinal dimension so far) for motion, then you can add another dimension: standard 2D+t video, 3D video, etc.

If we call for scale-space theory, where one bluntly considers that "objects are meaningful over a certain range of scales only", one can smooth a given image with smoothing kernels to get different images:

and pile them up to yield an additional "scale" dimension. Then, you can extend the discussion to images at arbitrary coordinates, like point clouds, with motion, etc.

Finally, I believe that image data (and their dimensions) can be considered as at least three objects:

how they should be in reality (if such a thing exists),
how they are modeled or described,
how we process them.

and their inherent dimension depends on that. Remember however that one can map two continuous spaces of different dimensions, albeit in a non continuous manner. For discrete spaces, like the ones for classical image/volume coordinates, mappings exist, like zigzag of Hilbert/Peano curves for image compression (JPEG, J2K), so dimensions are a matter of data structuration.


extracting ultra accurate phase and frequency with FFT with real world data

You only have 64 samples per FFT. That means your FFT bins are spaced f_sample/64; if your sample rate is 128·50 Hz, then your frequency resolution, without interpolation, is just 100 Hz. To achieve a frequency resolution of $0.01\,\text{Hz}=\frac{50\,\text{Hz}}{50\cdot 100} = \frac{\frac{f_{sample}}{128}}{5000}= \frac{f_{sample}}{640\cdot5000}$, you'd need an FFT of length 3.2 millions. You can't do that on an M4, obviously, so: the FFT is not a tool applicable to your problem here.

Two approaches:

Either, you first reduce the sampling rate (including sufficient analog anti-aliasing filters) or you decimate your signal (including sufficient digital anti-aliasing filters) to a rate that describes a much smaller bandwidth of potential frequencies that you actually care about. Afterwards, do the FFT. Notice that this always always trades length of observation for precision: it's really that simple, you can't get 1/100 Hz of resolution without 100 s of observation. If things change within that period: bad luck.

Or, you don't do that and instead use parametric estimators for frequency and phase. I kind of feel I recommend that here. Parametric estimators can be as simple as zero-crossing detectors, or as complex as Eigenvalue-based estimators, and what you use depends on a lot of factors that you don't discuss – it's hard to guess what you need.

However, I like to point people at Peter K.'s A Review of the Frequency Estimation and Tracking Problems. Many of these algorithms will deliver phase info.

Regarding an accuracy requirement of 0.1° or 0.01 Hz in general: You can't have arbitrary precision with limited observation if you have any noise in your signal. Estimation theory doesn't allow for that. I think you'll be running into theoretical limits rather fast, depending on which algorithm you end up (and using the lowest-variance method might not always be the best choice!).


Finding the Wiener filter transfer function

What you've listed is a classical set of linear equations known as the Wiener-Hopf Equations. There are many methods to solve this set of equations. These equations generalize the Yule-Walker equations for Autoregressive (AR) modeling. In fact, if the Wiener-Hopf equations are used to solve for a linear predictor of your signal, the equations are identical to the Yule-Walker equations.

So in reality, there are many ways to solve these equations depending on how you're modeling your signal. This is where the concepts of Wiener filtering, Wiener smoothing, Wiener deconvolution arise. In each of these, they can be segmented further into causal/anti-causal realizations, FIR/IIR, etc. I recommend searching for the Wiener-Hopf equations and understand what they mean, then solving them comes secondary.

As an example, consider if you assume that you have a signal with additive noise that is uncorrelated with the signal. Then the cross-correlation terms have no dependence on the signal. Then you need to model your noise and plug in the resulting autocorrelation sequence on the right hand side of the equation.

Cheers


This is termed 3D object retrieval. There are many possible approaches to follow. First, if your meshes are nice, clean and descriptive, you could then consult the academic literature of computer graphics, while computer vision deals with more noisy scenes, partial visibility, non-meshed data. As I am a computer vision oriented person, I would be biased towards such approaches. But keep in mind that Eurographics organizes a workshop on 3d object retrieval every year, called 3DOR. Each year different problems are addressed, but it is very common to find something meeting your requirements.
The following might be good references in the direction you would like to go:

From Low-Cost Depth Sensors to CAD: Cross-Domain 3D Shape Retrieval via Regression Tree Fields

Database-Assisted Object Retrieval for Real-Time 3D Reconstruction

Feature-based Similarity Search in 3D Object Databases

There are many many more works though. For example, here are some, incorporating deep learning:

3DNet: Large-Scale Object Class Recognition from CAD Models

3D Shape Retrieval Using a Single Depth Image from Low-cost Sensors

3D ShapeNets: A Deep Representation for Volumetric Shapes


No, it is not possible. A proof follows.

Let $\mathrm{E}(L \mid 0)$ be the expected codeword length for the prefix code $0$ that minimizes it. The answer to the question is "no" if it is always true that:

$$\mathrm{E}(L \mid 0) \le \mathrm{E}(L),\tag{1}$$

where $\mathrm{E}(L)$ is the total expected codeword length using a randomly chosen prefix code $i$ from a set of $N$ prefix codes denoted by $i\in\{1, 2, 3, \ldots, N\}$. We can divide both sides of Eq. 1 by $\mathrm{E}(L \mid 0)$, which is positive, so we could equivalently show that:

$$\Leftrightarrow 1 \le \frac{\mathrm{E}(L)}{\mathrm{E}(L \mid 0)}.\tag{2}$$

Using the law of total expectation:

$$\frac{\mathrm{E}(L)}{\mathrm{E}(L \mid 0)} = \frac{\sum_{i=1}^{N}\mathrm{E}(L \mid i)P(i)}{\mathrm{E}(L \mid 0)} = \sum_{i=1}^{N}\frac{\mathrm{E}(L \mid i)}{\mathrm{E}(L \mid 0)}P(i),\tag{3}$$

where $\mathrm{E}(L \mid i)$ is the expected codeword length for a prefix code $i$ with an associated non-zero probability $P(i)$. Because the prefix code $0$ minimizes the expected codeword length, we have $\mathrm{E}(L \mid 0) \le \mathrm{E}(L \mid i)$ $\Rightarrow$ $\mathrm{E}(L \mid i)/\mathrm{E}(L \mid 0) \ge 1$ $\Rightarrow$ $\mathrm{E}(L \mid i)/\mathrm{E}(L \mid 0) - 1 \ge 0$ for all $i$. It follows that:

$$\frac{\mathrm{E}(L)}{\mathrm{E}(L \mid 0)} = \sum_{i=1}^N P(i) + \sum_{i=1}^{N}\left(\frac{\mathrm{E}(L \mid i)}{\mathrm{E}(L \mid 0)}-1\right)P(i) \ge \sum_{i=1}^N P(i)= 1,\tag{4}$$

where the second sum is never negative. This proves Eq. 2. Therefore, the answer to the question is "no".


Mathematical Approach to Detect If a 2D Signal Is Separable

Let's assume our data is in finite dimension.
So $ x \left[ m, n \right] \in \mathbb{R}^{M \times N} $. So it can be written as a matrix $ X \in \mathbb{R}^{M \times N} $.

Using the SVD Decomposition the matrix can be written as:

$$ X = \sum_{i = 1}^{n} {\sigma}_{i} {u}_{i} {v}_{i}^{T} $$

Seprable 2D Matrix is defined as a rank 1 Matrix which can be composed by Outer Product of 2 vectors.
Looking at the SVD Decomposition of $ X $ we can conclude that $ X $ is separable operator if and only if $ \forall i > 1 \; {\sigma}_{i} = 0 $ and it is given by:

$$ X = {\sigma}_{1} {u}_{1} {v}_{1}^{T} $$

So, here is a simple recipe:

Write you data in Matrix Form.
Apply the SVD on the matrix.
In MATLAB it can be done using svd().
Look for the rank of the matrix. Namely how many non vanishing singular values it has.
If the number is 1 the matrix is considered separable. It can be separable by the its first left eigen vector and first right eigen vector.

For example on Separable Operators have a look at - How to Prove a 2D Filter Is Separable?.


Successive dynamic time warping with a set of window sizes

Hmm

There are few properties to help you.

Let Wi be a warping window of width i

It is the case that DTW_wi(A,B) >= DTW_wi+1(A,B)

Other than that, there is little to say.

If you have not already, you should read
http://www.cs.unm.edu/~mueen/DTW.pdf
and
http://www.cs.ucr.edu/~eamonn/CIKM_2016_DTW_clustering.pdf

If you are only interested in the smallest items in the matrix, then using lower bounds (the LB_Keogh) will give you an easy 100-fold speed up.