Mitigate floating-point numerical errors for very-low corner frequency low-pass filter with DSP

I think your issue might be coefficient quantization and filter topology. A direct form biquad has poor quantization effects about 0 and π radians. It’s easier to analyze such effects in fixed point, but even though floating point has a much larger range, it still has shortcomings. In particular, if you add a very small number to a very large number, the small number disappears because it can’t be aligned for the operation in the available number of mantissa bits. This can cause the order of operations to affect the result. For instance, where S is a small number and L is large, L - L + S = S, but L + S - L = 0.
Udo Zolzer covers the differences between several filter structures in his book Digital Audio Signal Processing. I borrow direct form quantization effects on pole locations from the book:

See how precision is lost near 0 and π. Other filter topologies might have higher precision near 0, while being much worse near π, which may be a good tradeoff for uses such as yours. The Gold and Rader form has a very even distribution, it looks like a perfect grid.
Another simple and popular filter that has good quantization characteristics at low frequencies is the "Chamberlin" state variable filter. There are improved versions of this filter, as it has problems at higher frequencies (from about one-sixth of the sample rate and up), but the plain Chamberlin is very good at low frequencies, where you need it.
See my article on the Chamberlin state variable filter here:
The digital state variable filter
Zolzer presents modified Chamberlin structures here:
The Modified Chamberlin and Zölzer Filter Structures
In particular, see the graph of quantization effect near zero for the Chamberlin structure—very dense near zero, at the expense of poorer performance at high frequencies, compared to the direct form graph:


Can't implement correctly Python IIR and FIR Filters

Right now I have an error that says: Error: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

    import myDSP
        import matplotlib.pyplot as plt
        import scipy.signal as signal
        import numpy as np

        fs, x = myDSP.readWav('Ring05c.wav')

        fc1=1200
        fc2=3000
        fc= np.arange(fc1,fc2+1)
        ft=150
        rp=0.1
        rs=40
        t = np.linspace(-1, 1, 201)
        #x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1)+0.1*np.sin(2*np.pi*1.25*t + 1) +0.18*np.cos(2*np.pi*3.85*t))
wp=np.array([fc-ft/2])/(fs/2)
ws=np.array([fc+ft/2])/(fs/2)
        L, Wn = signal.cheb1ord(wp,ws, rp, rs);
        b, a = signal.cheby1(L, rp, Wn,'low')
        f,H=signal.freqz(b,a,worN=256,plot=None,fs=fs)
        plt.plot(b,'.-b',label='b coefficients')
        plt.plot(a,'.-r',label='a coefficients')
        plt.legend()

        #h = signal.lfilter(b, a,x,axis=0)
        #w=signal.hann(65)


Connection between filter equation in frequency domain and difference equation in time domain?

Hint:

The transfer function $H(\omega)$ is equal to the ratio of output and input spectra:

$$H(\omega)=\frac{Y(\omega)}{X(\omega)}$$

So your first equation can be rewritten as

$$Y(\omega)=X(\omega) - 1.176X(\omega)e^{-j\omega} + X(\omega)e^{-2j\omega}$$

What do you get if you transform this equation to the time-domain?


The description complex helical sequence indicates that it is a sequence of complex numbers that, when plotted in three dimensions (real part, imaginary part, position in the sequence), resembles a screw, or helix.

See this illustration (by RobHar, on Wikimedia Commons):

From a mathematical point of view, all signals can be decomposed into complex sinusoids (sine waves), which are terms of the form $A(x) = e^{(\omega x +\phi)i}$. This is an often useful perspective. As $x$ is varied, $A(x)$ describes a circle in the complex plane; plotted in three dimensions, it is the helix from above.

Real-life signals normally consist of real numbers only, i.e. the imaginary part is always zero. This means that the complex sinusoids of which it is composed come in mirror pairs, where one runs clockwise and one anticlockwise; otherwise (frequency, amplitude, phase), they are identical. In such a pair, the imaginary parts cancel, and the result is a real-only signal.

An analytic signal is made from a real signal by removing one half from each of those pairs so that all remaining sinusoids move in the same direction. It is no longer real-valued, i.e. the imaginary part is in general non-zero, but the real part is identical to the original signal.

(Technically, the sense of rotation is considered the sign of the frequency $\omega$, and an analytic signal only contains sinusoids with non-negative frequency.)


Adaptive Predictive Coding: Transmit signal AND prediction coefficients

This sounds like a variant of Differential Pulse-Code Modulation. In this set of techniques, the correlation of adjacent values and the effect of quantisation are exploited for "compression" (whether this is for transmission or storage). That is, instead of encoding the sample values, you encode the difference between them which, within limits, could be described by a much smaller number than the values themselves.

However, another technique that might be useful to you here is run length encoding. In run length encoding, the encoded word describes either the sample value itself or how many repetitions of the same value are to be decoded.

In either case, you would require higher bandwidth around the signal transients (many messages sent between the encoder and the decoder) and much lower during the longer times of constant values.

Hope this helps.


Here's the outline of the argument, feel free to fill in the details.

The averaging operator is like a convolution with a "square" pulse of height $1/T_s$ supported on the interval $[-T_s/2, T_s/2]$.

You can express the square pulse as a sum of two heaviside step functions.

Finally, recall the Laplace transform of a step function $\mathcal{L} \{H(t-T_s/2)\}(s) = \frac{e^{-sT_s/2}}{s}.$


How does `envelope` use the sliding window to determine RMS in MATLAB?

Since your window is not centered, you receive a shift of half-length of the chosen window.


First, we define $h_3[n] = h_2[n] \star h_2[n] = [1, 2, 1]$. From this result, we know that $h_1[n]$ must have 5 elements (so that $h[n]$ ends up with 7 elements).

Let's define $h_1[n] = [g_1, g_2, g_3, g_4, g_5]$. We can find these as follows, using the definition of discrete convolution.

First, we know that $g_1 h_3[0] = h[0] = 1$, so $g_1 = 1$.

Then, the next value of the convolution is $h[1] =  g_2 h_3[0] + g_1 h_3[1] = 5$, from which we deduce that $g_2 = 3$.

Then, $h[2] = g_3 h_3[0] + g_2 h_3[1] + g_1 h_3[2] = 9$, or $g_3 = 2$.

You can stop here, since (by symmetry, in this particular example) $g_4 = g_2$ and $g_5 = g_1$.

In general, the problem can be seen as solving a system of linear equations. In this example, we need to solve the system:

$$\begin{bmatrix}
h_3[0] & 0 & 0 & 0 & 0 \\
h_3[1] & h_3[0] & 0 & 0 & 0 \\
h_3[2] & h_3[1] & h_3[0] & 0 & 0 \\
0 & h_3[2] & h_3[1] & h_3[0] & 0 \\
0 & 0 & h_3[2] & h_3[1] & h_3[0] \\
0 & 0 & 0 & h_3[2] & h_3[1] \\
0 & 0 & 0 & 0 & h_3[2]
\end{bmatrix}
\begin{bmatrix}
g_1 \\ g_2 \\ g_3 \\ g_4 \\ g_5
\end{bmatrix}
=
\begin{bmatrix}
h[0] \\ h[1] \\ h[2] \\ h[3] \\ h[4]  \\ h[5] \\ h[6]
\end{bmatrix}$$

Thanks to @Hilmar for pointing out a boneheaded mistake in my original answer!


How much of a problem is the Gibbs phenomenon and would this solve it?

I would say that this is interesting.  There has been a lot of work done regarding the study of Gibbs Phenomenon.  You should check out the following document to get a better understanding of how it comes up in practical DSP applications:

http://people.clarkson.edu/~ajerri/books/examples/Gibbs_Book.pdf

The typical way to manage Gibbs Phenomenon is to use time domain window functions that taper data at the start and finish.  The window functions reduce spurious contributions to frequency domain information that come from discontinuities at the edges of data sequences.

I haven't seen much application of generating signals by composing them with individual sine waves.  Generally signal generation is done directly in the time domain.  I'm not sure how the function constructions you've documented can be employed to solve a practical problem, but perhaps there is application if you can identify an appropriate problem.


What is the suitable mathematical form of the multipath channel impulse response?

Let's say we have a function that we transmit which is:
$$x(t)$$
Transmitting the function is kinda easy part. We assume there is just the function we are trying to send.

However things are a little bit ugly in the air. Let's look at the receiver part.

What we generally consider first is the AWGN. Which is:
$$r(t) = y(t) + n(t)$$

$r(t)$ is our received function. $y(t)$ is our function that reaches our hand and we assume we don't know what it is. In this case, for now, please consider $y(t) = x(t)$. Lastly we have gaussian distributed noise which is: $n(t)$.

Well, we have some fading types should be added. Let's combine it together and call them Channel Response. Channel responses should vary over time but for this situation I will assume it is constant.

Let's call this $h_k$.

The subscript of $h$ basically means we have different channel response for each channels. And hence we have:
$$h*y(t) + n(t)$$

But, hey! I don't see that at the receiver! There is some odd peaks at the signal! What is that thing?

Sadly, this is multipath and, at the receiver, you are staring at this signal:

$$(\sum\limits_{k=1}^K = h_k*y(t)) + n(t)$$

Aaaand no! They are not coming at the same time. Let's add a delay $\tau_k$

$$(\sum\limits_{k=1}^K = h_k*y(t - \tau_k)) + n(t)$$

$K$ is the number of multipath. Please notice that as the multipath occurs, each path has it's own channel. As you can see, there are different channel responses and delay times for each.

I hope it will be useful.

As always, have a nice day.


What are characteristics of different MR images?

Speaking from a medical imaging perspective, the SNR of MR images mainly depends on signal acquisition schemes, image reconstruction parameters, and the MR scanner hardware, such as, number of receive coils. To know about more about other parameters go here: http://wiki.ci.uchicago.edu/pub/HNL/DifficultQuestion/howdoesmriwork.pdf.

Parallel MR imaging is one of the acquisition technique which improves SNR. For example, I'm sharing with you the results from a paper published in journal of magnetic resonance imaging on "Measurement of Signal-to-Noise Ratios in MR
Images" by Olaf Dietrich and others below:

PS: Take the reciprocal to get noise to signal ratio.


Unit impulse response $h[n]$ of a discrete time system with multiple characteristic roots 0

If you have a difference equation
$$\sum_{k=0}^Na_ky[n-k]=\sum_{k=0}^Mb_kx[n-k]\tag{1}$$
and if $M>N$, the general expression for the corresponding causal impulse response is
$$h[n]=A_0\delta[n]+A_1\delta[n-1]+\ldots+A_{M-N}\delta[n-(M-N)]+u[n]\sum_{k}c_k\gamma_k^n\tag{2}$$
where $\gamma_k$ are the roots of the characteristic equation, and where I've assumed simple roots. Note that roots $\gamma_k=0$ do not contribute to $(2)$ (even if there are multiple roots equal to zero).
For the given example we have $M-N=3$ and one non-zero root $\gamma=\frac12$. Hence the solution has the form
$$h[n]=A_0\delta[n]+A_1\delta[n-1]+A_2\delta[n-2]+A_3\delta[n-3]+c\left(\frac12\right)^n\tag{3}$$
The constants $A_i$, $i=0,1,\ldots,3$, and $c$ must be determined by iterating through the difference equation with $x[n]=\delta[n]$, $h[n]=y[n]$, and $h[-1]=0$.
If you do that you'll find that $c=0$, which means that the given difference equation describes an FIR system. In terms of transfer functions this means that there is a pole-zero cancellation, leaving only poles at the origin of the complex plane.


Reason for >8 channels in an audio codec

Since it will most likely not hurt to have the possibility for more streams, this might for one be just to future-proof the system. Nothing would be more stupid than to specify for 8 channels and suddenly 9.1-systems become the next big thing and your codec is not suitable.

Other than that I can imagine bigger scenarios with much more channels such as recording all the voices in an orchestra. But I honestly don't know whether or not AAC is of any importance there.


Answer
If you need to guarantee operation within these ranges, then yes. (Assuming the 0.001 Hz is a typo and meant 0.1 Hz.)
Stochastic Modelling
That "1ppm" is typically a "condensed" number, e.g. the observed maximum over 24h of oscillations, or can be found as integral of phase noise power, which is often given in an oscillator's datasheet. Typically, small frequency errors are more likely than large ones.
It might make sense to model your system probabilistically: Your frequency errors in receiver and transmitter, $X_{RX}$ and $X_{TX}$ are independent and will have a probability density function $f_{X_{RX}}(x)$ and $f_{X_{TX}}(x)$, respectively.
So, the probability density function for the difference $\Delta=X_{TX}-X_{RX}$ is given by $f_\Delta(\delta) = f_{X_{TX}}(x)*f_{X_{TX}}(-x)$, $*$ being the convolution operator.
This will allow you to define something like "in x % of cases, it should work, so I only need to account for a range of at most this and that in frequency error".
For example: This makes a lot of sense, since frequency synchronization is not the only thing that can fail in a receiver; if your channel has a 0.1% chance of simply being too noisy, then building a PLL loop filter that is wide enough to work for 99.99996% of cases doesn't make much sense, and you should probably rather reduce your loop filter bandwidth to make the receiver's phase recovery less sensitive to noise.


DFT and then IDFT does not provide the same signal

If you want your IFFT result to match the input to the FFT, don't throw away either the "imaginary" (or sine correlation) components, or the "upper half" of the spectrum, before the IFFT.  Also make sure the input to the IFFT is exactly conjugate symmetric, so that all the imaginary components in the IFFT results get cancelled out (except for rounding errors or quantization noise).  The length of the IFFT has to be the same as that of the FFT.

Your 99 and 95 Hz components are also the -1 and -5 Hz components needed (in complex conjugate symmetry) to cancel out any rotation of the time domain data into the imaginary plane.  Since the FFT is a complex to complex transform, you can't just ignore these components even when dealing with strictly real data.

The 0 Hz or DC component of an FFT result appears double height because it includes it's own mirrored value.

Unless your sinewave's are exactly periodic within the FFT length, you will need to interpolate the true amplitudes, since their peaks will be between the FFT result bin centers.

The IDFT and IFFT can be used as signal generators, although great care must be taken with phase for any frequencies that are not exactly periodic in the FFTs length.  See Phase Vocoder analysis/synthesis for details.


Measure the length of curves via a smartphone's accelerometer

Accelerometers and gyroscopes provide corrupted measures of the acceleration $a$ and angular velocity $\omega$. So, they need to be preprocessed in order to use them.
Firstly, measure provided by the accelerometer ($a'$) is affected by Earth gravity $g$ (in the local frame reference), a bias $b_g$ and some random noise $\eta_a$. So if the true acceleration is called $a$, one has that measure is:
$$ a' = a + g + b_a + \eta_a$$
So in order to integrate accelerometer measures, one should know at least $b_a$ and $g$ in the accelerometer frame, and subtract them from $a'$. Moreover, in some datasheets you can find covariance of $\eta_a$. However, bias is time dependent term and it needs to be estimated each time we initialize the system, as well as the gravity in the local frame. In some cases, if the movement is not in gravity direction (2D movement) and gravity is aligned to some accelerometer axis (lets say $z$ axis), one can neglect $a_z$ and the knowledge of $g$ is not required (maybe it is your case).
Secondly, gyroscope is also corrupted by bias $b_g$ and noise $\eta_g$. So, if the true angular speed is $\omega$, one has that the measure is:
$$ \omega' = \omega +b_g + \eta_g$$
Again, you need to know at least the gyro bias. and subtract it from the measure.
Everything becomes easier if your device has no bias or it can be neglected w.r.t $a$ and $\omega$, and your movement is 2D (gravity is always oriented in the same way).
So once that you have your corrected measures, you can integrate them to obtain the trajectory:
$$ R(t+\Delta t) = R(t) \exp(\omega \Delta t)$$
$$ p(t+\Delta t) = p(t) +v(0) t + \frac{1}{2}R(t)a(t)t²$$
Where $R$ is the rotation matrix (see rotation matrix) and $\exp$ is the exponential map (see exponential map). $v(0)$ is the initial speed.
It can be concluded that having a known initial speed, a 2D movement, and a low corrupted gyro/accelerometer, it can be easily estimated the trajectory.
Edit: Calibration (i.e. getting $g$, $b_a$ and $b_g$) can be intermediately done if device remains static. In this case, accelerometer gives $a' = g + b_a +\eta_a$ and gyro gives $\omega' = b_g + \eta_g$


Denoising Oscilloscope data the right way [python3]

If you are simply interested in plotting the data then any data reduction technique would do,even if it appears to be crude.

Effectively, the plotting function itself will not plot all the data, because the space assigned to the plot has a finite number of $N_x \times N_y$ pixels assigned to it. For example, if your plotting area was $1024 \times 768$, then any waveform with more than 1024 samples would need to somehow be re-scaled or re-sampled. Plotting functions do that internally anyway but always retain a link to the original data to do the whole process again, usually onUpdate.

What you could try to do is to resample your waveform and specifically downsample it.

Downsampling is composed of two steps: The step of anti-alias filtering and the step of decimation.

If you are purely interested in plotting a waveform, you can try to simply decimate it with a low factor (say, retaining 1 sample out of every 5) and then send it for plotting. It will still have more data than required by the plot but it will be of more manageable size.

The reason for opting for pure decimation here is because the filtering step (although necessary) can be costly (in terms of computation) and take a bit longer. If you can afford the filtering time however, I would suggest that you do a proper downsampling.

If you are using Python, you can try this function (or this one) to handle the downsampling step.

It goes without saying here that if you wish to process the signal further, you must downsample it properly.

Finally, another reason why you are running out of memory might be the fact that you are retaining all data in memory even when you are done using it. Therefore, once you have sent a waveform for plotting, try del()eting it, possibly followed by gc.collect() to ensure that that part of memory will become immediately available.

Hope this helps.


Why does a synchrosqueezed wavelet transform show oscillating behavior?

This was interesting to figure out. The key lies in the phase transform, and how CWT interacts with own derivative upon insufficient component separation. Relevant are, and I'll be answering, the following:

What causes the wavy pattern?
Is the wavy pattern truly a sine, or a lookalike?
Why is only one band wavy, not both?
What can be done about it?

Answer split in two to separate meat from supplement, but much of insight's in latter; sections labeled in recommended reading order.
TL;DR: skip to Conclusion. Answer code.

0. How does the phase transform work?
A helpful requisite for following sections, covered here.

1. What causes the wavy pattern?
Begin by visualizing the phase transform, w:

The lower frequency's CWT and SSQ are nearly that of a pure tone, but the higher frequency is all wiggly. Let us plot one such wiggly row:

This sure looks sinusoidal. But is it? Stay tuned.
First: is the wavy SSQ explained by the w? Absolutely: w is wavy across timesteps, meaning for every column, we reassign CWT rows to a different SSQ row, in a wavy/sinusoidal pattern (see [complement] on reassignment). Our answer thus lies in understanding why w is wavy to begin with. Second answer begins here.

4. Inseparability: pre-transform perspective
Recall, the CWT and dCWT operations both take place in the frequency domain; we reach the lowest level of analysis by inspecting it directly. With xh = fft(pad(x)), psih = freq-domain wavelet, and dpsih = freq-domain wavelet with the derivative (* 1j * xi / dt), we plot both, overlapped and rescaled to fit in same graph, as well as the result of the multiplication - for our same row, and one on higher scale:

First, note for our original row scales[80] = 5.08, we indeed have a sum of two pure cisoids (analytic wavelet -> all negative frequencies are zero, and Morlet has pure-real DFT; cosine input -> pure real DFT; thus a DFT coefficient multiplies its complex exponential basis without cancellation), one at higher frequency having much greater amplitude.
Next, compare the cisoid amplitudes for psih vs dpsih; higher frequency dominates more for dpsih than for psih for both example scales. Coincidence? Well, they both do * xh, thus we expect the dpsih wavelet itself to skew toward right. Plotting,

this is confirmed. Does this happen for every scale? Yes: * xi = * linspace(0, pi, N) - taking the derivative involves multiplying by a line, thus skewing the entire wavelet toward higher frequencies, also breaking even symmetry (though not a lot).
All puzzle pieces are in. The skewness creates a difference in relative amplitudes at the two frequencies w.r.t. the non-derivative. The general form for dWx / Wx is then (still assuming purely real fft(Wx)):
$$
\frac{dW_x}{W_x} = \frac{a\exp{(j(\phi_1 + \pi/2))} + b\exp{(j(\phi_2 + \pi/2))}}{c\exp{(j\phi_1)} + d\exp{(j\phi_2)}}, \tag{4} \\
\phi_i = \omega_i \frac{n}{N}, \ \ n=[0, 1, ..., N - 1]
$$
($+\pi/2$ per the derivative). Let's first inspect interesting cases numerically, then derive an exact relation as proof (second answer).

8. What can be done about it?
Strictly speaking, nothing, unless the wavelet decays completely to zero eventually (untrue for many wavelets); as long as there's more than one frequency in the entire DFT of the input, there will be some form of interference. However, practically, such interference can be made negligible by:

Ensuring bands are sufficiently separated
Reducing wavelet frequency width (for Morlet, increase mu)

Let's inspect Morlet mu=5 vs mu=20:

mu=20 is clearly narrower; to get a better idea, single out a few rows:

Overlapped at same center frequency (~same row):

We thus predict mu=20 will show no wiggliness we can see. In fact, mu=20 is overkill, we'll settle for mu=10:

Can we also guess how to make both bands wiggly for small mu? Note, increasing scales lowers the wavelet's center frequency (shifts it left), but also reduces its frequency width; thus our wiggle odds improve toward lower scales (or higher frequencies). This turned out easier said than done; after some experimentation, I've had limited success, but maybe it can be done better:

9. Conclusion

What causes the wavy pattern? The frequency-domain wavelet and its derivative intercept the input bands at slightly different ratios per latter being skewed. The resulting phase transform is (imag. part of) a ratio of sums of complex sinusoids, which is of form $(\cos(\phi) + k_1)/(\cos(\phi) + k_2)$. The reassignment pattern dictated by w can then have a sine-like pattern across rows, which manifests as a sine-looking synchrosqueezed CWT.
Is the wavy pattern truly a sine, or a lookalike? Latter; it can never be exactly sine, but close to. The pattern is characterized by $(\cos(f_1 - f_2) + k_1)/(\cos(f_1 - f_2) + k_2)$, and has a period of $|f_1 - f_2|^{-1}$, where $f_1,f_2=$ input frequencies.
Why is only one band wavy, not both? With increasing scale (and decreasing center frequency), the frequency-domain wavelet moves to the left while shrinking in width; the resulting cisoid proportions work out such that dominantly sine reassignment pattern is toward the higher input frequency, where the wavelet is wider, but the lower one isn't devoid of such a pattern either.
What can be done about it? Separate input frequencies more, or reduce wavelet frequency width across all scales (for Morlet wavelet, raise mu).


This is a moving difference filter which acts as a discrete time derivative and is less sensitive to high frequency noise. When all the samples in the waveform are the same (DC) the last ten samples will cancel out the first 10 samples for DC cancellation but when a change occurs within a 10 sample interval compared to 10 samples prior the difference will be maximized.
This filter is identical to a 10 sample moving average followed by a difference over 10 samples and provides a filtered estimate of the time derivative of the waveform as shown in the diagram below ($z^{-10}$ in the diagram refers to a 10 sample delay and MAF refers to a "Moving Average Filter" which is the sum of the prior 10 samples):

This equivalence form helps provide a more intuitive understanding of the filter with further details provided below:
The basic form of a differencing filter is $y[n] = x[n]-x[n-1]$, with coefficients given as [1 -1] and is a discrete time approximation of a derivative (as the Forward Euler mapping of $s$ which is the Laplace Transform of a time derivative):

Also consider how the time derivative is $\lim_{T \rightarrow 0}\frac{x(t+T)-x(T)}{T}$ in comparison to the form given above.
This filter has a simple high pass frequency response with the magnitude response shown below, given as the response from DC to the sampling rate where the frequency axis has been normalized by the sampling rate. Like the derivative which has a frequency response that increases as a function of f (given the Laplace transform of a time derivative is $s$), this filter is most sensitive to the highest frequency components in the signal (at the Nyquist frequency at half the sampling rate given by the normalized frequency of $0.5$):

If we insert 9 zeros to upsample the filter by 10, to get coefficients given as [1 0 0 0 0 0 0 0 0 0 -1], (which is the differencing stage in the first diagram provided above) such a zero insert creates periodicity in the frequency response, repeating the same frequency response ten times over the frequency range from DC to the sampling rate as shown in the magnitude response below (creating the classic "comb filter"):

The final response of the OP's filter is found from the convolution (or cascade) of the interpolated differencing filter with a 10 sample moving average filter (which removes the high frequency sensitivity as a low pass filter).
A moving average filter over 10 samples with coefficients given as [ones(1,10)] would have a Dirichlet response which is an aliased Sinc function. This is clear when we consider the coefficients of the filter is the impulse response of the filter, and the Fourier Transform of the impulse response is the frequency response. The coefficients are a sampled pulse, and the Fourier Transform of a pulse is a Sinc function. The response of a filter with coefficients [ones(1,10)] is shown below:

And the final form of the filter as the convolution of the coefficients of the interpolated differencing filter and the 10 sample moving average filter is shown with the response below (omitting the scaling by $0.05$). Convolution in time is multiplication in frequency, and we see that the resulting response is the product of the responses given above:


How to remove certain background from image?

If the output of your Gábor-filter is reliable and the variation in your image-data is not so high that the outcome looks completely different, then you can use the following approach (parts of it where already mentioned):

Binarize your second image with any automatic thresholding algorithm. The range of  thresholds which will work is large as you will see.
Use a component labeling algorithm to label each connected region of pixel with a unique value.
Calculate for every component of your image a property which describes, how close your object is to a filled circle. For this, you can use for instance the compactness. I used below the fraction of pixels within the equivalent disk radius. This radius is the radius a disk with the same area as your object would have.

Area is simply counting the number of pixels of your object
Calculating the radius of a circle given its area should be manageable too
For the pixel inside this radius, you need the center of each image object, but this is just the mean of all positions of the object-pixels. If you are a physicist then you know the center of mass of several point-masses. This is equivalent.
Now you calculate for every pixel of each object the distance from its center and check whether or not it is smaller then the circle radius. Dividing the two numbers and you get the fraction of inside and outside pixel.

Take the object with the highest fraction. Note: objects with only one pixel will get a value of 1. So you should set a size-threshold and take only objects larger, say 10 pixel.

To see, that an automatic binarization should work, here are the results for a very low and a very high threshold:

Update Component Labeling

The choice of the component labeling is not critical for your application. I would suggest, if you have to implement it by yourself, to use a very simple method. The one pass version from the Wikipedia site is very easy. Basically you iterate through your binary image and when you meet a pixel which is white and is not already labeled you use a new label to for this object and start with this pixel.

The process of labeling this object with a label is basically similar to a flood fill. This is on the Wikipedia-site the inner steps 1-4 in the algorithm. You start with this labeled pixel and put all its neighbors on a stack (they used a vector). For a pixel on the stack you check whether it is foreground and not already labeled. If you have to label it, you put again all its neighbors on the stack. Do this until your stack is empty.

Then you continue your scan through the image. Unlike the description on the Wiki-site, you don't have to remove a pixel from your original image, you just skip when you have a value different from 0 in your label-image.


Amplitude estimation of sinusoid in known spiky spectral noise

The best fit time domain solution can be found by constructing two two basis vectors with your known frequency and calculate the coefficients directly.  The magnitude and phase can then be directly determined from these values.

Let C be a vector of cosine values over your frame and S be a vector of sine values.  You then want to find $(a,b)$ so that $aC + bS$ is a close as possible to your function.

$$ F = a C + b S $$

Dot this with the basis vectors.

$$ C \cdot F = a C \cdot C + b C \cdot S $$

$$ S \cdot F = a S \cdot C + b S \cdot S $$

Since the dot products are scalars, this is simply a linear system of two equations with two unknowns.  If you have a whole number of cycles in your frame then $ C \cdot S $ and $ S \cdot C $ become zero and the solution is trivial.

Note, this operation is the equivalent of a single bin in a DFT with the right parameters.  That is, pick a DFT frame with a whole number of cycles of your known frequency, and look at the corresponding bin.

The general case, without cycle alignment, can also be solved in the frequency domain with way fewer calculations.  See my blog article:

Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1

Do this for each of your peaks.  For better results, subtract out the peaks you have found before processing the others.


what are the Pros and Cons of ESPRIT and MUSIC algorythms?

The prime thing such algorithms aim to do is to make use of more information that you may have about the signal. In this case, the extra information is that you know the number of signals (sinusoids) present in your measurements.

One pro for both is, therefore, when your measurements match the assumption, you get a more accurate representation of the spectrum: it has the right number of lines. The FFT does not guarantee this.

A con is that, if your assumption is wrong (there are three lines instead of two), then these algorithms perform worse.

These algorithms do not give better resolution than the FFT. Only more data can give more resolution.

There is a good exposition here comparing the two to each other.


The complex conjugate is applied in the frequency domain, which corresponds to flipping the signal backwards in the time domain. When you have $$f(x)\psi^*_{s,u}(x)$$ this is a multiplication in the frequency domain. Multiplication becomes convolution in the time domain. So you should use conv(f,fliplr(psi))


What is the absolute phase and the relative phase, of a signal?

"Absolute phase" is really relative phase to a point in time at which the phase measurement is referenced (Noon GMT on Tuesday, the first sample in a time-domain window of sample data for instance, or maybe the middle sample in a vector).
Relative phase can also be measured with respect to another waveform of the exact same constant frequency, as the difference in phase measurement between the two waveforms at any one point will be the same at all other points, as long as the two frequencies stay identical.

Phase is important because it is one of the 3 parameters needed in an equation that completely describes a sinusoid.  Without phase information, the waveform is ambiguous.

$$y(t) = a\cdot \cos\left(b\cdot t + \phi_{t_0}\right), \quad \text{ where }\phi_{t_0} \text{is the phase at time } t = t_0.$$


In case your signal is fully stationary, just go ahead and use a "sinusoidal fitting" and take out the frequency that is not off your interest.
Then you have a disturbance free signal. Although since these signals are, I guess, not multiple integers of each other, your main frequency will also be low pass filtered. You can also use a notch filter, which will have similar behavior.

Please update us if you found a good answer.


Question about Fourier descriptor and the physical meaning of the coefficients

This question overlaps, to an extent, with this one. The best way to understand the Fourier Descriptors is through a basic understanding of the Discrete Fourier Transform.

The "curve" is basically a complex signal at the input of the DFT that is decomposed to a sum of wavy ellipsoid rings of increasing frequency.

Another convenient way to think about this is that what we usually see as "amplitude vs time" has now become "radius vs angle". Amplitude of a quantity versus time is how this quantity increases and decreases over time. In the radius vs angle way of thinking, you map the amplitude to the radius of an ellipsoid and the angle to a point on its circumference. Just as you read a diagram left to right, you "read" a curve beginning at some point and proceeding to the next either clockwise or counter clockwise.

  ...where this DC component is in the descriptor vector, since its a even number of values. Is it the term with the highest value?

The DC component is the first "bin" of the DFT. Given some signal $v[n]$, its DFT is $V = \mathcal{F}(v)$. But $V$ is usually indexed by $k$, as in $V[k]$. The $k^{th}$ value of $V$ is called a "bin". $V[0]$ is what is known as the DC component and it is basically a constant value. Because of the way the DFT is evaluated, the $V[0]$ is basically $\frac{1}{N}\sum_{v \in V} v$ and because $v$ is complex, the result can be complex too. But since you have encoded your $x,y$ curve as a complex signal of $v[n] = x[n] + j y[n]$, that $V[0]$ is basically representing the centroid of the curve.

  If I recalculate the shape with only 2 descriptors I obtain a circle around the mean position with a certain radius. Is the radius then the second highest value?

It is difficult to know what your reconstruction looks like without a specific example. The first component will be an ellipsoid which represents the mean radius through all points. If you keep adding components, all that you do is adding "squiggliness" to the average ellipsoid until a point where that original ellipsoid has become as squiggly as the curve.

Here is an excellent example:

From this paper.

Every value of $V$ is how quickly the radius varies versus the angle argument around the curve. So every value of $V$ is a "radius".

  I'm mainly interested in the ellipses that is created with recalculating the shape with 4 factors. But I would also like to understand the whole concept correctly. I read that the $A_{-1}$ and $A_1$ coefficients represent the major and minor axis of this ellipses. And that with $d = 2(|A_{-1}| + |A_{1}|) $ and $e = 2(|A_{-1}| - |A_1|)$ the axis can be calculated in a length unit. Do I have to scale these values somehow? Because if I pick the third and fourth highest values, the result is not what I expect.

What you are probably experiencing there is the scaling factor of the DFT.

I do not have Gonzalez' book handy right now but I will check and amend this part if it refers to something else, especially with respect to the way $A$ is indexed here.

In general there is a $\frac{1}{N}$ scaling factor for the DFT which is usually applied in the inverse transform (see for example here and here. So, the forward transform evaluates sums directly and this is why you end up with $V$s that are very high compared to the original signal. If you want to use those frequency coefficients directly, you have to divide them by $\frac{1}{N}$.

But, apart from the DC component, the rest of the ellipsoids do not have a "radius".

What are you trying to measure?

Hope this helps


Demodulation of OFDM symbol in LTE for multi-user case

The zero padding (or the cyclic prefix) takes care of inter-symbol interference that occurs due to delay spread. It is hence something you append to your signal in time domain. Hence the order is this:

Assemble all the data symbols along the subcarriers (can contain data for multiple users, as you said).
Compute DFT (via FFT, typically)
Add cyclic prefix (or zero padding)

Then transmit.


Harmonic partials have frequencies $f_i=if_1$, $i=2,3\ldots,$ where $f_1$ is the fundamental frequency. So the frequencies are not doubled, but they are integer multiples of the fundamental frequency. But if you talk about periods, then you have $T_i=1/f_i$, i.e. $T_i=T_1/i$, $i=2,3,\ldots$. So the period of the $i^{th}$ harmonic partial is the fundamental period divided by $i$.


Bandpass butterworth filter frequencies in scipy

You're probably running into numerical precision issues on the higher-order filters with sharp cutoffs. It could be a limitation of the butter function in SciPy, or it could be due to the filter structure that it's using when evaluating the frequency response. A better way of implementing high-order IIR filters is using a second-order sections (SOS) structure.

A little Googling suggests that SciPy doesn't yet support converting a filter to its SOS representation. Furthermore, there are apparently some issues with how filter parameters are handled internally that cause bad results with high-order filters. That sounds like what you're describing. These are apparently in the process of being addressed.

With that said, I came across this example that shows how you might implement an SOS structure using SciPy. I haven't tested the code at all, so there's no warranty, but it may be worth a look.


The Nyquist criteria refers not to the frequency, but to the bandwidth, which is related to information density in a signal.  A very high frequency signal, of approximately known frequency, with a sufficiently small bandwidth, will still be aliased or folded down with baseband frequencies by undersampling.  But if the bandwidth (or other known characteristics) of the signal is known to rule out other aliased frequencies (such as the non-existence of baseband spectrum content), that knowledge in conjunction with the undersampled samples may still provide enough information to reconstruct the signal.


Why do we use sinc interpolation rather than just multiply the frequency domain signal by rect?

The "two" operations you talk about are one and the same:  low pass filtering using the $\mbox{rect}()$ function in the frequency domain and $\mbox{sinc}()$ interpolation in the time domain are the same.


RMS value of a signal with a DC offset in MATLAB

RMS values don't simply add up. The RMS value of $N$ discrete samples $x[n]$ (assuming $0\le n<N$) is

$$x_{RMS}=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}x^2[n]}\tag{1}$$

Now if you have $y[n]=x[n]+c$ with some constant $c$, the RMS value of $y[n]$ is

$$\begin{align}y_{RMS}&=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}y^2[n]}\\&=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}(x[n]+c)^2}\\&=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}(x^2[n]+2cx[n]+c^2)}\\&=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}x^2[n]+\frac{2c}{N}\sum_{n=0}^{N-1}x[n]+\frac{1}{N}\sum_{n=0}^{N-1}c^2}\tag{2}\end{align}$$

If we define the mean value of $x[n]$ as

$$m_x=\frac{1}{N}\sum_{n=0}^{N-1}x[n]\tag{3}$$

then from $(2)$ the square of the RMS value of $y[n]$ can be expressed as

$$y^2_{RMS}=x^2_{RMS}+2cm_x+c^2\tag{4}$$

So if you assume that the mean of $x[n]$ equals zero, the squares of the RMS values add up:

$$y^2_{RMS}=x^2_{RMS}+c^2\tag{5}$$

You can easily verify Eq. $(5)$ using your example:

$$(282.7014)^2+(200)^2=(346.2948)^2\tag{6}$$


Relation between frequency spectrum and PDF of a random variable

Let's take this from the bottom:

  Is this a reasonable question?

Of course!

  In other words, would the PDF of a random variable constrain the possible frequency spectra the resultant random process may take?

That depends. For example, if the PDF of your RV was a Dirac impulse, i.e. your RV was actually a constant, then, obviously, yes, that will influence spectral aspects of your random process.

Now, in general, a random process is random, i.e. one can not say what its spectrum is going to be (or whether that is possible to define at all).

Hence, you can only do things like considering a specific statistical/stochastic property the random process has.

That, for example, and also most likely closest to what you're looking for, would be the Power Spectral Density (PSD).

Now, intuitively, one would just say that you could take infinitely many realizations of a random process, transform them to the frequency domain (using the Fourier transform), and calculate the expectation value of these transform's magnitude squares.

However, as Matt hinted at, it's really not easy to show that this transform, or the expectation value, are well-defined for some types of processes (or if you can't make some strong assumptions on the process' properties).

Hence, mathematically, you take the elegant detour of first calculating the autocorrelation, which, again, is a property of the process, and Fourier transforming that.

For particularly  well-behaved processes (mainly: Your process needs to be wide-sense stationary), this identical to the aforementioned expectation of the magnitude squared Fourier transform across all realizations; that's the Wiener-Chintschin theorem, which I mentioned in another answer that might be of relevance to you.

The point here is: Together with some statements on the autocorrelation of the process, and by observing the PDFs at fixed times, you might come to some statements on the process' PSD.

  When looking at the frequency spectrum of the generated data does the knowledge of the exact distribution function have any constraints/give any information about the frequency spectrum of the samples?

No. This is exactly the case where lack of ergodicity means that you can't derive much from a time ensemble.


Sinc Interpolation implementation in Matlab for reconstructing or interpolating?

I haven't done this in MatLab, but have done it using C++. So I am not in a position to provide you any matlab code example or snippet, however I will direct you to some links that will help you clear your doubts and concepts regarding the sinc function, filtering and reconstruction.

Sinc function is basically used for low pass filtering signals and for reconstruction of the signal using the previous or old samples. Check out these links, It will clear your concepts about sinc function such as how its used for low pass filtering and reconstruction.

concept of sinc function :- http://www.dspguide.com/ch11/2.htm

Low pass filtering :- http://www.analog.com/media/en/technical-documentation/dsp-book/dsp_book_Ch16.pdf

Reconstruction is done through interpolation, check :- What is an algorithm to re-sample from a variable rate to a fixed rate?

There are already loads of discussion here about sinc functions, sinc windowing filters, interpolation techniques, reconstruction etc. Reading them will certainly help you clear your doubts and concepts regarding these. Hope this helps.


Beginner level: deconvolving Gaussian noise to find initial PDF

You have a random variable $X$ whose pdf is $f_X(x)$. when another random variable $N\sim f_N(n)$ is added to $X$, then the pdf of the sum is just $f_Y=f_N*f_X$, where $*$ is convolution. This can be written in the Fourier domain as $F_Y=F_NF_X$. So if you know both $F_Y$ and $F_N$, and your $F_N$ is well-behaved, then $F_X=\frac{F_Y}{F_N}$.

After that, you can acquire $f_X$ by inverse Fourier transform.

Start by looking at the histogram of $y$, by which you can find $f_Y$. I think you can easily do the rest, given the above information.


Parameters of Gaussian Kernel in the Context of Image Convolution

When dealing with Gaussian Blur in the Image Processing context the following holds:

The Standard Deviation, $ \sigma $, is sometimes called radius. I think this goes back to Photoshop.
If you implement this using FIR Filter (Well, Gaussian Kernel is infinite so you approximate it) usually the radius of the filter will be something like ceil(4 * kernelStd). This is because the Kernel almost vanishes after $ 4 \sigma $.

Some in practice remarks:

The Gaussian Kernel is Separable hence if implemented as FIR filter it is implemented as 2 1D Convolutions - Along rows and along columns.
You may have a look at the Image Convolution GitHub Repository for efficient C Code.
There are many other way to approximate the Gaussian Blur which are much faster (Though sometimes there is accuracy vs. speed balance, most of the time negligible). You may have a look at Fast Gaussian Blur GitHub Repository for IIR and Box Blur implementations which are insensitive to the radius parameter.


Besides what is discussed below (see for instance the above Cris Luengo's comment): for reasonably-sized images, and modern computing techniques, I can see no reason in practice why the basic template match counting could be slow, and why one should look for faster, imprecise and incomplete alternatives.

[UPDATE 20180730]
You are trying to compare matches in a pair of matrices. Ontologically, the comparison operation  is not a linear operation. Hence, no convolution operation, naturally linear, seems be the easiest approach. Neither the FFT or Fourier stuff in general, as they are quite related to linearity.

IMHO, you ought to look at the concept of Template Matching, and avoid the notion of kernel, which is a debated notion.

Instead, a simple counting approach could do the job. For instance in Matlab:

%% The OP kernel and image
% dataKernel = [1 2 1;1 2 1;1 2 1;];
% dataImage = [1 1 2 1 2;1 1 2 1 2;1 1 2 1 2;];

%% A random kernel and image
dataKernel = floor(256*rand(64,64));
dataImage = floor(256*rand(300,300));

%% Compute sizes
[nRowKernel,nColKernel] = size(dataKernel);
[nRowImage,nColImage] = size(dataImage);

[nRowCrop] = nRowImage - nRowKernel+1;
[nColCrop] = nColImage - nColKernel+1;

lMatch = zeros(nRowCrop,nColCrop);
dataImageCrop = zeros(size(dataKernel));

tic;
for iRowCrop = 1:nRowCrop;
    for iColCrop = 1:nColCrop;
        dataImageCrop = dataImage(iRowCrop:iRowCrop+nRowKernel-1,iColCrop:iColCrop+nColKernel-1); % Crop a portion of the image
        lMatch(iRowCrop,iColCrop) = length(find(dataImageCrop - dataKernel == 0));
    end
end

Here, I just used an ad-hoc displacement for the crop, but than could be automated easily. No product, just differences. If your data are limited size integer, there might be more options for speed. But even in Matlab, this is relatively fast.

FFT-based cosine similarity for fast template matching seems to better explain what is alluded too in OpenCV performance on template matching (provided by the OP). If the image is of size $I_x \times I_y$, and the template is  $T_x \times T_y$, with $T_x < I_x$ and  $T_y < I_y$, you can expect a number of operations of:

$O(I_xI_y  T_x T_y)$ with the standard comparison,
$O(I_xI_y\log I_x \log I_y)$ with a Fourier method.

From afar, Fourier might seem to be beneficial whenever $\log I_x \log I_y <  T_x T_y$. But:

this is only an asymptotic result the constant behind the $O(\cdot)$ can make a huge difference
then you have to find maxima, which add to the overall cost
the hardware implementation (pipeline, CPU/GPU) can made a lot of difference in practice, compared to number of operations
indeed, the first one only require adds/subtracts and absolute values on integers,  the second complex operations on floats; thus, there might be huge variations between implementations
a difference of pixels, like the sum of absolute difference (SAD) does not translate easily in Fourier because it is orthogonal
implementing fast correlations, the FFT might be able to provide you with a "best match", but it is not local, so I cannot imagine right now how this can directly give you a count by patch location.

So far, to me:

if you only want a best match, then you can do Fourier on the whole image;
if you really want (even approximate) counts at all template locations, I do not see an easy and useful Fourier implementation for that so far.

However, you can speed the matching up, with several tricks, like using the symmetry of the template (if any), the integral image, multiscale pyramid-like matching, etc. A few additional pointers:

Fast Template Matching
Fast Template Matching With Polynomials
FasT-Match: Fast Affine Template Matching - The Computer Vision ...
Very Fast Template Matching


Can the magnitude of DFT of a triangular wave be max for at 0?

What you see in the red graph is not an impulse ($\delta$) at zero frequency (in such case we would have a constant zero frequency which is equivalent to a constant dc signal $x[n]=c$).

What you see is actually a squared $\mathrm{sinc}$ function. You can easily verify that the frequency representation of a triangular function (both in continuous and discrete time) is in the form of $\mathrm{sinc}^2$. Look at here, for instance for the continuous triangular signal to see that

$$\mathcal{F}\{\mathrm{tri}(t)\}=\mathrm{sinc}^2(f)$$

and $\mathrm{sinc}$ has a maximum at zero.


Hk are modulations of the low pass filter (band pass instead of low pass).

$$
(W_{4}^{k}) = e^{-2j\pi k /4}.
$$
For $z = e^{j\omega}$: $$H(W_{4}^{k} z) = H(e^{j(\omega-2\pi k/4)})$$
This means that the filters $H_k$ are shifted in frequency- these are the band pass filters you want to get using your filter bank.
For $k=0$ $H_0$ will pass the frequencies $[\frac{-\pi}{4} \frac{\pi}{4}]$, for $k=1,[\frac{\pi}{4} \frac{3\pi}{4}]$ etc..

In the DFT filter bank scheme,
$y_0[m]$ are the outputs from $H_0$ , $y_1[n]$ are the outputs from $H_1$ and so on.


How to estimate the number of taps needed for subsequent channel estimation algorithms?

The length of the impulse response is typically related to the frequency resolution of the channel transfer function. As a rule of thumb: the more detail there is in the frequency response, the longer the impulse response will be.

In practice there are a few things you can do: If you have full access to a similar, you can simply measure it with a very long impulse response measurement. Than you can truncate the impulse response and see what happens to the transfer function. The truncation will create errors and this way you can dial in the impulse response length to the point where the error is still tolerable.

You can also use physical knowledge about the channel. For example an audio amplifier has only a few electronic components, all of which are specifically designed to create a flat transfer function with little phase distortion. A handful of samples is fine for that. On the other hand  look at a loudspeaker in a room: the sound bounces around with multiple reflections until it finally dies off. In this case you would need many thousands of samples (not practical at all).

Many systems have bandpass or high pass characteristic: all acoustic systems are high pass because air can't transmit DC sound. Most communication systems are band pass since the information needs to stay away from the extreme edges of the band. In these case often the length of the impulse response is determined by the high pass roll off, i.e. the frequency and steepness of the high pass.


Implications of the Shannon-Nyquist sampling :

a. Sampling rate below Nyquist rate: Aliasing happens: shifted spectra of the sampled signal will overlap; i.e., high frequency signal content will be aliased into lower frequency regions of the spectrum and low frequencies into higher, therefore, irrecoverable misinterpretations happen about the authentic (true) content of the sampled signal spectrum;

b. Sampling rate equals Nyquist rate : Depends on the signal type. If there is no impulse at the Nyquist frequency, or the spectrum is zero at the Nyquist frequency, then this will in principle satisfy Nyquist sampling rate, but requires an ideal anti-aliasing filter to be practically implemented and thus usually avoided and a slightly higher sampling rate is preferred instead.

c. Sampling rate larger than Nyquist rate : Satifies Shannon-Nyquist sampling theorem, allowing perfect and unique recovery of the original continuous-time signal from its samples. However, going too high, unless necessary will increase system cost. So practically it's best to limit the sampling rate to a satisfactory minimum or an achievable maximum.


Does it really take 6 minutes and 20 seconds to properly fade-out a single precision audio signal?

You start off with wrong assumptions, so your result is nonsensical.

Your definition of good and bad should be based on some psychoacoustic understanding of what type if distortion and artifact is actually audible and/or objectionable
The 24-bit mantissa argument (it if were valid at all) would only hold a fixed point number. In floating point you can apply pretty much any gain you want. The smallest non-zero gain is actually 1.1755e-38.

There are various different ways to design a fade: linear, quasi-log, "half-window" functions, all kinds of continous derivative functions etc. They all sound slightly different and the best choice really depends on your application: musical fade, cross-fade between two signals, emergency fade out that needs to be fast, etc.


Are high pass/low pass filters FIR filters?

The expressions Finite Impulse Response (FIR) and Infinite Impulse Response (IIR) refer to the time-domain response of the filter. When an impulse is fed into the FIR filter in the time domain, it's response goes to zero within a finite amount of time. The response of the IIR filter does not.
Low Pass Filter (LPF), Band Pass Filter (BPF), High Pass Filter (HPF), Band Stop Filter (BSF) and Notch Filter (NF) are expressions that refer to the frequency-domain response of the filter - which frequencies it passes, stops or attenuates.
Both domains can be used to classify a filter's behavior. You can implement a LPF as a FIR filter. In the frequency domain, it will pass low frequencies and stop high frequencies. Simultaneously, in the time domain, the impulse response of the filter will return to zero after a finite amount of time. You can also implement a LPF as an IIR filter. It will have the same behavior in the frequency domain, but its impulse response will never reach zero in the time domain.
The filter's response in the time and frequency domains complement each other.
I hope this helps.


Firstly, you need to know that the FFT computes the DFT (discrete Fourier transform) in an efficient manner. Hence, the output of an N-point FFT and N-point DFT are exactly the same.

Now, especially, if N is a power-of-two, the FFT can be calculated very efficiently. Though, if N can be factorized into rel. small prime numbers (2, 3, 5), still the calculation can be done very efficiently. Strictly speaking, the Radix-2 FFT only works, if N is power-of-two. However, in common implementations in software, N can be any number and the implementation chooses the appropriate algorithms for most efficient calculation.

Hence, if N is power of two, the calculation by FFT is just a bit faster, compared to when N would be a large prime number.

Regarding sampling frequency: If $F_s$ is your sampling frequency and $T$ is the duration of the FFT window, then $N=TF_s$. It's hard to answer more detailed what you mean by relation to sampling frequency, because I dont know exactly what you are asking.

Regarding your comment: The DFT maps a subset of the frequency range $[0,\dots,F_s]$  to the $N$ DFT-bins, which are indexed by $k=0,\ldots,N-1$. Since all frequency bins in the DFT output have the same distance to each other (i.e. they are linearly spaced), you immediately end up with $f=\frac{k}{N}F_s$, i.e. the $k$th bin corresponds to the physical frequency $f=\frac{k}{N}F_s$.


wav files generally have uncompressed data. That means you can get the exact waveform of the speech signal in digital form. You can tell whether the signal has some speech or not by calculating energy and using threshold to make the decision.

Again, the above technique may not work if there is noise along with speech.

Speech signal is quasi stationary signal and same speech will have different waveforms based on who is speaking, his pitch, his slang, his emotion and many more factors. So, it is very difficult to match two speech files.

The two problems that you wanted to solve

% of match between two files.
Detect noise.

are special fields in signal processing by itself. First one can be achieved by doing speech recognition on both the files and get the exact text for both files and then decide some rules to get % of match. For example, even though waveform of both speech are different but content is same, only with speech recognition you can 100% match.

coming to noise in the speech signal, there is again separate field in signal processing called as noise reduction. even though the process of noise reduction is used to reduce noise, one of the steps in noise reduction is to detect noise. You can use that step to detect noise.

But, both these topics are pure signal processing topics and i don't think you can achieve these with some simple function call or some simple python code.


Is there an algorithm to compare several pictures?

It looks like you are looking for algorithms for Content-Based Image Retrieval. As long as you don't expect any miracles, you should be quite happy with a colour histogram approach. Have a look at MPEG-7 colour descriptors.


Can we really use equal gain combining in the receiver in the real world?

Because we have to know the phase of elements in the channel matrix first before we do the EGC,however,the channel should be randomly generated every second,so I wonder that if we can do EGC in the real world or not?

Yeah, so you continuously have to estimate the channel. You need to do that anyway, no matter whether you use any combining method or not, because phase is typically important to the receiver!

So, think of EGC as a kind of channel knowledge-based beamforming.

Therefore: Yes, EGC is used. The difference to MRC is that there's no need to estimate the receive path's SNR – and depending on the signal and receiver model, that SNR estimate might be non-trivial (e.g. when you don't know the different noise figures of your receivers, or when you can technically easily know the phase of your receivers, but not the strength of the signal, e.g. because it's very weak, before you used combining).


Envelop detection using matlab hilbert function

This must be an artifact of the implementation for the hilbert function, in that any implementation of finite length will have a finite ripple. It doesn't appear the underlying filter using in the hilbert function can be modified, so there are three suggestions to obtain a smoothed demodulated result:
Method 1: Implement the Hilbert using the FFT: take the FFT of the original signal, set all the negative frequencies to zero (the upper half of the FFT) and double the positive frequencies leaving bin 0 (DC) as is (and as Overlord pointed out in the comments, leave the Nyquist bin at $N/2$ as is when the FFT length $N$ is even). The IFFT will be the Hilbert Transform and taking the absolute value of this signal will be the desired envelope. (As @aconcernedcitizen pointed out in the comments, this IS the method specifically implemented by Octave using the hilbert function from the signal package.)
Method 2: Implement the Hilbert with quadrature phase tracking filters where the filter length versus ripple can be traded.
Method 3: Use a traditional AM demodulator by multiplying the modulated signal with the coherent carrier and then low pass filtering the result.
A demonstration of Method 1 is shown below, which is a suitable solution for the OP's case:

Zooming in on the start of the signal shows the distortion limited to the start-up condition and the smooth envelope after that.

The Matlab code for this is as follows:
N = length(sig);
sig_spectrum = fft(sig);
hilbert_spectrum = zeros(1, N);
hilbert_spectrum(1) = sig_spectrum(1);
hilbert_spectrum(2: ceil(N/2) - 1) = 2 * sig_spectrum(2:ceil(N/2) - 1);
if mod(N, 2) == 0
  hilbert_spectrum(N/2) = sig_spectrum(N/2);
endif
hilbert_time = ifft(hilbert_spectrum);
envelope = abs(hilbert_time);


What should be done if the mask-values are too small?

A few of the stuff that appears in this response were initially packed together in a compact comment. Inevitably, in this answer, I would have to provide code, which I wouldn't like in this case. You have to get what is going on. So, I am going to try and talk a little bit more "around" the code and hopefully by the end of this, the code will have written itself.

In just a few words: $H$ is your frequency response in the frequency domain, it is based on a bandpass pattern but it is supposed to be applied in a directional way. So, build a bandpass and rotate it.

Now then...Bandpass "pattern". You need something that goes up and then down and when it goes up it stays up as flat as possible. Take a plot of some $y_1 = \frac{x}{a}$. Well, obviously, that's a line with a slope of $\frac{1}{a}$ and it crosses zero at zero. Nothing fancy going on there.

Now, offset that zero a bit to the right, so $y_2 = \frac{b-x}{a}$. What did you just do? You brought zero bang on on $b$ but also you introduced some NEGATIVE values for $x>b$.

To see why this will come handy, do a $y_3=(\frac{b-x}{a})^2$. This is a parabola. It goes down and then up. It is also curved, let's remove this curve by undoing the square with a root, so $y_4=\sqrt{(\frac{b-x}{a})^2}$. What's going on now? We have two lines, one that arrives down to $b$ and one that departs upwards from $b$. Let's now reverse these trends.
$y_5=\frac{1}{\sqrt{(\frac{b-x}{a})^2}}$.

What is any $u=\frac{1}{v}$? A hyperbola. Remember those two lines that meet at $b$? If you were to substitute $v$ for $y_5$, then, when $y_5$ goes down, towards $b$, it tends to smaller values, therefore, $u$ goes upwards!. And when $y_5$ departs from $b$ towards increasing values then $y_5$ goes downwards.

What happens bang on on $b$? That's when $y_5$ hits zero and if $v$ hits zero then $u$ goes to infinity. We don't like that. So, let's make that $v$ a sum, to control the denominator and ensure it never hits zero. $y_6 = \frac{1}{1+\sqrt{(\frac{b-x}{a})^2}}$. So now, when that root term goes to zero, $y_6$ stays at $1$.

Let's deal with that stay flat business. Remember that $\sqrt{q^2}$ trick we did to get rid of possible negative values of $q$? Well, we could have done this with $|q|$, which means taking the absolute value of some $q$. But this power raising has an added benefit. So, in general, if you raise to some power, you increase the slope of a function. That is, how sharply the curve evolves. Take some $q^2$ and $q^4$, for any two $q, q+\delta_q$ where $\delta_q > 0$, then, $q^4$ grows faster than $q^2$.

BUT!, raising to a power has another property as well, raising large numbers to a power makes them larger and raising small numbers ($ < 1$) to a power makes them...even smaller!. Why is this important for us? Because as our $y_2$ approaches its $b$ limit, it grows smaller and smaller and smaller and we can accelerate this decrease by raising to some power. Therefore, we now go $\sqrt{q^n}$ and use $n$ to control how fast that function will approach small values and stay small leading to that flatness characteristic.

Now, here is the thing with multiplication, if you raise some $q$ to an odd power, you retain the original sign of $q$ and if that was negative then $\sqrt(q)$ will become complex. We don't like that in this case. So, what do we do? We ensure that some $q$ is raised to an even power by multiplying by two. All together now:

$$y_7 = \frac{1}{1+\sqrt{(\frac{b-x}{a})^{2n}}}$$

OK, fine, how do I apply $y_7$ as a filter to a signal? Pretend that your $x$ is now frequency, mirror $y_7$ around some $x_{\frac{N}{2}}$ (where $N$ is the full length of your impulse response), take its inverse Fourier Transform, throw away the complex part (which should have been zero anyway) and treat what you get back as your impulse response. That is your $h$.

Hold on, if $x$ now represents frequency then what are the $b,a,n$? Consider how $y_7$ goes up towards something and then comes down from that something. So, $b$, is your Center Frequency. The centre frequency of the bandpass filter. $a$ is now your bandwidth, that is how fat that lobe of values that go up and down is. And finally, $n$ is now your filter's degree, that is, how flat it tends to be. In fact, if you raise $n$ to some absurd value (which, remember, it will be doubled), then it will look very much like a square pulse which would be an ideal filter. But, the more your crank up $n$ the more Fourier coefficients you need to synthesize it accurately, and then the impulse response gets longer (that $N$ from earlier) and...anyway, just keep that in mind, it's not important for now.

WHAT THE?!?!? HOW DOES THIS APPLY TO MY ORIGINAL QUESTION?!?!!?!?!?!

We are very close now, all that we have to do is bring in another dimension, so that our $y_7$ now gets a $y$ as well. (I hope using $y$ and $y_7$ is not confusing, they are just symbols)

I am not going to go all over this again because we need to focus on what is really important here. When you bring in another dimension you can now define direction. I am not going to go through the full two dimensional Fourier spectrum explanation but I will just say this. When you look at the frequency spectrum of an image, the low frequencies are at the centre and the high frequencies are at the periphery (by convention).

Therefore, filters that treat all directions equally look like cakes or donuts. It's not hard to imagine why. Take our bandpass profile $y_7$ and mentally rotate it about the end where its zero lies. This will give you a ring (or doughnut) about the centre of the image. Pick any direction, what you see is the bandpass filter and its symmetric image, remember earlier that I said, "...mirror it around some $x_{\frac{N}{2}}$..."?.

But, hey, we now have the luxury to filter along a specific direction and not treat all of them equally. So, let's take our filter profile and apply it to just a specific direction.

How are we going to do this?

Remember our $y_7$? This is effectively an $f(x)$, a function of some $x$ variable which we can assign any meaning we like. It could be some generic $x$, it could be a frequency axis segment. But, at the end of the day, what we did there was to scan $x$ and assign a $y$ value, so the height of $y$ depends on the value of $x$.

What we are going to do now, is scan the frequency plane and assign a height to it, so imagine your $H$ as being a height map. But! We want to be able to control for the direction which translates to rotation.

So, while I am scanning some $u,v$ point on the frequency plane, I am assigning its height on some $u^*, v^*$ which are the rotated coordinates, around the centre $t_x, t_y$.

At this point, I have to say, I can see why the authors want to call $D_v, D_h$ their cutoff frequencies and $tx,ty$ their "sub-band centre frequency" but from the viewpoint of the above explanation, $D_v, D_h$ would correspond to the bandwidth and $tx,ty$ would correspond to the cutoff frequency.

Now, you work in C# and you probably have very good reasons for this but my quick test was performed in Octave (which is very similar to MATLAB). It is a piece of FOSS and that's where you could have just done a imagesc(H) and looking at your results we could more or less understand where the problem is.

I am not going to give you the full code (yet) because that's not going to do any good. But I will write here the line where I suspect you are running into trouble (straightforward translation to C#):

H(i+N2+1,j+N2+1) = 1./(1+(0.414*(sqrt(((su./Dh)+(sv./Dv)).^(2.*n)))));

I will also suggest an alternative way that you can achieve the same thing, much more straightforwardly, just re-using the stuff I have been talking about so far. Less lines of code, less bugs to worry about.

Now, you might notice that the authors in their paper, state that $H$ should be made symmetric around the origin along the orientation you are trying to filter. So, at the end of the code that is building $H$ you have an extra step (which I suspect you are not doing at the moment) to mirror it. If you try to mirror this function AFTER it has been captured in some $H$ matrix, it is going to be hard work.

For the following to work, I am assuming that in C# you are using some elementary image processing library.

So, write a function that builds you a $y_7$ from above. Nice, clean, simple, one dimensional bandpass filter in the frequency domain. Once you have that, MIRROR IT by something like $y_8 = [y_7(N_{y_7}:-1:1) y_7]$. The indexing convention is $StartValue:Step:StopValue$, so effectively, if your $y_7$ was "Y7VALUES" your $y_8$ is now "SEULAV7Y7VALUES", notice that we skipped the value of 0. Otherwise, the "Y" would appear twice and that would mess the spectrum.

So, your $y_8$ is now one row of your $H$. To turn it into an "image mask", replicate $N$ copies of $y_8$, one below the other. Once you are done, you have your $H$ which if you were to look at as if it was an image, it looks like two parallel lines, two fat bands, one somewhere to the left of the middle of the image and one to the right. This is your $H$ at zero degrees.

To generate the same $H$ at a different orientation. Simply rotate that initial $H$.

DON'T ROTATE A ROTATED VERSION!!!! (so if you did 45 degrees, don't rotate that another 45 to go to 90), because that brings in distortions due to interpolation. Your image processing library should have elementary rotation capabilities. If it doesn't....good luck, you have the Hough transform coming up in the same paper...

I hope this helps.


Why does SNR have an impact on the bit rate?

Yes, try to understand what is the concept of SNR  and it depends on what. Then, it will be easier to understand why SNR and BER increase and decrease in all cases.

In your case, OFDM, the relationship between the SNR and BER is given:

$\textrm{BER} = \frac{1}{2} \left(1 - \sqrt{ \frac{L/N \times \textrm{SNR}}{2 + L/N \times \textrm{SNR}}} \right)$

Where $N$ is the number of subcarrier and $L$ is the channel diversity. $\textrm{SNR} = p/\sigma^2$ .. where $p$ is the transmitted signal power, and $\sigma$ is the variance.

so now, you check the relationship between the type of used modulation with the above equation, (QPSK and N) , it's easier to explore the behavior of BER. remember also that probability of error when using QPSK is different about when using for example QAM and so on

good luck


Is it possible to simplify the convolution integral if the functions are non-zero in disjoint areas?

Link to the closed-form solution:

wolframalpha

From left to right: $f$, $h$, and $f*h$:


Why is CDMA not enough for wireless communication?

You're misunderstanding what CDMA, TDMA and FDMA do:

CDMA doesn't increase the channel capacity in any way. It's a MA, MA = Multiple Access mechanism.

In other words, it's just a way of dividing the spectrum among multiple users. No matter what you do, you can't get more data through a channel than the physics and math allows – and Shannon's channel capacity gives a pretty clear upper boundary for that.

All the mechanisms you quote hence only divide the available channel across some measure.

For TDMA and FDMA, as you notice, this is very intuitively understandable.

But to be honest, CDMA is not that different. Instead of sending a single symbol in a single symbol time slot, you have to multiply it onto a longer code sequence.

The maximum number of simultaneous users of CDMA is simply the number $N$ of orthogonal CDMA sequences.

$N$ orthogonal vectors (i.e. sequences) form a base to an $N$-dimensional vector space (which is the signal space). Vectors from an $N$-dimensional space (which is what you're going to transmit to transport one symbol) have inherently $N$ elements.

So, to make (up to) $N$ people talk at once, you need to "blow up" everyone's transmission's symbols (which are $1$-dimensional, each) to $N$-dimensional vectors, increasing the time it takes to send one symbol by a factor of $N$.

That's the same efficiency as using TDMA to split time into $N$ slots (everyone sends $1$ symbol every $N$ slots), or dividing the bandwidth into $N$ channels and using FDMA (everyone can only use $\frac1N$ of the available bandwidth, reducing the individual symbol rate by $N$).


Can FFT Substitute Correlation for This Case?

How did you plan to use the FFT in the first place?

If it is to naively look for peaks in the magnitude spectrum and deduce the frequency from it, this is a bad idea. As you have seen, this has a limited resolution in the lowest frequencies. Resolution limitations aside, looking for a peak in the FFT is a poor way of estimating periodicity since the true fundamental frequency might not be the one with the most salient amplitude.

However, you can use FFT as an intermediate step to efficiently compute the autocorrelation function (FFT with zero-padding of the input to avoid wrapping -> square of magnitude of coefficients -> IFFT). Looking for peaks in the autocorrelation function is a more reliable period estimator than FFT peaks.

Thus, given your requirements, you could use sliding windows of 256 samples, with a 61% overlap, compute the autocorrelation through the FFT, look for the autocorrelation peak in the [50, 125] samples range and use this as an estimate of the period. Because of the 61% overlap, this will give you one period value per second. Observe here that it is possible to use large analysis windows (256 samples, or more) and make them overlap so that we can still get a detection for every 100 sample. The downside is that if the period of the signal rapidly change, your detection will be "smudged".

I see a few other approaches for your problem that do not suffer from any resolution compromise on the temporal axis:

If the waveform of your input signal is known (square, sawtooth, sine...), you could try locking an oscillator of similar waveform to it through a PLL, and use the PLL frequency as your estimator.
If your signal is purely sinusoidal, you can fit a complex exponential to it (through eigenanalysis of the autocorrelation matrix).
An approach that would give a very high resolution (one decision per sample) would be to run your signal through a bank of 75 comb filters in parallel, with delays of 0.5, 0.51, 0.52 .. 1.25 seconds. Compute the energy of the output for each of these 75 signals ; the one with the highest energy gives you the period.

Note that your problem has very similar requirements to that of tempo detection (onset detection yields an envelope signal at 100 Hz or the likes, we aim to get a tempo map in the 60 BPM - 180 BPM, ie 0.33 - 1Hz range, with a temporal resolution of 1s). Historically, the comb filterbank approach was used ; while most recent solutions are based on autocorrelation or similar autocorrelation-based pitch estimators (yin) ; with analysis windows in the 5 - 10s range.


Two things come to mind.

First, let's address the low-amplitude noise problem. This one is a typical textbook case for a noise gate. To give you an extreme example, these are sometimes even used on guitar distortion effects which are essentially insane amounts of controlled overdrive. I have to admit that I failed to find any literature on creating digital noise gates, but they exist abundantly as plugins so maybe digging through some open-source projects will score you a good basis for learning how to implement them. In essence, you measure windowed energy of a signal, and if it falls below a certain threshold, you fade your gain to zero (usually decaying exponentials are used). In any case, I'm sure you know what a noise gate is and I'm beating around the bush here. Make one and use it.

Secondly, when you say roughly set the input gain for me automatically if it's way off, that sounds like a textbook case for using a fuzzy controller. I'm linking to the Wikipedia page because I'm feeling a bit lazy, but you can find tons of great papers out there just on the first ten pages of Google results. The concept is exactly as you described it. You'll have several (say five for the sake of a simple example) levels of how close the gain is to what you think it should be:

- Way too low
- Low
- About right
- High
- Way too high

Based on these you will have some sort of fuzzy response like:

- Increase at a high rate
- Increase at a low rate
- Do nothing
- Decrease at a low rate
- Decrease at a high rate

These are my thoughts on how I would do it. You still need some means of measuring your fuzzy control input (assessing how good current gain level is), but that's something for you to work out. I'd say that if you're riding very close to clipping, that's High, if you've clipped, that's Way too high, and so on. Depending on what the spec considers acceptable, fuzzy control can sometimes let you get away with very simple solutions.


Detection of signal having gaussian noise with non zero mean and a covariance matrix using Neyman-Pearson Lemma(LRT)

You're looking at a Generalized Gaussian Detector$^{\ \text [1]}$. The Neyman-Pearson detector decides $\mathcal H_1$ if:
$$
\frac{p(\mathbf x; \mathcal H_1)}{p(\mathbf x; \mathcal H_0)}> \gamma\tag{1}
$$
With $\mathbf w$ and $\mathbf x$ independent; you get here:
$$\mathbf x\sim
\begin{cases}
\mathcal N\left(\mu_w, \mathbf C_w\right) &\text{under} \mathcal \quad \mathcal H_0\\[2ex]
\mathcal N\left(\mu_s + \mu_w, \mathbf C_s + \mathbf C_w\right) &\text{under} \mathcal \quad\mathcal H_1
\end{cases}$$

The NP inequality in equation $(1)$ gives:

$$
\frac{\displaystyle\frac{1}{(2\pi)^{N/2}\mathrm{det}^{1/2}\left(\mathbf C_s + C_w\right)}\exp\left[-\frac{1}{2}\left(\mathbf x-\mu_s - \mu_w\right)^T\left(\mathbf C_s + \mathbf C_w\right)^{-1}\left(\mathbf x-\mu_s - \mu_w\right)\right]}{\displaystyle\frac{1}{(2\pi)^{N/2}\mathrm{det}^{1/2}\left(\mathbf C_w\right)}\exp\left[-\frac{1}{2}\left(\mathbf x-\mu_w\right)^T\mathbf C_w^{-1}\left(\mathbf x-\mu_w\right)\right]}>\gamma
$$
Taking the logarithm on both sides gives the needed test statistic $T(\mathbf x)$:
$$
\underbrace{\ln\left(\frac{\displaystyle\frac{1}{(2\pi)^{N/2}\mathrm{det}^{1/2}\left(\mathbf C_s + C_w\right)}\exp\left[-\frac{1}{2}\left(\mathbf x-\mu_s - \mu_w\right)^T\left(\mathbf C_s + \mathbf C_w\right)^{-1}\left(\mathbf x-\mu_s - \mu_w\right)\right]}{\displaystyle\frac{1}{(2\pi)^{N/2}\mathrm{det}^{1/2}\left(\mathbf C_w\right)}\exp\left[-\frac{1}{2}\left(\mathbf x-\mu_w\right)^T\mathbf C_w^{-1}\left(\mathbf x-\mu_w\right)\right]}\right)}_{T(\mathbf x)}>\ln(\gamma)
$$
You can continue from there. Also, have a look at this related question on NP detector.

$^{\ \text [1]}$: Kay, S. Fundamentals of Statistical Signal Processing: Detection Theory.


What happens to frequency domain components if I circular-shift the time domain samples?

If you circularly shift the array in the time domain, the DFT of the shifted sequence will have the same magnitudes in each bin, but the phases will be different, with $X[k]$ being transformed to $X[k]e^{jk\theta}$,
$k = 0, 1, \ldots, N-1$, for a fixed value of $\theta$ that I will leave for
you to figure out.


Distinguishing FIR and IIR from difference equation

Welcome to DSP! First of all I would recommend reading chapter 2, specially 2.2 to 2.5 of Discrete Time Signal Processing 3rd edition by Oppenheim and Schafer. Then 3.2 and 5.2 of that very same book.

About your question. Well, in general (see comments on Axel answer for an exception) if you see an equation which depends on past output terms, like $1.2*y [n−1] $ in your equation then it is an IIR system. If it only depends on current, past or future inputs, it is an FIR system.

If it is an FIR system, you may determine stability, causality, linearity and time invariance from the Linear Constant Coefficient Difference Equation. However, if it is an IIR system, unless initial conditions are given to you, there's no way you could prove any of these because you don't know the value of the past output the first time you are going to "use" the system and it could determine whether it behaves one way or another.

For a system for which the input and output satisfy a linear constant coefficient difference equation, if the initial condition is that the system is initially at rest, then the system will be linear, time invariant, and causal. For stability analyze the Region of Convergence of the $Z$ transform and make sure it includes the unit circle


One word for that technique is superresolution.

Robert Gawron has a blog post here and the Python implementation here.

Usually, this technique relies on each image being slightly offset from the others.  The only gain you'd get from not moving between shots would be to reduce the noise level.


Is there two Types of Scale Space in SURF algorithm?

They implemented each one separately, then compared their results to each other.

FH-9  --->  it's a detector that it's initial filter size is 9x9

FH-15 -->  it's a detector that it's initial filter size is 15x15


For any estimate of the coherence, the cross spectra and spectra need to be averaged before forming the estimate of the coherence.  If there is no averaging, the squared coherence will always be 1.  You can see this by plugging in a single Fourier coefficient for X and Y: $C_{xy}^2=(\alpha_x\alpha_y)^2/(\alpha_x^2\alpha_y^2)=1$.  If you instead use a sum or average of two Fourier coefficients for X and Y, the answer can be less than one.

In multitaper spectral analysis, this averaging is done by averaging spectral estimates made using the same time series but different taper windows.  If you only use one taper window, there is no averaging being done.


remove spatial sinusoid noise from 2d image

Here is one approach to try:

This problem is about finding the correct re-indexing algorithm for pixel values at every raw.

For each raw we can search for possible shifts optimizing some cost function that represent the correct alignment. The cost function will take the raw intensities and two neighboring raw's (upper and lower raw's) as input parameters.  The output of the cost function should represent a smoothness/continuity measure.

Possible cost function to experiment with:

Minimizing total variation
$l_{1}$ norm minimization in some transform basis.

Constraints on the sine function can be used to improve the optimization algorithm. If you assume a narrow range for pixel shifts you should be able to use a brute-force approach ( searching for all possible shifts).


why would you want to calculate the PSD of 'the random process' and not just a single OFDM symbol – surely that's the granularity of focus.

Because you usually don't care about the single symbol as much as you care about making statements on e.g. the channel in general. So, you estimate the PSD.


Re-implementing the "spectrogram" function from matlab

Differences between your implementation and what the spectogram function is displaying.

This looks wrong and you don't preallocate the mag_spectrum matrix (which is bad practice).

mag_spectrum(:,idx*1024) = abs(spectrum).^2;

mag_spectrum(idx, :) = abs(spectrum).^2; % Use this instead

mag_spectrum(:, idx) = abs(spectrum).^2; % Or this... depending on how you want the data to 'waterfall'
You don't have the proper scaling (see reference paper and code below). Also, I believe the spectogram is displayed in dB.
You're saving the two-sided FFT. This is unnecessary and the spectogram defaults to the one-sided FFT. Use the proper scaling accordingly.

Here's some code used to compute the average PSD of an input sequence using a provided window and some other parameters. You can modify the last few lines to return the entire spectogram by getting rid of the call to 'mean' and properly handling the last step of reindexing.

Performing the overlap is pretty easy if you convert the data vector into a matrix using reshape, and then proceed to concatenate that matrix with a shifted version of itself.

A good explanation of some of the subtleties of the PSD is presented in this paper, especially regarding scaling.

[px f] = welch(randn(1,32768), 1024, ones(1,1024), 512, 1024, 0);
plot(f, 10*log10(abs(px)))

function [Px, f] = welch(x, fs, win, overlap, Nfft, oneside)
% x=vector of data samples
% fs=sampling frequency (in Hz)
% win=vector containing the window
% overlap=number of points to overlap
% Nfft=length of FFT to use
% oneside=a logical variable indicating whether
%    scaling is for 1-sided or 2-sided PSD (true if 1-sided)
% Px=estimate of the power spectral density
% f=vector of frequencies (in Hz) for plotting
if nargin < 6
  error('Insufficient input arguments')
end

x = x(:);
win = win(:);

Nw = length(win);
Ns = length(x);

% Compute the frequency vector and scale factor
if oneside == true
  sfactor = 2 ./ (fs * sum(win.^2));
  f = fs * (0:Nfft/2).'./Nfft;
else
  sfactor = 1 ./ (fs * sum(win.^2));
  f = fs * (-Nfft/2:Nfft/2-1).'./Nfft;
end

% Clip excess data
x = x(1:floor(Ns/(Nw-overlap)) * (Nw-overlap));

% Perform overlap by appending a shifted version
% of the data matrix
y = reshape(x, Nw-overlap, []);
y = [y; circshift(y(1:overlap, :), [0 -1])];
y = bsxfun(@times, y, win);

% get rid of last column
if size(y, 2) > 1
  y = y(:,1:end-1);
end

% Compute the PSD
Px = sfactor * mean( abs(fft(y, Nfft)).^2, 2);

% Return full or onesided result
if oneside == true
  Px = Px(1:end/2+1);
else
  Px = fftshift(Px);
end


Without disagreement to Dilip's valid comments, this is to show the corner case as to when you can detect the presence of DSSS. This is not demodulating the data but detecting the presence of a repeating DSSS signal that is buried well below the noise. Specifically this would occur when the sequence repeats such as done with the GPS C/A code signals (the PRN sequence is repeated 20 times for each data symbol). Such repetition can then be detected by computing an autocorrelation on the received sequence over the duration of repetition, revealing both the presence of a repeating signal and its repetition rate.

The plot below shows the auto-correlation for a GPS signal that is 14 dB below the noise floor, which is a strong but realistic terrestrial GPS signal. Further carrier offsets would need to be removed similar to as done in pre-correlation acquisition search algorithms. The impulse in the center is the impulse you would see if the GPS signal wasn't present (an impulse auto-correlation as expected for white-noise). It is the additional correlation peaks that reveal the presence of something more than white noise.

The time domain signal as received prior to the auto-correlation is plotted below. The GPS signal is only +/-1V in this plot, showing how "buried in noise" it is.

Certainly correlating to the actual sequence would provide a much clearer result as well as the opportunity for data demodulation. This is just showing that the presence of the signal can indeed be detected under this condition of repeating data bits.

In typical DSSS for covert communications, each epoch of the PRN sequence is encoded with an equiprobable data bit making such pattern detection impossible. If one were to send a long series of ones or zeros in the transmitted data prior to spreading, this would be a vulnerability in allowing for the presence of the signal to be detected. Further DSSS on its own, with codes generated from linear code generators (such as C/A code) are not to be considered secure.


Discretizing a Controller with the Backward Difference Method

This is kind of hand-wavy, but you can look at this from two different perspectives:

One, you can look at $z^{-1}$ as a "back-step" operator; i.e. if $X(z) = \mathcal{Z}\lbrace x_n \rbrace$, then (with a few 'i's left undotted and 't's uncrossed) $\frac{X(z)}{z} = \mathcal{Z}\lbrace x_{n-1} \rbrace$.  You can also look at $s$ as a derivative operator: if $X(s) = \mathcal{Z}\lbrace x(t) \rbrace$, then $s X(s) = \mathcal{Z}\lbrace \frac{d}{dt}x(t) \rbrace$.

And -- the time derivative of $x(t)$ is approximately $\frac{d}{dt}x(t) \simeq \frac{x(t) - x(t - T)}{T}$.  If I throw caution to the wind and just sample the right side, that translates to, again roughly, $\frac{d}{dt}x(t) \simeq \frac{x_n - x_{n-1}}{T}$.  Now the Laplace of the left and the Z transform of the right gives $s X \simeq \frac{1 - z^{-1}}{T} X$.  Note that I have left off all pretense of a frequency domain variable on this -- I think that if there's an antithesis of "rigor" it must be "limpor", and we've just found it.

So you can take a system transfer function in $s$, expand it into a differential equation, then you can use the approximation $\frac{d}{dt}x(t) \simeq \frac{x(t) - x(t - T)}{T}$, then you can take the Z transform, then you turn the resulting difference equation into a transfer function in $z$ -- and you have just taken the very long way around to simply substituting $\frac{1 - z^{-1}}{T}$ for $s$.

The other way you can look at this is to observe that with a continuous-time system, the delay operator has a Laplace-domain transform: $\mathcal{L}\lbrace x(t - T)\rbrace = e^{-sT} X(s)$.  So you just set $z^{-1} = e^{-sT}$.  There's ways you can make this more rigorous, but compared to my above treatment this is already granite.

For $sT \ll 1$, $e^{-sT} \simeq 1 - sT$.  Substitute in $z^{-1} = e^{-sT}$ and solve for $s$ and you get $s \simeq \frac{1 - z^{-1}}{T}$.  Voila.  To get the first forward difference start with $z = e^{sT} \simeq 1 + sT$ and you find that $s \simeq \frac{z - 1}{T}$.  Use $z = e^{sT} = \frac{e^{\frac{sT}{2}}}{e^{-\frac{sT}{2}}}$ and you get the Tustin approxmation.

Note, when you actually go to do these, the word approximation.  You can generally get by with control systems because you need to sample way faster than the plant dynamics anyway, but -- approximation.


Choosing Codes or sequences with excellent Auto-Correlation properties

There is a vast literature on this subject. In particular, deleting one bit
(the leading bit in one specific implementation) from each of the Walsh-Hadamard sequences and applying a permutation to the remaining $2^n-1$ bits will result
in $2^n$ sequences of length $2^n-1$ that consist of

(i) the all-zeroes sequence

(ii) the $2^n-1$ cyclic shifts of a pseudonoise (PN), or maximal-length linear  feedback shift register (LFSR), sequence

A PN sequence has ideal autocorrelation properties, but the cross-correlation
between a sequence and its cyclic shift has a peak value that is the same as
the autocorrelation peak.

For more than you probably want to know, I refer you to the paper
D.V. Sarwate and M.B. Pursley, "Cross-correlation properties of pseudorandom and related sequences," Proc. IEEE, vol.68, pp.593-619, May 1980.  It includes a
detailed discussion of the Gold sequences mentioned by JasonR as well as the
small sets and large sets of Kasami sequences that you will encounter in some
applications.


Fixed-point Exponential Moving Average filter with given cutoff frequency

You should scale your analog input to match your desired resolution, it's better to have bigger resolution (specially when dealing with IIR filters), I don't know the reason of returning y>>N, probably the scaling problem is there, I would implemented it like this:

int low_pass_EMA_fixed_point(int x, int y, int alpha_scaled, int Q, int N)
{
  int32_t y_temp;
  y_temp = (alpha_scaled*(x<<N))+((1<<Q)-alpha_scaled)*y;
  return (int)(y_temp >> Q);
}

Also a good practice is test your filter with an impulse to check if it was implemented correctly before using the actual data.

EDIT: tested this program with onlinegdb.com (the printf are used to copy the result into octave):

#include <stdio.h>
#include <stdlib.h>

#define NUM 1024

int low_pass_EMA_fixed_point(int x, int y, int alpha_scaled, int Q, int N);

int main()
{
    int alpha_scaled = 1097; //Q12
    int Q = 12;
    int N = 2;
    static int x[NUM];
    int y = 0;

    x[0] = 1024; //input in Q10
    printf("y=[");
    for(int i=0;i<NUM;i++)
    {
        y = low_pass_EMA_fixed_point(x[i], y, alpha_scaled, Q, N);
        if(i == NUM-1)
        {
            printf("%d];",y); //print the last output
        }
        else
        {
            printf("%d,",y);
        }
    }
    return 0;
}

int low_pass_EMA_fixed_point(int x, int y, int alpha_scaled, int Q, int N)
{
  int32_t y_temp;
  y_temp = (alpha_scaled*(x<<N))+(((1<<Q)-alpha_scaled)*y);

  return (int)(y_temp>>Q);
}

And I'm getting this frequency response in octave:


The obvious answer: define a mathematical model of how the temperature affects your data, and use that model to remove the effect of temperature on your data.

Sorry, but with the amount of information you give on your data, the way your temperature affects it or any approach, this is the precisest possible answer.

I'd like to add that what signal processing defines as "filter" is probably not the tool you're looking for, unless you now the spectral properties of your data and your temperature.


Window functions
Different window functions introduce different amounts of scalloping noise/artifacts. Scalloping noise occurs when the frequency bins of a DFT don't optimally align with the DTFT.

The shape of the DTFT corresponds to the Fourier transform of the window function, and thus window functions whose Fourier transform has shorter side-lobes will have less scalloping noise. However, window functions with less scalloping noise are narrower in the time domain (the Rectangular window has the most scalloping noise) and thus, without window overlapping, will lose the frequencies at the edges of the window.
Terminology
Since the frequency content of overlapping windows are displayed sequentially rather than averaged, I find talking about overlapping windows confusing. I would rather divide the time domain into Frames for which we try to obtain a picture of the frequencies in their range as accurately as possible. This may involve extending the window function past the range of the frames to minimize the scalloping noise/windowing loss tradeoff. However, this will introduce leakage of the frequencies between neighboring frames (reducing time resolution).
Window lengths and sample rates
Edit:
My answer here was bad and I replaced it with this:

Frequency resolution is how similar the window function's magnitude in frequency space is to the Dirac delta function. This is because the product of the window and the signal in the time domain becomes convolution in the frequency domain (and a convolution with the Dirac delta function is a sampling which would give perfect frequency resolution) The fatter the mainlobe (quantified by it's variance), and the higher the sidelobes, the worse the frequency resolution. Additionally, Time resolution can be quantified a the variance of the window function in the time domain.

From wikipedia https://en.wikipedia.org/wiki/Short-time_Fourier_transform#Resolution_issues:

The product of the standard deviation in time and frequency is limited. The boundary of the uncertainty principle (best simultaneous resolution of both) is reached with a Gaussian window function, as the Gaussian minimizes the Fourier uncertainty principle.

The key to note is that Praat's Gaussian Window is twice as long as it's other window functions. In the STFT you can extend you can have an extremely long window function compared to the time interval you want a view of the frequencies for, so long as the window function is close enough to zero for outside of your desired interval. This is optimal. Use the adjustable windows https://en.wikipedia.org/wiki/Window_function#Adjustable_windows
Zero padding
Zero padding doubles the number of bins while keeping the maximum frequency the same, interpolating the DTFT:

From Dan Boschen's answer What happens when N increases in N-point DFT


Line based Extended Kalman Filter for localization : Measurement Jacobian

How are you modeling the state transitions? How are you modeling the output equation i.e. How does the state get to the measurements? Without that information, we can't tell you much about how to get your Jacobian right.

Usually, with multiple measurements, you just stack the state / observations so things go from $2 \times 3$ to $2\times 24$ in the case of 8 measurements, instead of augmenting in another dimension.


What's the best method to extract vocals from music?

The example sounds like pitch tracking combined with comb filtering. Using method like this could be implemented by "drawing" as well, as in SPEAR analysis and re-synthesis program.

If the vocal is in the middle and other instruments are not, then you can also try center channel isolation:

http://www.virtualdub.org/blog/pivot/entry.php?id=102

The methods however, sound horrible in addition to being inaccurate and full of artifacts. Better solution to the problem would be a program that could effectively guess which part of the signal belongs to the voice and which doesn't. Much more difficult for sure though.