I can't see any match with the initial $y_n$ in equation $(1)$, (see update). With a slight change to the notation, you have the following results in the first two cases, for Notation 1 and 2 respectively:

Notation 1:
\begin{align}
y_{n_1}&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=0}^{\frac{N}{2}-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}\tag{$\mathbf{1a}$}\\
&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=0}^{\frac{N}{2}-1}Y_k^* \overbrace{e^{j2\pi n}}^{=1} e^{-j\frac{2\pi nk}{N}}\\
&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\left(\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}\right)^*\\
&=y_{n/2}+y^*_{n/2}\tag{$\scriptstyle{\text{$y_{n/2}=y_n$ for $n=0, \ldots, N/2 -1$}}$}\\
&=2y_{\frac n2}\tag{$\scriptstyle{\text{if $y_n$ is real, then $y_{n/2}$ is real}}$}
\end{align}
Notation 2:
\begin{align}
y_{n_2}&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=1}^{\frac{N}{2}}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}\tag{$\mathbf{1b}$}\\
&=\underbrace{\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+Y_N+\sum_{k=1}^{\frac{N}{2}-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}}_{\textrm{same as notation 1}}-Y_N-Y_{N/2}\\
&=2y_{n/2}-Y_N-Y_{N/2}\\
\end{align}

For Notation 3 however, the equation would be identical to equation $(1)$ if you changed the upper limit of the second summation to $N-1$ instead of $N$, otherwise you have an additional term as a result of the $N^{\rm th}$ term as shown below:
\begin{align}
y_{n_3}&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=\frac{N}{2}}^{N}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}\tag{$\mathbf{1c}$}\\
&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\left(\sum_{k=N/2}^{N-1}Y_k e^{j\frac{2\pi nk}{N}}\right)^*+Y_0\\
&=y_{n/2}+\left(y_{n/2}^{'}\right)^*+Y_0\tag{$\scriptstyle{\text{$y_{n/2}^{'}=y_n$ for $n=N/2, \ldots, N-1$}}$}\\
&=y_{n/2}+y_{n/2}^{'}+Y_0\\
&=y_n + Y_0
\end{align}

So, the correct notation should be:

$$
y_n=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=\frac{N}{2}}^{N-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}\tag{$\mathbf{2}$}
$$

UPDATE:

Regardless of equation $(1)$, to test Notation 1 I will maintain my last notation and go proving if the equality below holds, i.e. if equation $(\mathbf{2})$ equals equation $(\mathbf{1a})$ in Notation 1.

\begin{align}
\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=\frac{N}{2}}^{N-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=0}^{\frac{N}{2}-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}\\
\implies\sum_{k=\frac{N}{2}}^{N-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}&=\sum_{k=0}^{\frac{N}{2}-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}\\
\end{align}

One way could be testing if the following equalities are satisfied:

\begin{align}
Y_{N/2}e^{j\frac{2\pi n(N/2)}{N}}&=Y_Ne^{j\frac{2\pi n(N)}{N}}\\
Y_{N/2-1}e^{j\frac{2\pi n(N/2-1)}{N}}&=Y_{N-1}e^{j\frac{2\pi n(N-1)}{N}}\\
Y_{N/2-2}e^{j\frac{2\pi n(N/2-2)}{N}}&=Y_{N-2}e^{j\frac{2\pi n(N-2)}{N}}\\
\vdots&=\vdots\\
Y_2e^{j\frac{4\pi n}{N}}&=Y_{N/2+2}e^{j\frac{2\pi n(N/2+2)}{N}}\\
Y_1e^{j\frac{2\pi n}{N}}&=Y_{N/2+1}e^{j\frac{2\pi n(N/2+1)}{N}}\\
\end{align}

It's neither the periodicity nor the Hermitian property of the DFT (for a real-valued signal i.e $Y_{N-k}=Y_{-k}=Y_k^*$) that we're seeing in the above equations. But let's check the first two and the last one for instance:

For the first equation:
\begin{align}
&Y_{N/2}e^{j\frac{2\pi n(N/2)}{N}}=e^{j\frac{2\pi n(N/2)}{N}}\left(\sum_{n=0}^{N-1}y[n]e^{-j\frac{2\pi(N/2)n}N}\right)=- \sum_{n=0}^{N-1}y[n]e^{-j\pi n}=\sum_{n=0}^{N-1}y[n]=Y_0\\
&Y_{N}e^{j\frac{2\pi n(N)}{N}}=e^{j\frac{2\pi n(N)}{N}}\left(\sum_{n=0}^{N-1}y[n]e^{-j\frac{2\pi(N)n}N}\right)=\sum_{n=0}^{N-1}y[n]e^{-j2\pi n}=\sum_{n=0}^{N-1}y[n]=Y_0\\
&\implies Y_{N/2}e^{j\frac{2\pi n(N/2)}{N}}= Y_{N}e^{j\frac{2\pi n(N)}{N}}\\
\end{align}
For the second equation:
\begin{align}
&Y_{N/2-1}e^{j\frac{2\pi n(N/2-1)}{N}}=e^{j\frac{2\pi n(N/2-1)}{N}}\left(\sum_{n=0}^{N-1}y[n]e^{-j\frac{2\pi(N/2-1)n}N}\right)= e^{-j\frac{2\pi n}{N}}\sum_{n=0}^{N-1}y[n]e^{j\frac{2\pi n}N}=e^{-j\frac{2\pi n}{N}}Y_1^*\\
&Y_{N-1}e^{j\frac{2\pi n(N-1)}{N}}=e^{j\frac{2\pi n(N-1)}{N}}\left(\sum_{n=0}^{N-1}y[n]e^{-j\frac{2\pi(N-1)n}N}\right)=e^{-j\frac{2\pi n}{N}}  \sum_{n=0}^{N-1}y[n]e^{\frac{j2\pi n}N}=e^{-j\frac{2\pi n}{N}}Y_1^*\\
&\implies Y_{N/2-1}e^{j\frac{2\pi n(N/2-1)}{N}}=Y_{N-1}e^{j\frac{2\pi n(N-1)}{N}}\\
\end{align}
For the last equation:

\begin{align}
&Y_1e^{j\frac{2\pi n}{N}}=e^{j\frac{2\pi n}{N}}Y_1\\
&Y_{N/2+1}e^{j\frac{2\pi n(N/2+1)}{N}}=e^{j\frac{2\pi n(N/2+1)}{N}}\left(\sum_{n=0}^{N-1}y[n]e^{-j\frac{2\pi(N/2+1)n}N}\right)=e^{j\frac{2\pi n}{N}}  \sum_{n=0}^{N-1}y[n]e^{\frac{-j2\pi n}N}=e^{j\frac{2\pi n}{N}}Y_1\\
&\implies Y_1e^{j\frac{2\pi n}{N}}=Y_{N/2+1}e^{j\frac{2\pi n(N/2+1)}{N}}
\end{align}

So, it looks like the initial notation in equation $(\mathbf{1a})$ is indeed identical to equation $(\mathbf{2})$. In conclusion, Notation 1 is correct for a real-valued signal. And Notation 3 need a slight change in the summation limit of the second term to be true. That's my two cents.


Fast Morphological Erosion / Dilation with Variable Window Size

The trick when using masks is that if the window is shifted one to the left, only one row is added and one removed. Therefore one can e.g. calculate the mean of a 5x5 window quicker, as if all 25 values needed to be added and divided again.

In erosion or dilation this fact could be used to, if one would store the max of each row. Then only the maxes of one new row would be needed and a comparison between all rows. However I guess nobody does that, since the overhead is quite big and big windows would need additional space for each row's max...

Your problem however couldn't even use this, since the window size is changing. I assume the changes are random. Therefore there will be no quicker way than moving from pixel to pixel and apply the different window sizes.

But why would someone use different windows for different pixels?


Find filter coefficients to model a device using its measurement

We can represent a digital filter as:

A difference equation, relating input and output at every time step $n$:
$y[n] = \sum_{k=0}^{K-1} b[k] x[n-k] - \sum_{m=1}^{M} a[m] y[n-m]$. The filter coefficients are $b_k$ and $a_m$. In case of an FIR, $a[m] = 0\,\forall\, m $. This form is basically how the filter can be implemented (in the discrete time domain).
By z-Transform of the difference equation, we can also represent the above equation as: $Y(z) = \frac{\sum_{k=0}^{K-1} b[k] z^{-k}}{1 + \sum_{m=1}^{M} a[m] z^{-m}} X(z) = H(z)X(z)$.
The frequency response is simply this z-Transform evaluated at the unit circle $H(z = \mathrm{e}^{j\omega})$ (this is often simply written as  $H(\omega)$).
The impulse response is the inverse transform of $H(\mathrm{e}^{j\omega})$. By convolution with the impulse response, we also can get the output of the filter, as in $y[n] = h[n] * x[n]$. However, $h[n]$ is only finite for an FIR (hence the name, obviously), where we can easily check that in fact it is simply the coefficients $b[k]$.

I hope this helps you with understanding the big picture. So to answer your questions:

The IFFT of the frequency response yields the impulse response. This corresponds directly to the filter coefficients only in the case of an FIR filter. In case you want to model your system as an IIR, you can use an algorithm such as the kind implemented in MATLAB's invfreqz, which is essentially an optimization to fit your frequency response to the complex rational function $H(z)$.
Yes, since $H(z)$ is a representation of the filter's difference equation, the coefficients of which are what we call filter coefficients.

In any case, these are important basics you can learn from e.g. the classic book on DSP. Even the MATLAB documentation (which im sure you have already been checking) is a decent resource.


How is a Particle Filter used to Estimate Parameters of a State Transition Function?

I suppose I can treat the parameters of the state transition function as the space in which I want to generate the particles ... Is this how it's done with a particle filter?

More or less, yes.  You'll find the richest sources of information on this if you look for books on system identification and on adaptive control (of which one important branch is doing system identification on-line, then adapting one's controller to match the identified system).
I can't recommend a good system identification book, but Karl Åström and Björn Wittenmark's "Adaptive Control" is a pretty good book on that topic.


What Is Meant by Clustering of ICA (Independent Component Analysis) Components?

Since ICA is a generalization of PCA (Principal Component Analysis) one could use it for intuition about the process.

In PCA we create a new coordinate system where we can represent each data sample.
So if we have 10 Dimensional data for each sample we have 10 components on the new coordinate system.
Classic clustering approach is to use less than 10 components (Process called Dimensionality Reduction). Then one could apply clustering algorithm of his choice with K-Means Clustering being the classic choice.

The same idea goes for ICA. Once you have a representation of the data in the new system of representation you're applying on it (Or part of it) a clustering algorithm.

Some hands on examples are given by:

Kaggle - Unsupervised Learning: PCA and Clustering.
SciKit Learn - Unsupervised Learning: Seeking Representations of the Data.


1 As defined, the device-dependent fingerprint offset should be constant for a given device so if you have nominal_carrier_freq = 1000; then  you could say

unique_dev_offset = 1;
unique_dev_carrier_freq = nominal_carrier freq + unique_dev_offset

2 Using the Gaussian distribution is common for modeling any deviation from nominal. Think of the context. Our goal is to have devices that are identical, but in real life they have some slight differences. Take a resistor for example, with a nominal value R. If we have a set of resistors we will see some are a little different than R but the majority of them are centered around R. So our model single modal distribution around R. With decreasing probability of being far away from R. This could be a lot of functions (triangle for instance meets those requirements) but the Gaussian is used because it has nice mathematical properties and usually fits real world data very well

3 I think you have to say the channel is constant for this to work. now imagine we have two signals. both are supposed to be the same but they have some unique characteristics. in the images below, I have two signals, one at nominal frequency (0 carrier freq offset), and one with a nonzero carrier frequency offset. We can see that even after adding noise the signals are still distinguishable from each other. Noise won't affect the characteristic frequency offset of the signals. It may modify the appearance of the signals, but not the underlying characteristics

hope that helps


Identifying the age of a speaker using speech recognition

Have a look at the acoustic scale website for some inspiration

http://www.acousticscale.org/wiki/index.php/Main_Page

You can estimate age based on a combination of vocal tract length (vtl) and pitch. Both of these attributes can be extracted from vowel sounds. Children will have short vtls and high pitch, adult males will have low pitch and long vtl, adult women will have fairly long vtls and higher pitch. Dwarfs can have short vtls and low pitch. The castrato will have long vtl and high pitch.


Yes, we use windows to reduce spectral leakage and you show an example where using a window increases it.

There is more than one way to explain it.

If you express a DFT as a matrix vector product

$$ \mathbf{y}=\mathbf{Wx}
$$
where
$$
\mathbf{WW^H}=c\mathbf{I}
$$
where $c$ is a constant that depends on how you define your DFT. The point is that $\mathbf{W}$ is an orthogonal matrix. Your 3 original sine waves correspond to 3 rows of $\mathbf{W}$. Your 3 sines are perpendicular to all the other $N-3$ rows, so there is no leakage.

When you uses a window like a Hamming Window, it’s equivalent to multiplying each row element by element by the window. This modified matrix $\mathbf{\tilde{W}}$ isn’t orthogonal any more.
$$
\mathbf{\tilde{W}\tilde{W}^H}\ne c \mathbf{I}
$$
The resulting matrix is nearly equal to a constant times $\mathbf{I}$ but the off diagonals are nonzero.  Your 3 original sines now project on all the rows of $\mathbf{\tilde{W}}$. These nonzero projections are another way to say leakage.

The nonorthogonality of the widowed matrix is easy to verify in matlab using the dftmtx command.


Generating Kawase Blur Kernels to Approximate a Gaussian Blur on an Image

Looking at Intel - An Investigation of Fast Real Time GPU Based Image Blur Algorithms By Filip Strugar it looks like the Kawase kernel is just a way of implementing a linear kernel quickly, but in a way that constrains the kernel somewhat.

This means that you could make such an algorithm.  Either choose a set of spreads and adjust their weights (if that is possible) for a best fit, or search a number of possible spreads (with the best weights) and find the ones that fit best.

Given the ad-hoc nature of the kernel I suspect you'll end up with a table of best fits for each set of spreads, and then you (or anyone using the table) would want to make a judgement call about effect vs. clock ticks.

The kernel is designed with some specific hardware in mind, so you'll need to check to see if you can adjust the weights at all, or if you can only adjust the spacing and the repetition at each spacing.  At any rate, because the spacing is discrete, its a fairly highly nonlinear optimization problem.  It may be best to just do a brute-force search at each size of Gaussian you want to emulate.  Just try out a bunch of different spacing combinations, and see which ones give you the closest match to a Gaussian of a certain size for the number of iterations.


How does shifting a signal in the time domain affect its frequency domain?

According to the Fourier's transform properties, the Fourier's transform of $f(t+dt)$ would be $F(f)e(j2\pi fdt)$
So, basically, the spectrum (Fourier's transform magnitude, representing the frequency content of your signal) does not change. However, f Fourier's now has a phase-shift proportional to the frequency.

To answer your question : the phase-shift's sign is the time delay's sign.


comparison of speech signal processing and financial data

Financial market trading data is thought to be fractal in dimension and thus not sufficiently band-limited before sampling.  So smaller windows might not have any effect on improving pseudo-stationarity.

Whereas the physical processes behind speech production are better understood and more physically limited (the rate of human muscle contractions controlling lung, larynx, and head cavity resonances, etc.), as well as the sound capture being low-pass filtered below typical human hearing thresholds before sampling above the Nyquist rate.


solving recurrence relation with Z-transform (LTI and causal)

With such simple recursions it's often easiest (and instructive) to actually go through the recursion and find a pattern. With $y[0]=20$ and using $c=35$ we get

$$\begin{align}y[1]&=c+\frac{y[0]}{2}\\y[2]&=c+\frac{c}{2}+\frac{y[0]}{4}\\y[3]&=c+\frac{c}{2}+\frac{c}{4}+\frac{y[0]}{8}\\\vdots\end{align}$$

So we can conclude

$$\begin{align}y[n]&=c\left[1+\frac12+\frac14+\ldots+\left(\frac12\right)^{n-1}\right]+y[0]\left(\frac12\right)^n\\&=c\sum_{k=0}^{n-1}\left(\frac12\right)^k+y[0]\left(\frac12\right)^n\\&=c\cdot\frac{1-\left(\frac12\right)^{n}}{1-\frac12}+y[0]\left(\frac12\right)^n\\&=2c+\big(y[0]-2c\big)\left(\frac12\right)^n,\qquad n>0\end{align}$$

which is the same expression you obtained.


Oversampling then filtering and downsampling

As you've figured out, for the filter specifications you gave, you need a very high filter. To summarize, you're asking for:

Passband edge: 30 Hz
Stopband edge: 32 Hz
Minimum stopband attenuation: 68 dB

The main determining factor for the filter order that you need to meet the specifications is the width of the transition band as a fraction of the sample rate. In your case, you want a transition band that is 2 Hz wide, at a sample rate of 4096 Hz. That's a very sharp transition!
A better approach is to perform the decimation in multiple filter stages that are cascaded in series. This allows you to meet the overall frequency response requirements that you have without using a single very-high-order filter. You need an overall decimation factor of 64 to be implemented by the cascade of filters. It can be hard to determine the optimal (in terms of total overall computations to apply the filters) arrangement of stages, but I would recommend the following approach as a start:

Stage 1: Design a filter with a passband edge at 30 Hz and a stopband edge of (512 - 32) = 480 Hz. Decimate the output of the filter by 8. Thus, at the output of this filter, you will have a signal sampled at 512 Hz.
You might notice that this approach yields some aliasing, as the transition band of the stage-1 filter is allowed to extend beyond the new Nyquist frequency of 256 Hz. However, this isn't a problem, because none of the aliasing occurs within the desired passband of the overall cascade. You're allowing some energy to alias through into what will become the stopband of stage 2, where it will be eliminated. Stretching out the transition band of the stage 1 filter allows you to save some computations, because the wider transition will require a lower filter order. Since stage 1 operates at the highest sample rate, this can be a significant benefit.

Stage 2: Design a filter with a passband edge at 30 Hz and a stopband edge of 32 Hz. Decimate the output of the filter by 8. Thus, at the output of this filter, you will have a signal sampled at 64 Hz with a 30 Hz passband.
While this is still quite a sharp transition band, it should be easier to design because the transition band width is a factor of 8 larger when expressed as a ratio of the filter's sample rate. In addition, the work required to apply the stage-2 filter is reduced further because it operates at the lower 512 Hz rate instead of 4096 Hz.

Using MATLAB's fdatool, the above specifications yield the following filters:

Stage 1: 29th-order FIR, 0.1 dB passband ripple, 80 dB stopband attenuation

Stage 2: 852nd-order FIR, 0.1 dB passband ripple, 80 dB stopband attenuation

I just used the default attenuation settings for illustration. You'll notice that with this setup, the second-stage filter is still quite complex. While I've used filters of this length before, it's likely that you might want to break the problem down even further to avoid such long filters. To do so, you would add more stages of decimation. For instance, you might try three stages of decimation by 4, or 6 stages of decimation by 2. The best way forward is dictated by your system's constraints, but you should see the idea of the approach.


Kalman filter, defining the measurement model

Your observation are the $y_k$ your observation model is defined by $Y_k=HX_k+\eta_k$, where $\eta_k$ if a sample from a normal distribution.

If you want your kalman filter to estimate position, velocity and acceleration your state vector is : $X_k=[x_k,\dot{x}_k,\ddot{x}_k]^T$ and therefore your observation matrix is : $H=[0,0,1]$.

As stated by Rhei, your observations are what you get from the system at each time step, i.e. the measured acceleration.


Can we determine if a filter is butterworth or chebyshev from its physical topology?

As pointed out in a comment by Jason R, you can obviously come up with some random coefficients for the numerator and denominator polynomials of the filter's transfer function, and you'll end up with arbitrary filters that do not satisfy the definitions of any of the filter types you mentioned (even if you make sure that the filter is causal and stable).

Concerning your second question, the answer is also no. A Butterworth and a Chebyshev (type I) filter have the same structure, just the element values are different. The structure of an elliptic filter is different because it has zeros at finite frequencies. For passive filters, these are achieved either by parallel circuits of L and C in the signal path, or by series circuits of L and C towards ground. However, this is also true for type II Chebyshev filters, which also have zeros at finite frequencies. So if your passive filter topology has LC resonators (parallel in the signal path, or serial towards ground) you know that the filter can't have a Butterworth, Bessel, or type I Chebyshev characteristic. It could be an elliptic or a Chebyshev (type II) filter, but it needn't be. If the structure has no LC resonators, the opposite is true: it can't be an elliptic or a Chebyshev (type II) filter, but it could be a Bessel, Butterworth, or Chebyshev (type I) filter. But again, it could also be some non-standard filter. It all depends on the element values.


How to plot magnitude and phase response by hand if I have the Transfer Function?

Use the transformation $z = e^{j\omega}$, you will get

$$H(e^{j\omega}) = \frac{1-e^{-j\omega}}{5\left(1+2e^{-j\omega}\right)}$$

Solving this (using $e^{j\omega} = \cos \omega + j \sin \omega$), you should get something like (please double check):

$$H(\omega) = \frac{\left(\cos\omega-1 \right) + j4\sin\omega}{5\left(5+4\cos\omega\right)}$$

Here,

$$\text{Re}\left(H(\omega)\right) = A = \frac{\left(\cos\omega-1 \right)}{5\left(5+4\cos\omega\right)}$$

and, $$\text{Im}\left(H(\omega)\right) = B = \frac{4\sin\omega}{5\left(5+4\cos\omega\right)}$$

Now, the magnitude response will be

$$\left|H(\omega)\right| = \sqrt{A^2+B^2}$$

and the phase response will be

$$\angle{H(\omega)} = \tan^{-1}\frac{B}{A}$$

You will get both responses as a function of $\omega$, just vary $\omega$, calculate the values and plot the response.


How do I resample an image to a rotated grid?

You should read about image transformations. Anyhow, the idea is quite simple:

Compute the inverse transform $T(x,y)$.
For each point $(x,y)$ on the re-sampled grid:

Calculate $T(x,y)$, and find the 4 closest points on the original grid.
Apply interpolation technique (For example, bi-linear) on the neighbors, and put the value in the re-sampled grid.

Some references:

http://www.cis.rit.edu/class/simg782/lectures/lecture_02/lec782_05_02.pdf
http://www.imageprocessingbasics.com/geometric-transforms/


Obtaining beat histogram using Tzanetakis algorithm

The firsts steps by Tzanetakis can be described as an Scheirer algorithm(Tempo and beat analysis of acoustic musical signals by Eric Sheirer), you can use DTW, filterbanks or FFT to split your signal in N subbands!

I'll try to help you with what I did in the past, my steps are described below.

I can show a simple example of 6(six) (200, 400, 800, 1600, 3200, 6400) decomposed subbands aka "Scheirer 1998″, I used this sample to test:

After DTW or filterbanks you need extract the envelope of each subband:

Apply Full Wave Rectification
Apply Low Pass Filter
Apply Downsample
Apply Mean Removal

My envelopes from the six subbands above are:

OK now do you have six downsampled envelopes (Downsample can help you gain speed to compute the autocorrelation), you need sum all envelopes like:

Almost there, you need apply autocorrelation and to find bpm values between 60 and 200, you need compute the START and END position to be computeded from your summed envelopes:

start = (60 * downsampled_value) / (200)

end = (60 * downsampled_value) / (60)

Now you know the start and end to do the autocorrelation :-)

The next step is find the peak from your autocorrelation and convert to BPM:

BPM=60*downsampled_value/PEAK

And finally ride a histogram with the results of BPM's from each audio piece analyzed !

PS: Tzanetakis and Sheirer algorithm is not perfect its works well in strong rhythmic songs.


Measuring the similarity of two waveforms

It seems the term (and process) I am looking for is "autocorrelation".


Filters using equiripple approximations (different from type I lowpass)

The only place where I found the answer to this question was in Theory and Application of Digital Signal Processing by Rabiner & Gold.

In its chapter about theory and approximation of FIR filters, there is a section where all of my questions are explicitely answered:


Which image transformations are considered continuous transformations?

I think "continuous" is meant in the strict mathematical sense here, i.e. no abrupt changes in the coordinate mapping.

Most non-rigid image alignment algorithms (try to) produce continuous mappings. Distortion filters in popular photo-editing applications (like Photoshop's Warp effect) are usually continuous mappings, too.

ADD: Here are two examples of transformations that are continuous but not affine or perspective:

How to flatten image label on a food jar

Flatten a sudoku grid


How to get the conjugate symmetry of the complex matrix using MATLAB?

This question is off-topic, but:

At = A';

At is the conjugate transpose of A.


Finding the equivalent filter H(u,v) in the frequency domain of a 3x3 spatial mask

Hint: $g(x,y)$ averages all pixel values around  $f(x,y)$, at 8-neighbour connectivity distance 1, except $g(x,y)$. The classical uniform $3\times 3$ average is:
$$ m(x,y) =  \frac{1}{9}\sum_{i\in[-1,0,1]}\sum_{j\in[-1,0,1]}f(x-i,y-j)\,.$$
Thus: $$9m(x,y) = 8g(x,y)+f(x,y)\,.$$
Knowing the Fourier transform of the uniform $3\times 3$ mask, and that of the identity, and using linear superposition, you should be on track.


Fourier transform of certain noisy function

After thinking more about this, I am pretty sure that what you really want to do, and what your MATLAB probably in fact is doing, is study the a function like

$\displaystyle f(t)=\int^tdt^\prime (t-t^\prime)\eta(t^\prime)$

with $\eta(t)$ having the correlation property you've already mentioned.  This will give an FT that is definitely frequency dependent and, I think, looking like your MATLAB results.


It's not clear what you're asking, but I think you want to know why the DFT of a signal has extra conditions for it to real.

For the continuous case: the 2 is there because we need a function to be symmetric to have a real Fourier transform.

In the discrete case, the DFT obey the symmtery:

$$
X_{N-k} = X_{-k} = X_k^*
$$
i.e. conjugate symmetry (half the values are just complex conjugates of the other half). Thus $X_0$ and $X_{N/2}$ are real numbers whilst the remainder of the DFT is specified by $N/2-1$ complex numbers (hence the 2 in the second line).

Also, the answers to this question.

offer counterexamples to your first statement - that analytic functions have vanishing negative frequencies.


Finding the transfer function of a discrete signal described by two equations

HINT:

(because it's a homework problem)

Apply the $\mathcal{Z}$-transform to both equations. Use the first equation to express $Q(z)$ in terms of $X(z)$, and plug that into the second equation to obtain $Y(z)$ in terms of $X(z)$. From this you can easily get the transfer function $H(z)$.
Compute the poles and zeros of $H(z)$ and determine the constant $k$ such that the system is stable and/or minimum-phase. You should know what stability and minimum-phase mean in terms of poles and zeros (if not, it should be easy to find that online or in a textbook).


Analytical expression for the eigenvectors of a 3x3 real, symmetric matrix?

Wolfram|Alpha has made such things easy:


MATLAB: Fourier transform vs $\tt lsim$, results don't match

To perform the FFT approach correctly, you will need the FFT length to be at least $N_U + N_H - 1$ where $N_U$ is the length of $u$ and $N_H$ is the length of $h$.

It looks like $h$ is an IIR filter, so the FFT approach probably won't work well.

The Mathworks has a nice example comparing circular with linear convolution.

IIR filters are best implemented in the time domain. For their specs, they are usually much lower order than FIR filters and so we can afford to implement them directly.


Why does Simulink generate this code for a PID controller?

I'm not sur exactly where you went wrong because your haven't fully explained your approach but for the derivate component of the PID, you should convert the C++ lines to this :
$$ D_{out}[n] = N(err[n]D  -F[n-1]) \\
F[n] = F[n-1] + T_s*D_{out}[n] \\\\$$
In the z-domain : $$
F(1-z^{-1}) = T_s*D_{out} \\
F = \frac{T_s*D_{out}}{1-z^{-1}}\\
D_{out} = N(err * D-Fz^{-1}) \\
D_{out} = N(err*D-\frac{T_sD_{out}}{1-z^{-1}}z^{-1}) \\
D_{out}(1+\frac{NT_sz^{-1}}{1-z^{-1}}) = N(err*D) \\
\frac{D_{out}}{err} = \frac{N*D}{1+\frac{NT_sz^{-1}}{1-z^{-1}}} \\
\frac{D_{out}}{err} = \frac{N*D}{1+\frac{NT_s}{z-1}}
$$
Edit :
You posted C++ code. Why is x_last only updated once? It needs to be updated every period. Try my solution below
double y0=0.0;
double Ts=0.005;
double PID_DifferentialPart(double x,double D,double N){
  static double x_last,y_last;
  static int first=1;
  double y;
  if(first){
    first=0;
    x_last = 0;
    y=y0;
  } else {
    y=D*N*(x-x_last)-(N*Ts-1.0)*y_last;
  }
  y_last=y;
  x_last=x;
  return y;
}


How to segment multiple overlapping coins (ellipses)

Segmentation is generally a process that is very susceptible to noise. I would better use a detector, especially for geometric shapes like coins. Remember, if you have a good detection, you also ease the segmentation problem dramatically. For the example of coins, a good model would be to use an ellipse: every circle/ellipse appears to be an ellipse under perspective projection, up to a negligible error.

Luckily, there are many good methods that can detect ellipses without segmentation. ELSD is one, Hough transform is another. In a previous post, I have explained how to use ELSD to spot spherical particles. Now I will take the Hough transform approach - this results in a minimal code that works and which I will share below. This might also be more intuitive to understand.

I will now explain a hypothesize-merge algorithm. Assume that we are given a grayscale image $I$, with gaussian smoothed edges $E$. The idea would be to apply a Hough ellipse detector on $E$ to get a large set of ellipse hypotheses $\{\mathbf{\Psi}_{1..n}\}$ and then cluster them to a reduced subset $\{\mathbf{\Psi}_{1..k}\}$ where $k<<n$. If $k$ is known (for instance this image with six coins has $k=6$) we can make use of this knowledge. Though, I will not assume this. Using the referred ellipse detector (see above) and an available mean shift clustering, I have devised the following MATLAB script (all of this can be implemented in OpenCV with reasonable little effort):

%% PARAMETERS
% blur the image to suppress some noise
gaussSigma = 1.0;

% parameters of the canny algorithm
edgeThMin = 0.03;
edgeThMax = 0.3; % usually set to 10*edgeThMin

% override some default parameters - tune them to the size of a coin
minEllipseMajorAxis = 45;
maxEllipseMajorAxis = 110;

% coins are circle like but in case of projectivity set to a lower value
minEllipseAspectRatio = 0.7;

% bandwidth of the meanshift filter. this is a quite critical parameter
% that determines the number of detected ellipses and their grouping.
meanShiftBandwidth = 25;

%%
% some more parameters do exist within the code below, but
% they should be somewhat less important.
I=imread('XqEYO.jpg');

% smooth the image a bit - increases reliability
Ig = imgaussfilt(I, gaussSigma);
G = rgb2gray(Ig);
E = edge(G,'canny', [edgeThMin, edgeThMax]);

params.minMajorAxis = minEllipseMajorAxis;
params.maxMajorAxis = maxEllipseMajorAxis;

params.minAspectRatio = minEllipseAspectRatio;

% there can be multiple correct hypotheses (filter later)
params.numBest = 200;

% pair with k other points N^2 -> k*N
% less randomization means more speed
params.randomize = 5;

% note that the edge (or gradient) image is used
% code from: https://www.mathworks.com/matlabcentral/fileexchange/33970-ellipse-detection-using-1d-hough-transform
bestFitHyp = ellipseDetection(E, params);

% a fit looks like: (ra,rb,ang,x0,y0,C,Nb)
fprintf('Output %d best fits.\n', size(bestFitHyp,1));
%
figure; image(I);
%ellipse drawing implementation: http://www.mathworks.com/matlabcentral/fileexchange/289
ellipse(bestFitHyp(:,3),bestFitHyp(:,4),bestFitHyp(:,5)*pi/180,bestFitHyp(:,1),bestFitHyp(:,2),'g');
title('Hough Ellipse Hypotheses');

% now cluster the ellipses that belong together using x0,y0 only.
% we use mean-shift clustering : https://de.mathworks.com/matlabcentral/fileexchange/10161-mean-shift-clustering
[~,idx,~] = MeanShiftCluster(bestFitHyp(:,1:2)',meanShiftBandwidth, 0);
numClusters = max(idx);
bestFits = zeros(numClusters, 6);
for i=1:numClusters % for each class
    ind = find(idx==i);
    members = bestFitHyp(ind, :);
    for j=1:6 % for each component of an ellipse record the median
        bestFits (i, j) = median(members(:,j));
    end
end

figure; image(I);
%ellipse drawing implementation: http://www.mathworks.com/matlabcentral/fileexchange/289
ellipse(bestFits(:,3),bestFits(:,4),bestFits(:,5)*pi/180,bestFits(:,1),bestFits(:,2),'g');
title('Clustered Hypotheses');

To determine the cluster centers I use a median as this might be more robust than the mean alternative - you can try different things though. Let's look at the outputs of different steps.

The edges of the smoothed image:

Initially we allow for a large number of hypotheses to survive (~120). This initial result, together with the clustered hypotheses (final detections) look like:

Of course this result is a rough detection. In other words, it is an outcome of voting and not an outcome of non-linear parameter refinement. To further improve these results, you can write an optimizer to snap the ellipses on top of the image edges. See here on descriptions of how to do that. You can also use ELSD edges instead of the Canny to ease the job of the Hough transform and to make the algorithm more robust.

It might not be the best but I hope this gives you a direction. Here are more sample outputs of running the method on coin photos fetched from Google Images:


Signals Cross Correlation from Their Power Spectral Densities?

No,

For instance 2 instances of AWGN with the same STD have the same PSD while being totally uncorrelated (Moreover, Statistically Independent).

Also, pay attention that Correlation is not commutative:

$$ f \ast g \neq g \ast f \Leftrightarrow F \left( k \right) \overline{G \left( k \right)} \neq G \left( k \right) \overline{F \left( k \right)} $$


How to perform carrier phase recovery in software?

To demodulate a phase-shift keyed signal, of which BPSK is the simplest, you have to recover the carrier frequency, phase, and symbol timing.

Bursty Signals
Some signals are bursty and provide a known data sequence called a preamble or mid-amble (depending on whether it shows up at the beginning or middle of the burst).  Demodulators can use a matched filter that "looks" for the known data sequence and use it to figure out the burst's frequency, phase, and symbol timing.  They do this for every burst and don't need to bother with "drift" (the gradual distance that builds up between the demodulator and the received signal as small errors accumulate) because the bursts are generally short enough that the drift isn't a problem.

Continuous Signals
Then there are continuous signals.  They are in many ways tougher to recover than bursty signals because you normally do not have a known data sequence to aid in locking up with the signal, and you do have to worry about drift even after you lock onto the signal.  I'll try to describe at a high level the main steps that are typically used to recover continuous signals.

Carrier Recovery

Usually you know what frequency the signal you are looking for will be at, or at least what set of frequencies it could be at.  Even with this knowledge, though, you usually have to be able to correct frequency offset because no two transmitters transmit at the same frequency.  There is always some error.  The usual method, then, is to mix the frequency that you think that the signal will be at, and then correct the residual frequency error.  This can be done with a Costas Loop, or by taking the fourth power of the baseband signal data and looking for a frequency spike.  There should be a frequency spike at the carrier offset * 4 (e.g. if you take an FFT of the data to the 4th power and see a frequency spike at 8300 Hz, that means the carrier offset is 8300/4 = 2075 Hz).  This is a very effective means of getting an initial lock on the frequency offset.  You can also use it to compensate for drift if you re-do it occasionally.  There is another way to compensate for drift which I will touch upon later.

Carrier Phase

At this point if you plotted your complex data in the complex plane (x-axis is real, y-axis is imaginary) it should look something like the following-

If you look closely you can see two dense areas towards the ends of the blurry line.  Those are the BPSK constellation points.  The points in between are the transitions between the constellation points.  Those will clear up once we get the symbol timing.  The reason that the line is at an angle is because of the carrier phase.  This can be measured by mirroring the signal by multiplying all points that have negative real values by $e^{j*\pi}$

and then taking the mean of the angles of the points.  Once you have calculated that, subtract that angle from all points that have had the carrier offset removed by multiplying the points by $e^{j*-\omega}$.  You can also compensate for carrier offset drift with this technique by dynamically updating the phase offset.  Once the data is phase corrected it should look something like this-

Once the data is phase corrected you can drop the imaginary part of the data because it does not add any information.

Symbol Timing

You should normally know, a priori, the symbol period of the signal you are trying to demodulate.  If you do need to determine the symbol period/frequency though, you can do it in a similar way to how the carrier offset was detected.  You can square the data which will cause a frequency spike at twice the symbol frequency.

Like with the carrier you will have to get the phase (timing) right and then compensate for drift.  The usual method for both of these problems is to look for the zero crossings.  Unless the noise is quite bad it should only cross the zero-point in the middle of a symbol transition from -1 to 1 or from 1 to -1.  Even if noise does cause this to happen in the middle of a symbol it won't happen very often.

The picture above is usually called an "eye diagram" or "eye pattern".  It is two symbol periods wide, and has many symbols "stacked" on top of one another.  I don't know if you are familiar with oscilloscopes or not, but you can get an oscilloscope to show a picture like this.  Anyway, the two "X's" are symbol transitions.  The high to low lines in the X's are when the symbol transitions from a 1 to a -1, and the low to high lines in the X's are when the symbol transitions from a -1 to a 1.  The point in between, where the pink line is, is the optimum place to sample the data to see if the symbol is a 1 or -1.

This same technique can be used to handle symbol timing drift.  Take a running average of the distance from the prior zero crossing and following zero crossing.  If the two averages are about the same, all is well.  If one is bigger than the other then you need to shift where you are grabbing your sample.

Once you have sampled the symbols at the correct points your data points should get something that looks like the classic BPSK constellation points.

I hope this helps.


Zero peaks in fft data even after subtracting mean

You must be assuming you have a dominant peak at zero from only looking at the same graph you shared, but if you really did remove the mean, then value at bin 0 will be 0 (as bin 0 is directly proportional to the mean). Inspect the data carefully as their does not appear to be anything wrong with the code with respect to the FFT.

What is likely occurring is the bins close to bin 0 are non-zero, and these represent the lowest frequency fluctuations in your data. Each of these is zero-mean, but represent the slowest changing frequency components. In this case you are getting the exact results expected, representing what is in your actual data.


My proposal would be:

Filter noise
Create image without the lines
Subtract the image without the lines from the image with the lines
Apply threshold
Postprocess the image

Code:

close all;
% acquire input image
img = imread('ridges.bmp');
img = rgb2gray(img);
% remove some noise
img_gauss = imfilter(img,fspecial('gaussian',10,2));
figure();
imshow(img_gauss);
% remove all the lines from the image
img_filtered = medfilt2(img_gauss,[31,31],'symmetric');
figure();
imshow(img_filtered);
% subtract image without lines from image with lines
img_subtracted = uint8(double(img_gauss)-double(img_filtered)+127);
figure();
imshow(img_subtracted);
% threshold
img_thresh = img_subtracted>130;
figure();
imshow(img_thresh);
% use open to remove some noise
figure();
img_opened = imopen(img_thresh,ones(3));
imshow(img_opened);


Based on your observation, that spectrum of the signal is sufficiently distinguishable, you can use this as a feature to classify laughter from speech.

There are many ways you can look at the problem.

Approach #1

In once case, you can just look at vector of the MFCC. and apply this to any classifier. Since you have many co-efficient in frequency domain, you might want to look at Cascade Classifiers structure with boosting algorithms such as Adaboost based on this, you can compare between speech class vs. laugh class.

Approach #2

You realise that you speech is essentially a time varying signal. So one of the effective ways to do it is to look at time variation of the signal itself. For this, you can split signals in batches of samples, and look at the spectrum for that time. Now, you may realise that laugh might have more repetitive pattern for a stipulated duration where as speech inherently posses more information and hence the spectrum variation would be rather bigger. You can apply this to HMM type of model to see if you consistently remain within the same state for some frequency spectrum or you continuously keep changing. Here, even if occasionally the spectrum of speech resembles to that of laughter it will be more time changing.

Approach #3

Force to apply LPC/CELP type of coding on the signal and observe the residue. CELP Coding makes a very accurate model of speech production.

From the reference here: THEORY OF CELP CODING

  The redundancies in the speech signal are almost removed after the
  short term prediction and long term prediction of the speech signal
  and the residual has very little correlation left in it. Then an
  excitation is searched which synthesizes the speech and the codebook
  index and gain are searched from the fixed codebook. The optimum
  codebook index selection criterion is based on MMSE between the
  locally synthesized speech and the original speech signal.

To put it simply, after all the speech which is predicted from analyser is removed - what is left is the residue which is transmitted to recreate exact waveform.

How does that help with your problem? Basically, if you apply CELP coding, speech in the signal is mostly removed, what remains is residue. In case of laughter a majority of the signal might be retained because CELP will fail to predict such a signal with vocal tract modelling, where as individual speech will have very little residue. You can also analyse this residue back in frequency domain, to see if it is laughter or speech.


MP3 is a lossy compression which exploits -psychoacoustic masking.

Basically, the human ear has a perception limitations in identifying sounds due to active processing as well as due to non-linear frequency response.
Also note, that sound pressure level - the scale that corresponds to perceived loudness is actually logarithmic not linear.

OK - so this only indicates that we should "chop -off" frequencies that are not perceptible to human ear and keep those for which human ear has a high response. this simply can be achieved in 2 ways - either simply apply a filter to the sound and then encode it, or convert the entire signal into frequency domain and drop those frequencies that are on the edge.

NOW: MP3 does something much more than this.

Refer to the image here (Credit wikipedia):

What we see here, is that in the presence of frequency 250Hz, the 175Hz sound is not-audible at all to the human being unless that frequency also has dB level above a certain threshold mask. This frequency is called masking frequency. This is not a dominance of high frequency or low or otherwise. This is simply a function of which frequency produces sound at higher dB (at that time). the threshold of the audibility is also not static but depends on the dB level of the masking frequency. So in order to identify which all frequencies will go un-audible in the presence of other masking frequency one needs to first analyze the signal.

In actual effect the frequency masking effect is continuous in nature - but in practical codec - this is applied over bands.

Hence, we first divide the signal in frequency bands and then identify if any of the frequency bands will go un-audible in the presence of frequency bands which are louder. This is impossible to achieve through linear filtering of any kind! Hence, the filter-bank stage is employed.

UPDATE:

The different implementations of filter banks:
MP3 uses 2 stage filter bank - PQMF which produces 32 sub-bands followed by MDCT filter-bank which produces 575 bands. This is called hybrid filterbank. In general, MDCT can achieve higher resolution of filter bank which - in general PQMF can't achieve. Hence, Advanced Audio Coded (codec usually found in m4a and mp4 files) actually removes the redundant PQMF stage and uses MDCT event at a higher resolution. See Ref 3.

However, so the question is why is MP3 stuck to 2 stage filter bank?
As it turns out, the reason is more to do with compatibility with MPEG Layer I and Layer II. As indicated in the document here MP3 AND AAC EXPLAINED Ref 4
right from Fraunhofer, it says the real reason why MP3 remain stuck with 2 stage filter bank is due to backward compatibility with MPG 1 Layer I & II. From the paper (3.3.1):

  The polyphase filterbank has the purpose of making Layer-3 more
  similar to Layer-1 and Layer-2.

Ref [1]: Wikipedia - Psychoacoustics
Ref [2]: Theory behind MP3
Ref [3]: Modified Discrete Cosine Transform-Its Implications for Audio Coding and Error Concealment
Ref [4]: MP3 AND AAC EXPLAINED also available here.


Pitch detection using Cepstrum in matlab

I wrote a very very basic cepstrum pitch track in the past, well we can try using a pure senoidal signal to show how its works... let me see

For Cepstrum I always have used this steps:

Apply hamming windows in the signal
Apply FFT
Get magnitude
Convert to log scale
Apply IFFT

The equation for cepstrum:

IFFT(log(abs(FFT(s))))

But you can use FFT or IFFT, take a look:

IFFT(log(abs(FFT(s)))) == real(FFT(log(abs(FFT(s)))))

The difference is the scale representation, if do you end using FFT you need extract just the real information, for both above equations you will get the same shape:

For IFFT(log(abs(FFT(s)))):

For real(FFT(log(abs(FFT(s))))):

This is a cepstrum shape example from 4096 points sine in 440hz sampled at 44100hz

Now try get the maximum from the half cepstrum data to find the periodicity...

see my extremely basic code of how to get a sine frequency, just to exemplify to you:

Fs=44100;

%creating one simple signal

f=440; %we need find this
frame = sin(2*pi*f/Fs*(1:2048*2))';

win = hamming(length(frame));

%samples multplied by hamming window
windowedSignal = frame.*win;

fftResult=log(abs(fft(windowedSignal)));

ceps=real(ifft(fftResult));

nceps=length(ceps)

%find the peaks in ceps

peaks = zeros(nceps,1);

k=3;

while(k <= nceps/2 - 1)
   y1 = ceps(k - 1);
   y2 = ceps(k);
   y3 = ceps(k + 1);
   if (y2 > y1 && y2 >= y3)
      peaks(k)=ceps(k);
   end
k=k+1;
end

%get the maximum ...
[maxivalue, maxi]=max(peaks);

result = Fs/(maxi+1)

my result was close 441hz for better results maybe do you need apply parabolic interpolation ....


Result of conjugate symmetry property of DFT

The textbook probably even has the proof; but it's not hard to do: Simply write down the sum that gives you the $k$th bin of the transforf $X[k]$ of the time series $x[n]$. Then, flip $k$'s sign and apply the conjugate to finish your proof.

To show:

\begin{align}
X[k] &\overset{\text{claim}}= X^*[-k]\\
X[k] &= \sum_{n=0}^{N-1}x[n]e^{-j2\pi k\frac nN}\\
X^*[-k] &= \left(\sum_{n=0}^{N-1}x[n]e^{-j2\pi (-k)\frac nN}\right)^*\\
&= \sum_{n=0}^{N-1}\left(x[n]e^{-j2\pi (-k)\frac nN}\right)^*&\text{Sum operator linear}\\
&= \sum_{n=0}^{N-1}x[n]\left(e^{-j2\pi (-k)\frac nN}\right)^*&x[n]=x^*[n]\text{ since real}\\
&= \sum_{n=0}^{N-1}x[n]e^{-(-j2\pi (-k)\frac nN)}\\
&= \sum_{n=0}^{N-1}x[n]e^{-j2\pi k\frac nN}\\
&=X[k] &\rule{1ex}{1ex}
\end{align}

Now, having showed that half of the coefficients are just the conjugate of the others: It's obviously sufficient to only consider one half; the other can be directly calculated from that.


Reconstructing Tree Limbs that are Occluded by Leaves

The masked area hidden by leaves is very large and while the border and shading seems to imply what is hidden to a human being it is unsuited to "inpainting" which relies on pretty consistent boundary, such as where a foreground object temporarily obscures something continuing behind.

The branch orientation is a start but not enough data on its own to imply how the twigs extend without having some external concept of an oak tree versus a willow, poplar or pine.  That is more about 3d fractal modeling than about image restoration.

The least obviously visible algorithms for in-painting involve texture synthesis called "Healing" in Photoshop.  G'MIC and Resynth are popular open-source tools for that available as plug-ins for GIMP.  While visual continuity is quite good they have no concept of a tree or depth, just 2-d patch statistics so would need input from an artist as well in the example given.  The mask is white background on the outer edges and thick branches on the inside .. so texture interpolation gives a weird unnatural result.


Autocorrelation of the product of deterministic and random signal

As correctly pointed out in the comments, in general the process $Y(t)=x(t)M(t)$ is not wide-sense stationary (WSS), i.e. its autocorrelation function depends not only on the time difference parameter $\tau$, but also on the absolute time $t$:

$$R_Y(\tau,t)=E[Y(t+\tau)Y^*(t)]=E[x(t+\tau)M(t+\tau)x^*(t)M^*(t)]=\\
=x(t+\tau)x^*(t)E[M(t+\tau)M^*(t)]=x(t+\tau)x^*(t)R_M(\tau)\tag{1}$$

In your case you just need to evaluate (1) with the given function $x(t)$. However, as expected you will end up with a function of two variables because $Y(t)$ is not WSS.

There are a few special cases in which the resulting process is indeed WSS or can be easily made WSS. The first case is modulation by a complex exponential:

$$x(t)=e^{j\omega_ct}$$

in which case

$$x(t+\tau)x^*(t)=e^{j\omega_c(t+\tau)}e^{-j\omega_ct}=e^{j\omega_c\tau}$$

only depends on $\tau$ and not on $t$. Another case of interest is the case where $x(t)$ is $T$-periodic. In this case the process $x(t)M(t)$ can be made WSS by introducing a random phase epoch $\Theta$ with uniform distribution in the interval $[0,T]$ which is independent of $M(t)$:

$$Y(t)=x(t+\Theta)M(t)$$

The autocorrelation function of $Y(t)$ is then given by

$$R_Y(\tau)=R_M(\tau)\frac{1}{T}\int_0^Tx(\alpha+\tau)x^*(\alpha)d\alpha$$


Very basic question about how we define frequency in signal processing

For the tangent function a cycle is indeed done in $\pi$ radians, because you want

\begin{align}\tan(x) &= \tan(x+T)\\
\implies \frac{\sin(x)}{\cos(x)} &= \frac{\sin(x+T)}{\cos(x+T)}
\end{align}
$$
\implies \sin(x+T)\cos(x) - \sin(x)\cos(x+T) =0
$$
$$\implies \sin\big[\left(x + T\right)-x\big] =0\equiv \sin\left(T\right)=0$$
$$\implies T = \pi k, \quad \text{with}\quad k\in \mathbb Z$$

Smallest strictly positive $k$ (i.e. $k= 1$) gives us $T_0 = \pi$. Anyway, we know that.

Regarding the interpretations, let's use $\Omega \ [\text{rad/sec}]$ to avoid ambiguity with $\omega \ [\text{rad/sample}]$ mostly used in the discrete case.  From the definition of the angular frequency $\Omega = 2\pi F$ with the fundamental period $T_0 = 1/F$, as you have already described, we have
$$
\Omega_0 = \frac{2\pi}{T_0}\tag{1}
$$
The interpretation should go from $(1)$ and read as the number of cycles of the periodic signal in $\mathbf{2\pi}$ radians. As examples, look at the plot below,

As seen above,

$\tan(x)$ has $T_0 = \pi$ implying $\Omega_0 = 2\pi/\pi=2 \ \text{rad/sec}$ or 2 cycles in $\mathbf{2\pi}$ radians (or equivalently one cycle in $\pi$ radians)
$\sin(x)$ has $T_0 = 2\pi$ implying $\Omega_0 = 2\pi/2\pi=1 \ \text{rad/sec}$ or 1 cycle in $\mathbf{2\pi}$ radians
$\sin(2\pi x)$ has $T_0 = 1$ (see $x$-axis) implying $\Omega_0 = 2\pi/1=2\pi \ \text{rad/sec}$ or $\mathbf{2\pi\approx 6.2832}$ (see $x$-axis) cycles in $\mathbf{2\pi}$ radians.

See it this way, the ratio in equation $(1)$ is simply stating how many  $T_0$ (i.e. fundamental periods or full cycles) you have in $2\pi$. You have $T_0 = 2\pi \implies \Omega_0 = 1 \ [\text{rad/sec}]$ meaning you have $1$ fundamental period (full cycle) in $2\pi$. Similarly, when $T_0 = 1 \implies \Omega_0 = 2\pi \ [\text{rad/sec}]$ meaning you have $2\pi$ fundamental periods (full cycles) in $2\pi$. Etc. for other cases.


Confused with filtering linear chirp signal using band-pass filter

You are seeing the difference between "wide band FM" and "narrow band FM" in your modulated waveform. Observe your modulation index for each case and then review the sideband levels versus mod index using Bessel functions. (The modulation index is the ratio of the frequency deviation over the frequency modulation). The first case is wide-band FM as your deviation (approximately +/-500 Hz) is significantly larger than your mod rate (50 Hz), so the mod index is >>1. The second case is narrow band FM as your deviation is significantly smaller than your mod rate (9850) and the mod index is <<1.  Narrow band FM will have only one pair of significant sidebands that will be at the mod rate offset from your carrier; while wideband FM will have a significant number of sidebands each spaced at your modulation rate and going out slightly beyond your deviation before becoming insignificant.

So when you are doing narrow band FM with one pair of sidebands, note that your signal will be dominantly your carrier with two sidebands that are each 9850 Hz away! You are seeing leakage products at your filter output; notice how small the amplitude is compared to the first case.  (Given your modulation is a ramp and not a sine wave, you will actually have multiple sidebands for each of the Fourier series components of your ramp signal, but the point is these will all be spaced by 9850 Hz from your carrier frequency which will also be dominant given the small moduation index).

While in your wideband case you will have spectral energy at spacing every 50 Hz with approximately 1KHz total bandwidth (review Carson's rule on this) and thus will provide tones that will pass through your filter.

Below is a plot of the Bessel function of the first kind from this site:
http://www.statsref.com/HTML/index.html?key_functions.html

This shows the relative level of the carrier (given by $J_0(x)$) and the first five sidebands $J_1(x), J_2(x),J_3(x),J_4(x), J_5(x) $ at a given modulation index x.


How to filter a regular ticking out of a digital (on/off) signal?

Are you processing this data offline or in real time? If you're doing offline processing, or can live with a significant delay, then you could apply a simple filter as follows:

From each input sample, check to see if the sensor reported "on".
If it did, search for a "few" 16-minute periods in each direction and see whether the sensor reported "on" there as well. This is a noncausal operation; if you're processing data in real time, you can search in both directions by first applying a large enough delay to the signal.
If there are "enough" consecutive "on" samples at the spurious period in either direction, suppress the "on" indication from the sensor.

This should suppress the unwanted square wave of "on" indications, but as indicated in the comments, it also could possibly suppress valid "on" reports. The "few" and "enough" parameters you'll have to select yourself based upon what is reasonable given the overall behavior of your system. Basically, you want to set "few" so that you're looking at enough spurious periods to make a good decision. You'll want to set "enough" small enough to suppress short bursts of the spurious period, but large enough to avoid suppressing too many valid "on" indications that happen to be separated by 16 minutes.


the matter is due to this one:

product = ft1[p] * conj(ft2[p])

You can change to:

product = ft2[p] * conj(ft1[p])

I guess input1 is the reference image.


The problem is not the spectrogram parameters, these are correct since they only depend on what resolution you want in time and frequency domain. Also, the spectrogram interpretation is correct, there are multiple frequency peaks. The problem may be:

  I expected to see one high power frequency after pressure rise, instead of multiple frequencies

Why? If you are not in free-field conditions (and you mention a cylinder) there will be multiple wave reflections, resonances, standing waves and 3D pressure patterns. A sudden pressure step will distribute its power over the spectrum, but frequency peaks will appear because you are essentially "tuning" your cylinder.

Recall that pressure waves travel at sound speed. Considering for example the habitual 340 m/s; if your cylinder measures 1m, the wave front will travel back and forth each 0.0029s, adding multiple frequency components to the measurement.

As for standing wave patterns, we can consider the case of a closed cylinder, where we can compute the frequency of the fundamental frequency as $f_1 = \frac{a}{4L}$ with subsequent harmonics at $nf_1$ where $n$ are odd integers (if your cylinder is indeed closed, does it have a length of ~20cm?).

Regarding the problem of measuring higher frequency waves in ducts, consider that according to Eriksson, it may be assumed that waves in a circular duct propagate in a planar fashion (i.e., pressure is constant in every cross-sectional area of the duct) only below a cut-off frequency of:
$$
f = 1.84 \frac{a}{\pi D}\sqrt{1-M^2}
$$

Where a is the sound speed, D the duct diameter and M the Mach number of the mean flow. Above this frequency, the pressure distribution in a given section of the duct will not be uniform, it will have nodes. See the following animation for an idea on how pressure varies in a given planar section, noting that only the Z(0,0) is a typical planar wave.

Beyond the planar wave frequency it is difficult to correctly measure the pressure spectrum with only one transducer. In any case, I believe that the signal processing that you are performing is correct; the only problem is that you are witnessing the acoustic effects of your cylinder.


Computing real signal with minimum absolute values from even magnitude spectrum

Turns out this problem is known as crest factor minimization. I'm using an algorithm by Yang et al. (Phyiol. Meas. 36, 2015). Crest factor minimization algorithms typically are not proven to work on all possible magnitude spectra and only rarely proven to work on magnitude spectra with certain characteristics. Also, some iterative algorithms like Yang et al.'s have been shown in later publications to get stuck at a local minimum, not the global minimum. However, Yang et al.'s algorithm is satisfactory for the magnitude spectra I've tried so far.


I assume this is homework, but transforming a difference equation to the $z$-domain is simple; just recall the time-shifting property of the transform.

$$
x[n] \Leftrightarrow X(z) \rightarrow x[n-k] \Leftrightarrow z^{-k}X(z)
$$

So then we have:

$$
y[n] = \frac{1}{2}x[n] + x[n-1]
$$

$$
Y(z) = \frac{1}{2}X(z) + z^{-1}X(z)
$$

The transfer function can be written as:

$$
H(z) = \frac{Y(z)}{X(z)} = \frac{1}{2} + z^{-1}
$$


How to understand the energy centralizes at low frequency after making DCT on an image?What does energy mean?

It's just a fancy way of saying "amplitude" or "signal strength".  The DCT splits the original signal into a bunch of outputs, each of which gives the amplitude for a particular frequency.  Add them all back together again, through an inverse DCT, and you get the original signal back.


How to validate an estimated model in case of output-only data (in frequency domain)?

I don't know of a standard way to do this, but I think taking a leaf out of the training of artificial intelligence algorithms might help:

divide your output data into two disjoint sets,
perform whatever identification procedure you want to obtain the ARMA parameters, and
then compare the two transfer functions using whatever distance metric makes sense for your application.

In the figure below, I've decided that my distance metric is that the poles and zeros of each ARMA model are "not too far" from each other (or the ARMA poles and zeros found using all the data).

The $\color{black}{\bf black} $ markers are the pole ($\times$) and zero ($\circ$) positions for the ARMA model calculated from first half of the data ($G_1$).
The $\color{red}{\bf red} $ markers are the pole ($\times$) and zero ($\circ$) positions for the ARMA model calculated from second half of the data ($G_2$).
The $\color{grey}{\bf grey} $ markers are the pole ($\times$) and zero ($\circ$) positions for the ARMA model calculated from all of the data ($G_{\rm tot}$).

Another way to do it might be to look at the transfer function
$$
E = G_1 - G_2
$$
That is plotted below using the freqz command in R.

R Code Below

#30332
library(signal)
library(tseries)

# Create a time series that is not just white noise
T <- 1000

do_plot <- 1

#for (k in seq(1,100))
{
  x <- rnorm(T, 0, 1)
  bf <- butter(4,0.41)
  y <- filter(bf$b, bf$a, x)

  # Split data in two and compare ARIMA models
  L2 <- length(y)/2
  y1 <- y[100:L2] # Start away from 1 to avoid startup transients
  y2 <- y[(L2+1):length(y)]

  orders <- c(3,0,3)

  # First half
  tryCatch(
  {arima1 <- arima(y1, order = orders)},
  warning = function(w) {  },
  error = function(e) {  },
  finally = {});

  all_1 <- coefficients(arima1)
  num_1 <- all_1[4:6]
  den_1 <- all_1[1:3]

  # Second half
  tryCatch(
  {arima2 <- arima(y2, order = orders)},
  warning = function(w) {  },
  error = function(e) {  },
  finally = {});

  all_2 <- coefficients(arima2)
  num_2 <- all_2[4:6]
  den_2 <- all_2[1:3]

  # All data
  arima_tot <- arima(y, order = orders)
  all_tot <- coefficients(arima_tot)
  num_tot <- all_tot[4:6]
  den_tot <- all_tot[1:3]

  # Plot  boundaries
  xs <- c(-1,1)
  ys <- xs

  if (do_plot == 1)
  {
    par(mfrow=c(1,1), pty='s')
    plot(roots(den_1), xlim=xs, ylim=xs, lwd = 6, pch=4, xlab="Re(z)", ylab="Im(z)")
    do_plot <- 0
  }
  else
  {
    points(roots(den_1), lwd = 6, pch=4)
  }
  points(roots(den_2), lwd = 4, col="red", pch=4)
  points(roots(den_tot), lwd = 2, col="grey", pch=4)

  points(roots(num_1), lwd = 6, pch=1)
  points(roots(num_2), lwd = 4, col="red", pch=1)
  points(roots(num_tot), lwd = 2, col="grey", pch=1)
}

par(pty="m")
subtract_num <- conv(num_1, den_2) - conv(num_2, den_1)
subtract_den <- conv(den_1, den_2)

freqz(subtract_num, subtract_den)


I don't fully understand your question, since you correctly state

  $\tau_i$ are delays of the multipath

therefore, $\tau_1$ is simply the delay incurred by the shortest path.

And no, $\tau_{i\ge2}$ are not relative delays; look at the formula: you just add up different independent paths with different delays; every path has its own delay, which is totally independent from $\tau_1$.


How do I implement IFFT to filter low frequencies from a group,after Freq detection by FFT?

Frequency domain processing (i.e FFT, manipulating signal, IFFT) is actually pretty complicated and requires careful choice and a thorough understanding of things like framing, step size, windowing, overlap, wrap-around, circular vs. linear etc. Simple examples for these type of algorithms are overlap-add or overlap-save for linear time invariant filtering.

If the only thing you want to do is getting rid of the low frequencies in your signal, filtering with a high pass filter is much easier and quicker.


Rather than generating a longer block of data, i would like to make a longer file by concatenating the block of 4096 samples.

Bad idea. That means your noise becomes perfectly correlated with a period of 4096 samples, and that's definitely not white noise anymore, and you'll stand a realistic chance of noticing that audibly; depending on the sampling rate, that $\frac{f_\text{sample}}{4096}$ might be a directly audible frequency.

A lot of signal processing algorithms are sensitive to periodicity, so: really, don't do that.

By the way, you might be clipping. Gaussian noise amplitude is unbound, but your sound card's acceptable numeric range very much is.

  The output file is not smooth enough. There is an audible periodic pattern that indicates the multiple 'Y' segments.

Which was predictable, because that's exactly what you've been asking it to generate.

  I believed this could be due to the imperfect combination at the beginning and End of each Y segment.

No, this was you; you introduced periodicity, and your human sound system notices that. Congratulations! Your hearing is functional :)

You using a "blend over" technique on white noise says you  haven't really thought about what white means: it means that there's no correlation between successive samples. Thus, if you concatenate two different, uncorrelated white noise sequences, you get a white noise sequence.

If you concatenate the same sequence, you get something periodic.

What you need to do is  generate white noise of the desired length instead of concatenating short sequences. You forgot to mention why you would even want to do that, because (haven't benchmarked matlab) white noise generation should really be fast.


What Distance Measure Is Proper for Image Patches?

Well, You're asking on one of the most researched topics.

There is no 1 "Right" answer.
Many models have been suggested and many of them works pretty well.

In the Non Local Means they use Euclidean Distance.
Later revisions on the idea used Mahalanobis Distance.

Others are using different norm, other use Learned Dictionaries.

Some use features or transformation of the data.

I think the best method depends on the image (Or set of images) and the usage (Denoising or something else).


determining type of filter given its pole zero plot

You'd have to figure out the frequency response of the filter. Here are two methods. I prefer Method 2 because it's quick and dirty, and you don't really care about the exact gain values in the frequency response, just the general shape to figure out the type of the filter.

Method 1: Brute Force/Computer Assisted

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

b = [1,0,1]
a = [1,0,0.25]
w,h = sp.freqz(b,a)
plt.plot( w/np.pi, abs(h) )

Clearly there's a stop band around $\pi/2$, so this is a band-stop filter.

Method 2: Hand Computation

The frequency response can be obtained by moving along the unit circle in $z$-domain. The amplitude of the frequency response can be calculated by dividing the product of the distances from the zeros by the product of the distances from the poles.

In your example the pole-zero plot looks like this:

You can trace out the amplitude of the frequency response by moving the point $P$ along the unit circle from $0$ to $\pi$ and computing $\frac{\mbox{length(PO)}}{\mbox{length(PX)}}.$

Initially,  $\mathrm{length}(PO)>\mathrm{length}(PX)$, so the amplitude response is larger than $1$. Then gradually as we move closer to $\pi/2$, $\mathrm{length}(PO)$ keeps getting shorter until it becomes exactly zero at $\pi/2$.


Is there a way to do a DFT with low latency and fine frequency resolution at low frequencies?

The answer is no.  You've explained some of the reasons why a DFT alone is not a good first step when trying to do low latency pitch estimation on low pitched human voices.

Around the first DFT bin, the effects of the window and the complex conjugate image can interfere with just spectral frequency estimation, which a pitch of that frequency may even be missing.

But it may be possible to do pitch estimation with better resolution and latency than a bare DFT by using other methods (time domain, parametric, or composite, look for many other similar Q&As about pitch estimation on this site).


Differentials - Differences: Non causality in the system

Replacing the derivative by forward differences is just a way to define a transformation from the $s$-plane to the $z$-plane, and it has nothing to do with the causality of the resulting discrete-time filter. As an exercise you could try to do the same with backward differences.

Note that this method for transforming a continuous-time filter to a discrete-time filter is useless in practice because stability is generally not preserved, and the transformed frequency response may be drastically different from the original response.


A DCT is roughly equivalent to a DFT of a vector after it is doubled by mirroring by a symmetric reflection.  This produces FFT input that does not have a discontinuity either in the middle or circularly.  A DST is roughly equivalent to a DFT after an antisymmetric mirrored extension.  This antisymmetric addition can easily result in discontinuities both in the middle and around the circle.  Discontinuities are represented by energy in the high frequency bins in the FFT results.  These high frequency artifacts are usually undesirable when using a transform for compression.

Since the DCT does not have this potential high frequency content due to circular discontinuities (as does a DST or FFT), the same total energy is thus spread lower in frequency, which potentially allows for greater compression of the high frequency DCT bins, while remaining below some visible threshold.


Difference between band pass vs. low pass + high pass?

Is there any difference to simply using a band pass filter of 0.3 – 40
  Hz?

No - but some filter design tools might not provide an option to design a band-pass filter with a unity gain "plateau" over a range of frequencies (instead, the response of the band-pass filter will be a bump at a given frequency).

  why the notch at 50 Hz when everything 40+ Hz is thrown away, anyway?

50 Hz is dangerously close to 40 Hz, and might still be in the transition band of the filter. The 40 Hz filter would have to be very steep to entirely attenuate the 50 Hz. It would be much easier to remove the 50 Hz component separately with a steep notch filter.


Scratch Detection On Noisy Textured Surface

I tried PHOT (suggested by @applesoup) using the unofficial implementation given under https://github.com/thinkng but it did not work for these image sets. Maybe one needs to further investigate or tweak the algorithm a bit.

First, if you have sufficient amount of data, I do not believe that one could easily outperform a good deep architecture in the task of texture defect classification. The literature is full of these data driven methods. Some examples are here: [1], [2], [3], [4], [5].

However, due to limited amount of data, I tried a classical approach which seems to output somewhat acceptable results. In essence, the method involves:

Compute the multi-scale difference of Gaussian (DoG) image. Call this $D=\{D_i\}$, where $D_i$ is a single DoG:
$$
D_i = I * g_i^1 - I*g_i^2
$$
where $I$ is a single channel input image and $g_i^1$ and $g_i^2$ are two Gaussian kernels with $\sigma_2>\sigma_1$. $*$ denotes convolution. As your defects are line-like, such edge enhancement should ease the job of any subsequent analysis. The resulting DoG images look like:

You could already see that the defective areas are a bit emphasized. I show here only one of the scales, whereas multi-scale convolutions form a scale pyramid (see scale-space theory). That is what $D$ is.

Perform a curve (line) detection on this DoG pyramid $D$ by using the algorithms provided in one of my previous responses. The raw response of such a subpixel line detector on a single level of the hierarchy:

Even though the images are messy and cluttered by many false positives, the lines we are looking for seem to be there. I avoid showing all the detected lines as this would severely harm the perception.

The final stage is to filter these curvilinear structures to arrive at the desired scratches. While many filtering techniques can be used, I just let the lines that are long enough survive. The output after all these stages looks like:

This is not ideal, but probably close to what you have been looking for. I also did not spend much time to tune the parameters.

Note that this approach is solely based upon a geometric modeling of the problem. For different types of scratches these assumptions won't hold. Nevertheless, I thought it might be a good starting point.


the first thing that one must understand when doing fixed-point arithmetic is that, at the basic level, it is integer arithmetic with some scaling factors applied.  perhaps just one scaling factor applied.  this scaling factor is directly related to the position of the binary point.  unless the chip is a fixed-point DSP or similar, the binary point is implied, not explicit.  it, and the associated scaling factor, exists only in the mind of the programmer.  the CPU only thinks it is working on integers.

as an example an unsigned, fixed-point, 16-bit format, that i would call "1.15" or, perhaps "Q1.15" (i think the "Q" is superfluous and notation like "Q16" does not supply enough information) has 1 bit to the left of the binary point and 15 bits to the right.  but the 16 bits, sitting all by themselves, is an integer.  the range of that unsigned integer, $n$, is $0 \le n \le 2^{16}-1$.  but the range of the unsigned fixed-point number implied $\hat{n}$ is $0 \le \hat{n} \le 2 - 2^{-15}$ = 1.99996948242188 .

the relationship between the two is:
$$ \hat{n} = n \cdot 2^{-15} $$
or
$$ n = \hat{n} \cdot 2^{15} $$ .

note that the "15"(or "-15") in the exponents is the same as the "15" in the "Q1.15" notation.  the scaling factor is directly related to the number of implied bits to the right of the binary point.

if, instead it was a 16-bit signed, we might use "s.15" or "Qs.15" notation to refer to the format, the range of the signed integer $n$ is $-2^{15} \le n \le +2^{15} - 1$.  and the range of the implied fixed-point value, $\hat{n}$, is $-1 \le \hat{n} \le 1 - 2^{-15}$ = 0.99996948242188 .

and the mapping between the two is the same as above:
$$ \hat{n} = n \cdot 2^{-15} $$
or
$$ n = \hat{n} \cdot 2^{15} $$ .

the only difference is the use of unsigned int vs. int in your C code.  they both have the same number of bits to the right of the binary point and they both have the same mapping.

now when C does addition or subtraction or multiplication or division of these values, it is operating on the integer representation of those values.  but the integer representation is the same as your implied value, except for the scaling factors.

so if you want to divide some numerator $\hat{n}$ by some denominator $\hat{d}$ to result in some quotient $\hat{q}$, you must first establish what the ranges of all three numbers are.  then establish how many bits are to the right of the binary point for each number.  then map those numbers $\hat{n}$ and $\hat{d}$ to their integer values $n$ and $d$, and determine what the range of the integer quotient is before it is adjusted to represent $\hat{q}$.  there will be a scaling factor that will pop out at you (a power of 2) which will imply shifting of either $n$ or $d$ (most likely $n$) before the integer division is done.


Yes, this is (kind of) expected:

No semiconductor device in itself is truly linear over its full damage-free range; nor is it memoryless.

You're very likely seeing spikes whenever your input approaches zero, so that your division becomes unstable (mathematically), considering you'd still be amplifying noise and an offset current.

Also, you're likely seeing sub-linear amplification when your instantaneous signal becomes large. For a proper discussion I'd just refer you to a research on "intermodulation products", which you'll find in every real (non-perfect) amplifier.

To give you a short glimpse: Let's assume your amplifier is basically a transistor with support circuitry.

A transistor's output current (assume that's emitter current) over input voltage (assume that's base-emitter voltage) curve is (simplifiedly, not considering temporal effects etc):

$$I_\text{E} = I_\text{ES} \left(e^\frac{V_\text{BE}}{V_\text{T}} - 1\right) \tag1\label 1$$

Where $V_\text{T}$ and $I_\text{ES}$ are temperature- and semiconductor dictated constants.

You'll notice this is nothing close to a linear function. Now, the job of an amplifier designer is to operate the components (mainly, transistors) in the amplifier in a way that looks as linear as possible (for power amplifiers, that is a competing requirement to efficiency and large power capability, so these aren't typically all to linear or low-noise).

Imagine you do a Taylor series expansion of $\eqref 1$. You'll notice that the exponential function has the series representation

$$e^x = \sum\limits_{n = 0}^{\infty} \frac{x^n}{n!}\text.\tag2\label2$$

Transistor amplifiers aren't really perfectly represented by that exponential function, as you can do things like influence "weights" on individual series elements within bounds by careful semiconductor design:

$$f(x)=\sum\limits_{n = 0}^{\infty} \tilde a_n \frac{x^n}{n!}=\sum\limits_{n = 0}^{\infty} a_n {x^n}\label 3$$

What you don't want is the $x^0$ term, because that's just a constant bias current ($x^0 \equiv 1$), so you'd typically try to have a very small $a_0$.

What you do want is the $x^1=x$ term, because that's your linear amplification – and typically, you'd want that to be large, so if you can, you see towards that $a_1$ becomes large.

$a_2$ and higher lead to intermodulation, and you want them to be small – but compensating them technically leads to limited $a_1$, and that means your amplifier doesn't amplify much.

Now, your amplifier manufacturer won't offer you a specific input-output curve, often, but other things more important for people that work their amplifier with signals of different frequencies:

Second and Third order intercept points. I can only recommend reading these articles closely, and possibly get a book on semiconductor practice (Tietze-Schenk would be the go-to bible where I'm from), as you're working with amplifiers and need to know what you're doing


Basics of waveform, how to define the zero point

You can't produce a DC component with speech, so the line $y=0$ is defined by the speech signal having zero mean.


First, let's start with the wrong part of your post:

  to apply fft, images must have dimensions which are the multiple of power of 2.

This is incorrect. FFTs do not have to be a power of 2.  The only requirements for applying an FFT are that you can decompose the signal length into its prime factors and that you have an FFT that can deal with these prime factors.

  to apply the fft, the dimensions of the image and the mask must be equal.

This is correct. To do the point-wise multiplication in the frequency domain, it's easier to do when both image and kernel FFTs are the same size.

  we can take the max of the width and heights of both image and the mask, round them to next power of 2 and, pad them.

The power of 2, as pointed out before, is incorrect.

the width and heights of both image and the mask, round them is also incorrect.

Your second option:

  we can take the sum of the widths and heights of the image and mask, round them to next power of 2 and, pad them.

is closer to the truth, but still incorrect.

Suppose your image is $N\times M$ and your kernel is $K \times L$. To use the FFT to do linear convolution you need pad both of them to be $(N + K - 1) \times (M + L - 1)$. If you don't you'll get time aliasing from the circular convolution.

Now, to answer your questions:

  So, in both of the cases we are actually obtaining oversized output images with padding. Am I right?

Yes, to avoid circular convolution results (and the associated "time" / "spatial" aliasing), you need to pad.

  Suppose, I pad an image and get it FFT-filtered (Sharpened/Blurred/... and so on). So, I obtain a filtered image with padding. Now, I need to remove the padding to obtain the filtered image without padding. Am I correct?

Yes. If you're interested in an in-place replacement (so that the resulting image aligns with the original unfiltered image), then you need to remove the padding.

  If YES, should I just crop the image? Or, is there any better way possible altogether from the ground up?

Yes, just crop the image at the right place.


AWGN and Rayleigh fading in slow flat fading channels in MATLAB

As I described in this answer, it is generally assumed in the Rayleigh fading model that the energy is conserved on average (i.e. that the fading coefficient $R$ is such that $E\left\{R^2\right\} = 1$). As a result you would need to generate the Rayleigh distribution from two normal distribution with standard deviation of $1/\sqrt{2}$ (instead of with unit standard deviation). This can be achieved with:

rayleigh = sqrt( (random('norm',0,1/sqrt(2)))^2 + (random('norm',0,1/sqrt(2)))^2 );

or alternatively, since you are already using random from the statistics toolbox, with:

rayleigh = random('Rayleigh',1/sqrt(2));

Other than that, the expressions you have properly represent a slow flat Rayleigh fading channel for the transmission of a one-dimensional signal (such as BPSK, M-PAM, ...)


Improving spectral response with lost downsampled data

must maintain a specific FFT length and small FFT bin width, down at baseband.

Well, the FFT bin width is defined as $\frac{f_\text{sample}}{N}$, with $N$ being the FFT length and $f_\text{sample}$ the sampling rate of the signal that undergoes the FFT. So, if both these parameters are fixed, you have no freedom in choosing the bin width whatsoever.

Note that this describes a pure view on the "tool" FFT. My guess is that you want to build some spectrum estimation – and you can certainly build a "finer" estimator.

  Is there a way to use the lost samples to improve the spectral response of the FFT and/or the SNR, given the constraints?

Since if done without aliasing, the downsampling removed no information from the signal you're interested, no.

What you can, however, do is build a better frequency estimator; that might be using the same FFT length multiple times (or, in fact, even build a larger FFT from your small FFT used multiple times), but it's hard to advise here: your resource constraint doesn't give us anything to go on – and FFTs are typically highly optimized, so if a larger FFT itself won't fit/run on your system, then chances are that complex estimators won't, either.


The odd harmonics are what makes a square wave different from a sine wave in the frequency domain.  Otherwise the FFT would lose information about the wave shape of the time domain input.

The higher the fundamental frequency, the farther apart are any multiples of that frequency, such as overtones or harmonics.


So Robert Wannamaker did a real bang-up job parameterizing the 0 dB curve and came up with, what he called F-weighting.  It's in this AES paper from the 90s:

you can see that this can be turned into an $H(s)$, and with the bilinear transform, into an $H(z)$.

here is the curve on a dB vs. log/linear freq scale.  the log frequency plot is the upside-down version of your 0 dB curve above but is more updated from Fletcher-Munson to the ISO standard which cam from i dunno where.

here's a table of poles and zeros:


How is non-integer upres downres (1080p game on a 1440p screen etc) calculated

Assuming modern 3D games here.
most games you are aware of are not fixed-resolution – there's very rarely a 1080p game these days! The GPU only renders things onto a buffer the size of your screen's resolution when asked to (in the last step of the processing chain, basically), so games can simply have a fully relative-coordinate-system going on till the end.
So, the answer is: I bet what you think is a 1080p game isn't. Other than that, the physical dimension of the screen make absolutely no difference to the image generation – it's only the last step, where, no matter what the output resolution is, an operation is applied to get pixel values out of a 3D scene definition, consisting of polygons in a 3D coordinate system that has nothing to do with pixels.
You of course need to use appropriately many polygons for high resolutions; the textures you map onto these (and thus, things like static elements, menus etc) need to be right for the rendered resolution etc.
But: there's nothing "special" about the resolution; in any case, there's a projection from one coordinate system to the screen coordinate system (pixels) in a last (or close to last) step.

You could still interpolate the 1080p pixel values based on a 1440p rendered image, but intuitively, this seems like a huge disadvantage, but when you do this in practice, the image indeed looks better than when you have 100% render scale where every displayed pixel has a 1:1 correspondence to the rendered image.

That very much depends on the game you're looking at. In effect, shrinking with an appropriate scaling algorithm is a form of antialiasing setting. That's probably what you perceive as quality advantage: smoother edges.

If I rendered 2160p and displayed this on a 1080p screen, you could average every 2x2 pixels and display it on a single pixel

This should not be done (and rarely if ever is done today); it's a bad way of downscaling things. Signal theory tells us that when filtering with a boxcar filter (in this case, 2 samples wide) in spatial domain, you get ringing, and things look jagged.
GPUs have scalers, but even as simple shader, better scaling algorithms are easy to code.


The Fourier transform of a damped cosine and the units of the result

if $v(t)$ is an electromotive force (a voltage expressed in units of volt) that is a function of time (let's say expressed in units of second), then the Fourier Transform (or the Laplace Transform):

$$ \mathscr{F} \Big \{ v(t) \Big \} \triangleq V(i \omega) = \int\limits_{-\infty}^{\infty} v(t) e^{-i \omega t} \, \mathrm{d}t  $$

will have units of volt-second.  And $\omega$ will be in radian/second.


Removing a non linear noise trend from data

The trend can be perceived as some low frequency signal, so a simple highpass filter should remove it. It will naturally slightly decrease the peaks, but by very little, because such narrow peaks are mainly composed of high frequencies.

Here is some concept code in Python (using the NumPy, SciPy and Matplotlib libraries)

import numpy as np
import scipy.signal as scisig
import matplotlib.pyplot as plt
# assumes it has just two columns with f, r values
f, data = np.loadtxt('CQDataAbs.dat', unpack = True)
# order4 IIR Butterworth highpass filter coefficients,
# the critical frequency is relative to the Nyquist frequency,
# determining the right value might require some experimentation
# note that using a higher order with such low cuttof will likely
# be numerically unstable, if you need higher orders,
# split it into several lower-order stages applied sequentially
b, a = scisig.butter(4, 0.0003, btype='highpass')
# forward-backwards filter application with edge padding (mitigates edge defects),
# results in zero phase lag and an effective filter order 8 as it is applied twice
filtered = scisig.filtfilt(b, a, data)
plt.plot(f, filtered, 'b.')
plt.show()

And here is the result

The shoot-up at the end may indicate that there are some higher frequencies in that area. Perhaps you can just disregard that region during fitting if possible.

If you were interested in the trend itself, you would obtain it by using a lowpass filter instead of a highpass filter. With zero phase lag application the result of highpass filtering should be almost identical (within numerical accuracy) to subtracting this lowpass trend from the original data.

Note that you may find that a different filter type might give you better results. A Butterworth IIR filter (with forward-backwards application to remove phase lag) is usually a good choice when the shape of the signal is important, because the frequency response within the pass-band is maximally flat. You could of course use some large FIR filter for better performance if you have enough processing power. With FIR filters you also don't have to worry about numerical instability that can occur in IIR filters of high orders.

Just in case you think about using a moving average (which is a FIR lowpass filter): That will likely perform terribly in this case, because in each window containing a peak the average will shoot up and subtracting that from the original data will greatly decrease the peaks. From a frequency perspective, the moving average filter does not mitigate higher frequencies very well, which are very significant within the region of the peaks.


Fast / Efficient Way to Decompose a Separable integer 2D Filter Coefficients without the SVD

I have taken @Phonon's answer and modified it somewhat so that it uses the GCD approach on just the top row and left column, rather than on row/column sums. This seems to handle pathological cases a little better. It can still fail if the top row or left column are all zeroes, but these cases can be checked for prior to applying this method.

function [X, Y, valid] = separate(M)    % separate 2D kernel M into X and Y vectors
  X = M(1, :);                          % init X = top row of M
  Y = M(:, 1);                          % init Y = left column of M
  nx = numel(X);                        % nx = no of columns in M
  ny = numel(Y);                        % ny = no of rows in M
  gx = X(1);                            % gx = GCD of top row
  for i = 2:nx
    gx = gcd(gx, X(i));
  end
  gy = Y(1);                            % gy = GCD of left column
  for i = 2:ny
    gy = gcd(gy, Y(i));
  end
  X = X / gx;                           % scale X by GCD of X
  Y = Y / gy;                           % scale Y by GCD of Y
  scale = M(1, 1) / (X(1) * Y(1));      % calculate scale factor
  X = X * scale;                        % apply scale factor to X
  valid = all(all((M == Y * X)));       % result valid if we get back our original M
end

Many thanks to @Phonon and @Jason R for the original ideas for this.


Here is the procedure to follow:

Threshold your image
Obtain connected components - each distinct blob in the image gives you a component.
For each component, find the boundary pixels. Each boundary is represented as a set of points.
Fit one circle per each boundary point set. You could use a simple least squares fitting.

This should give you pretty good circles. You could do this in many languages, using the standard image processing libraries.


Comfort Noise Insertion in VAD algorithms

You may want to try a two-pronged approach.  First, have a standard pre-packaged CN to inject (like jan and Jason R suggest) at the start of the conversation/exchange.

As the conversation progresses, gather low-volume data, perhaps as part of your noise reduction algorithm I assume you already have, by sampling the lowest volume signal coming in between VAD end and start points. Compare your low-volume samples to your standard CN, and replace (or convolve to simply dilute your synthetic CN) if it matches some threshold or characteristic you established. All you need is some rand and normalization functions, and perhaps FFT, all of which are readily available in C (and assembly).


I am trying to convert a signal into another unit, but the error keeps increasing during each iteration

If I understand your notation, you have $y[n] = \frac{1+x[n]}{1-x[n]}y[n-1]$, where $y[n]$ is the output signal.

This can be rewritten as $y[n] = y[n-1] + \frac{2x[n]}{1-x[n]}y[n-1]$. This is an integrator ($y[n] = y[n-1] + w[n]$) where $w[n]$ is some non-linear function of $x$ and $y$.

This integrating component is what gives you the drift in output. Note that the small positive offset in your input signal will result in a drift in your output.

If you could explain your application, and why you came up with the original difference equation, maybe some of us could come up with some insights for a solution.


I wrote code in Java based on that article but it isn't working for me. Alpha (I-data scalar) calculates to be around $100$ and $cos(\omega)$ often calculates to an imaginary number. Did anyone have any luck with this or could point me to a better IQ Correction algorithm?

(To be clear: calculateMatrix, a, c, and d are all instance variables)

public void correctIQ(int[] data) {

    int length = data.length / 2;

    // ( I'(t) & Q'(t) )
    int[] iData = new int[length];
    int[] qData = new int[length];

    // Deinterleave IQ data
    for (int i = 0; i < length; i++) {
        iData[i] = data[i * 2];
        qData[i] = data[i * 2 + 1];
    }

    // ( I''(t) & Q''(t) )
    double[] qDataOffset = new double[length];
    double[] iDataOffset = new double[length];

    // Calculate DC offset
    double qOffset = mean(qData);
    double iOffset = mean(iData);

    // Apply DC offset
    for (int i = 0; i < length; i++) {
        iDataOffset[i] = iData[i] - iOffset;
        qDataOffset[i] = qData[i] - qOffset;
    }

    // Provided equations calculated on frequency change
    if(calculateMatrix) {
        double alpha = Math.sqrt(2.0 * multiplyMean(iDataOffset, iDataOffset));
        double sinw = (2.0 / alpha) * multiplyMean(iDataOffset, qDataOffset);
        double cosw = Math.sqrt(1.0 - (sinw * sinw));

        a = 1.0f / alpha;
        c = -sinw / (alpha * cosw);
        d = 1.0f / cosw;
    }

    // Replace into interleaved array
    for(int i=0;i<length;i++) {
        data[i*2]=(int)((iData[i]-iOffset)*a);
        data[i*2+1]=(int)((iData[i]-iOffset)*c+(qData[i]-qOffset)*d);
    }
}

double mean(int[] values) {
    long sum = 0L;
    for(int i=0;i<values.length;i++) {
        sum+=values[i];
    }
    return ((double)sum)/values.length;
}

double multiplyMean(double[] values1, double[] values2) {
    long sum = 0L;
    for(int i=0;i<values1.length;i++) {
        sum+=values1[i]*values2[i];
    }
    return ((double)sum)/values1.length;
}


You may need to ignore the transients between stable pitches (for example, most consonants when singing lyrics) as these may be un-pitched sounds or noises, and thus confusing your pitch estimation algorithm.

The singer(s) may also be producing some vibrato, which should correctly show as pitch variations towards the end of some held notes.


Question about a transient response statement (Oppenheim-Schafer"s discrete-time signal processing book, 3rd ed)

It's true that Eq. $(129)$ just presents an upper bound, from which it can't directly be concluded that the transient response $y_t[n]$ decays to zero. However, stability requires that
$$\sum_{k=0}^{\infty}|h[k]|<\infty\tag{1}$$
which can only be true if $\lim_{n\to\infty}h[n]=0$, i.e., if the impulse response decays to zero. If this is the case, then the transient response
$$y_t[n]=-\sum_{k=n+1}^{\infty}h[k]e^{j\omega (n-k)}$$
must also decay to zero for $n\to\infty$.


The temporal neighborhood refers to data nearby in time. That is; data located in the frames just before and just after the current frame.

The spatial neighborhood on the other side refers to the data nearby in space. In this case this is the pixels nearby in the same frame.

Spatio-temporal filters (which is 3D filters) takes both nearby pixels in the current frame and nearby pixels in nearby frames into account. In this way you can take advantage of information from more than one frame, which is often a good idea since most video sequences does not change much from frame to frame.


The Karhunen–Loeve Transform is the equivalent of PCA analysis for continuous signals, you could seek more informations on this type of Feature extraction.

1/The idea is to compute the covariance matrix on known signals (i don't know, maybe the ECG of a person suffering from a particular heart disease).

$C = (x-\bar{x})(x-\bar{x})^T$

where X is your dataset Matrix (idk, maybe k pattern each composed of N samples), $\bar{x}$ the mean value vector (size : k,1).

2/ Decompose this matrix in eigenvector ($V$) and eigenvalue ($D$):

$V^{-1}CV = D$

3/ Extract the main features (direction of highest variablility), depending on the eigenvalues, let's imagine that the first 7 eigenvalues will represent 99% of the energy used to represent the ECG signals.

4/ Projected data $= [V^T(X-\bar{x})^T]^T$ in a more meaningful feature space

see : https://machinelearningmastery.com/calculate-principal-component-analysis-scratch-python/


Parallel equalizer, working for boosts but not cuts

Cuts can simply be derived from boosts like so: p2 = 1 / p1

This is intuitive as this is how you would invert a bandshelf.

However, I'm not sure it leads to an elegant equation as inverting the sum of filters and recovering a rational filter expression might be messy, so there's potential a better way.


How to make the $\ell_2$ norm of all columns and rows of an $n \times n$ matrix equal to $\sqrt{n}$?

HINT If we have the diagonal matrix: $$ D = \left[\begin{array}{cccc}
d_1&0&0&0\\
0&d_2&0&0\\
0&0&\ddots&0\\
0&0&0&d_n
\end{array}\right]$$

Multiplying another matrix $$M_r = \left[\begin{array}{c}
r_1\\
r_2\\
\vdots\\
r_n
\end{array}\right]$$to the left with it multiplies each row like this:

$$DM_r = \left[\begin{array}{c}
d_1r_1\\
d_2r_2\\
\vdots\\
d_nr_n
\end{array}\right]$$

If the matrix has the columns $$M_c = \left[\begin{array}{cccc}
c_1&
c_2&
\dots&
c_n
\end{array}\right]$$

Then multiplying from the right with D will give:

$$M_cD = \left[\begin{array}{cccc}d_1c_1&d_2c_2&\dots&d_nc_n
\end{array}\right]$$