Moving average vs. Moving median I have read in many places that Moving median is a bit better than Moving average for some applications, because it is less sensitive to outliers.
I wanted to test this assertion on real data, but I am unable to see this effect (green: median, red: average). See here:

Here's the sample audio data test.wav. Here's the Python code:
import numpy as np
from scipy.io.wavfile import read
import matplotlib.pyplot as plt

def median(lst): return np.median(np.array(lst))
def mean(lst): return sum(lst)/len(lst)

(fs, x) = read('test.wav')

x = abs(x)
env = np.zeros_like(x)
env2 = np.zeros_like(x)

for i in range(len(x)):
    env[i] = median(x[max(i-1000,0):i+1])
    env2[i] = mean(x[max(i-1000,0):i+1])

plt.plot(range(len(x)), env, color = 'green')
plt.plot(range(len(x)), env2, color = 'red')
plt.show()

I have tried with various values for Window width (here in the code : 1000), and it was always the same: the moving median is not better than moving average (i.e. not less sensitive to outliers).

The same with Window width = 10000 (10000 >> the spike width) :

Question:
Can you provide an example showing that moving median is less sensitive to outliers than moving average? And if possible using the sample .WAV file data-set (download link).
i.e. is it possible to do a moving median on this data such that the result is like this yellow curve? (i.e. no more spike!)

Is moving median always less sensitive to outliers? Sometimes. It will work if you have a very short spike (preferrably shorter than the median/average sample size). However, if you have a large spike, then taking the median won't help eliminate the spikes.

Long spike

I have illustrated this using a dataset (with sample size = 5, taking into account 2 left entries and 2 right entries):

Orig. Median Avg.
Data  Filter Filter
-----------------
1
2
3      3    3
4      4    4
5      5    11.6
6      6    19
40     40   26.2
40     40   33.2
40     40   33.4
40     40   27
7      9    20.8
8      9    14.8
9      9    9
10     10   10
11
12

Notice how the the median of the all the 40s is 40. For example, take the 1st 40. The left values are 5,6 and the right values are 40,40, so we get a sorted dataset of 5,6,40,40,40 (the bolded 40 becomes our median filter result).

Short spike

You also wanted an example for the median filter to work. So, we will have a short spike. Try this:

Orig. Median Average
Data  Filter Filter
-----------------
1
2
3      3    3
4      4    10.8
5      5    11.6
40     6    12.4
6      7    13.2
7      8    14
8      8    8
9
10

You can see that the median filter did filter out the single large spike, while the spike skewed the average filter's results for entries 4,5,40,and,6.


Must OFDM signal be time-division in a multiple-transmitter and SFN system when the signal is free of spread-spectrum? The OFDM baseband signal can be demodulated since  each subcarrier is orthogonal to another. But In SFN system, i.e., Single Frequency Network, all OFDM signals from several transmitters are in the same frequency domain, which breaks the orthogonality.

So without some technology like CDMA, the OFDM sigals from each transmitter must be time-divided.

Is that right?

In a multi-user scenario where every transmitter uses the same center frequency and bandwidth, but transmits different data, time-division multiplex has to be applied. Otherwise, interference between the multiple transmit signals will make successful reconstruction of a specific transmit signal at the receiver virtually impossible.

However, in typical SFNs, all transmitters are part of a broadcast network and transmit identical data. Examples are Digital Audio Broadcast (DAB) and terrestrial Digital Video Broadcast (DVB-T). In such networks, the guard interval of OFDM symbols can compensate for the propagation delay between different transmitters and there is no need for time-division multiplex.

Also, I don't think CDMA and OFDM can be combined. As a multiple-access scheme, OFDM(A) allows sharing of the same frequency band by frequency multiplex while CDMA applies code multiplex.


Defining the convolution for finite-length signals How can the linear convolution be defined rigorously for two discrete signals $x = [x(0), …, x(N_x-1)]$ and $h = [h(0), …, h(N_h-1)]$ of different finite lengths $N_x$ and $N_h$ respectively?
Let's assume for simplicity that $N_x = 3$ and $N_h = 2$, i.e. $x = [x(0), x(1), x(2)]$ and $h = [h(0), h(1)]$. Does the sum go from $0$ to $N_x - 1$ or to $N_h - 1$? However,
$$\sum_{i = 0}^{N_x-1} x(i) h(n - i) \neq \sum_{i = 0}^{N_h-1} h(i) x(n - i)$$
i.e. the commutativity cannot hold, and also both sums will have "illegal" terms containing $h(-2)$, $h(-1)$, $h(2)$, $h(3)$, $h(4)$ and $x(-1)$, $x(3)$, $x(4)$ (which I guess can be set to zero, but that somehow seems not right since that redefines my initial signals).

The OP says that sequences

$x = [x(0), …, x(N_x-1)]$ and $h = [h(0), …, h(N_h-1)]$ (are) of different finite lengths $N_x$ and $N_h$ respectively

He does not claim that $x$ and $h$ are finite-length segments of possibly longer (maybe even infinitely longer) sequences, and so the standard assumption is that finite length means what it says. In particular, $x[n]$ has value $0$ whenever $n < 0$ or whenever $n \geq N_x$, and similarly, $h[n]$ has value $0$ whenever $n < 0$ or whenever $n \geq N_h$.
The general definition of the convolution of sequences $p$ and $q$ is that result of the convolution is another sequence, which we denote as $\big(p\star q\big)$ whose $n$-th term is given by
$$\big(p\star q\big)[n] = \sum_{k=-\infty}^\infty p[k]q[n-k] = \sum_{k=-\infty}^\infty p[n-k]q[k]\tag{1}$$
subject to the usual shibboleths about convergence of the sums and the like.  So, applying this to the case at hand, we have that
\begin{align}
\big(x\star h\big)[n] &= \sum_{k=-\infty}^\infty x[k]h[n-k]\\
&= \sum_{k=0}^{N_x-1} x[k]h[n-k] & \scriptstyle{\text{because $x[k]=0$ for $k<0$ or $k\geq N_x$}} \tag{2}
\end{align}
Using the alternative formula for $\big(p\star q\big)[n]$, we have
\begin{align}
\big(x\star h\big)[n] &= \sum_{k=-\infty}^\infty x[n-k]h[k]\\
&= \sum_{k=0}^{N_h-1} x[n-k]h[k] & \scriptstyle{\text{because $h[k]=0$ for $k<0$ or $k\geq N_x$}} \tag{3}
\end{align}
Both $(2)$ and $(3)$ are valid expressions for $\big(x\star h\big)[n]$ and for each integer $n$, both give the same value for $\big(x\star h\big)[n]$ even though the expressions are seemingly quite different.. The OP claims that
$$\sum_{i = 0}^{N_x-1} x(i) h(n - i) \neq \sum_{i = 0}^{N_h-1} h(i) x(n - i)$$
but this claim is false; for every choice of integer $n$, the two sums result in the same expression (except for a trivial re-arrangement of the terms).  Don't believe it?  Let's try a few special cases keeping firmly in mind that $x[n]$ and $h[n]=0$ if $n < 0$.
For $n = 0$,
\begin{align}\require{cancel}
\sum_{k=0}^{N_x-1} x[k]h[0-k] &= x[0]h[0] + x[1]\cancelto0{h[-1]} + x[2]\cancelto0{h[-2]} + \cdots\\
&= x[0]h[0]\\
\sum_{k=0}^{N_h-1} x[0-k]h[k] &= x[0]h[0] + \cancelto0{x[-1]}h[1] + \cancelto0{x[-2]}h[2] + \cdots \\&= x[0]h[0]
\end{align}
For $n = 1$,
\begin{align}
\sum_{k=0}^{N_x-1} x[k]h[1-k] &= x[0]h[1] + x[1]h[0] + x[2]\cancelto0{h[-1]} + \cdots \\&= x[0]h[1]+x[1]h[0]\\
\sum_{k=0}^{N_h-1} x[1-k]h[k] &= x[1]h[0] + x[0]h[1] + \cancelto0{x[-1]}h[2] + \cdots \\&= x[1]h[0] + x[0]h[1]
\end{align}
More generally, assuming that $N_h < N_x$, for each $n$, $0 \leq n \leq N_h -1$,
$$\sum_{k=0}^{N_x-1} x[k]h[n-k] = x[0]h[n] + x[1]h[n-1] + \cdots + x[n]h[0]\tag{4}$$
while
$$\sum_{k=0}^{N_h-1} x[n-k]h[k] = x[n]h[n]0] + x[n-1]h[1] + \cdots + x[0]h[n]\tag{5}$$
which is the same sum as in $(4)$ except for being in reverse order! Fortunately, addition is commutative.....
That the sums in $(2)$ and $(3)$ are equal (except for being written in reverse order) is true not just for $0 \leq n \leq N_h -1$ but for all $n$. I won't bother typing out the gory details, but each sum in $(2)$ and $(3)$ contains
at most $\min(N_x,N_h)$ nonzero terms and they occur in reverse order one from the other.  $(2)$ and $(3)$ are the same sum even though the expressions look different at first glance, regardless of the OP's claims, and convolution is indeed commutative: for every integer $n$,
$$\big(p\star q\big)[n] = \big(q\star p\big)[n].$$


Proof of First Difference Property for Fourier Series I am having trouble with deriving a proof for the first difference property for the  Fourier Series.

Here is my attempt at the derivation:

$$
y[n] = x[n] - x[n-1]
$$

Fourier Series Representation:

$$
a_k - a_ke^{-jk\omega_0}
$$

Fourier Series:

$$
y[n] = \sum_{k=<N>}(a_k-a_ke^{-jk\omega_0})e^{jk\omega_0n}
$$

I have set up the summation for the Fourier Series, however I have been having difficulty to compute the summation.

The summation should equal this:

$$
a_k(1-e^{-jk\omega_0})
$$

How would I evaluate the summation?

You have the result already written down in your question. If $a_k$ are the Fourier coefficients of $x[n]$, then your third formula is the Fourier series of $y[n]=x[n]-x[n-1]$. So the Fourier series coefficients of $y[n]$ are

$$b_k=a_k-a_ke^{-jk\omega_0}=a_k(1-e^{-jk\omega_0})$$

meaning that

$$y[n]=\sum_{k=0}^{N-1}b_ke^{j\omega_0 nk},\qquad\omega_0=2\pi/N$$


Demodulating a satellite BPSK signal I am currently working on a BPSK demodulator and decoder for a satellite communication system using GnuRadio.
The system uses a BPSK modulation at roughly 2GHz with a symbol rate of 200 ksymb/s. I recorded the signal directly from the transceiver module connected to my Ettus b200mini via a coax cable and some attenuation in-between using 2.4 MS/s and therefore 12 SPS.
I create the following graph for demodulating the signal, which works while using it between different SDRs for sending and receiving. But when I try to demodulate the signal from the transceiver module, I get a really noisy constellation figure (see screenshot).
I already tried different ways of filtering, different levels of attenuation between the sender and receiver, and many different settings for the used blocks, so far nothing helped.
Can anyone help me here? Even some hints on what I could try would already be nice.
Kind regards
valkyrie

It is difficult to advise about the noisy constellation.
If it is a channel with no added noise, the constellation looks bad.
Not easy to investigate but here are some suggestions to explore:

The signal power might be too low, too close to receiver sensitivity level, so the noise you see on the constellation could be the one of the receiver.
Try to increase the power of the signal to be out of the receiver noise zone.

Automatic gain control issue: Try to freeze the power and AGC or bypass it if possible check the result

Synchronization recovery issue: Synchronize the transceiver and the receiver with their 10 MHz (so you have no carrier and sample offset). Then you can bypass symbol sync module in GNU radio

RCC filter: Maybe your filters do not have the same response in the transceiver and receiver. You could check their roll-off


Simple FFT filtering vs. e.g. butterworth filtering I am currently working on functional connectivity analysis of EEG, and need to bandpass filter my data into different frequency bands (Delta, Theta, Alpha and Beta). An important thing is that the filter doesn't cause phase distortion. I know a butterworth filter and processing the data in both the forward and reverse directions (e.g. with filtfilt in Matlab) is a good approach, however I was still wondering about the simple FFT/IFFT filtering approach. What are the exact disadvantages of this method? Does this method cause phase distortions? What are 'edge effects' that occur with this method?

Edge effects due to circular convolution can occur if you don't zero pad the data sufficiently before FFT/IFFT filtering.  Fast convolution filtering using shorter sequential segments or windows from a longer data stream requires both padding of the FFT and overlap add/save processing of the overlapping IFFT results.

Just naively zeroing FFT bins instead of using the transform of a properly designed time domain filter can lead to extreme ripples in the filter's frequency response and/or ringing in the FFT/IFFT filtered result.

No phase distortion should occur if the filter kernel has appropriate symmetry in the time domain.

The computational cost of an FFT (in CPU cycles or latency time) can be higher than with IIR filtering.  But that depends on the length of the data and the filter complexity (in filter order and/or number of passes).


Characteristic function of a random Gaussian variable I have to find the characteristic function of a random Gaussian variable of $$ \sigma_z (w)  = E e^{i w z } $$. This is the variable and I know , from the theory that the characteristic function of this variable is the Fourier transform of the probability density. From theory I read that the probability density is $$\frac{1}{\sigma \sqrt 2 \pi } e^{- \frac{(z - \mu)^2
}{2 \sigma^{2}}} $$. Now I’m blocked because I don’t know how to calculate this :/
Thank you so much

In general, the characteristic function of a random variable is related to the fourier transform of the distribution as follows:

$$\varphi_Z(-\omega) = \mathscr F \big\{ f_Z(z) \big\}$$

Why?

Because, fourier transform of a PDF $f_Z(z)$ is:

$$\mathscr F \big\{ f_Z(z) \big\} = \int^{\infty}_{-\infty} f_Z(z)\cdot e^{-j\omega z} \ \mathrm dz$$

And Characteristic function of a PDF is given by:

$$\varphi_Z(\omega) = \mathbf E_Z \big\{e^{j\omega z} \big\} = \int^{\infty}_{-\infty} f_Z(z) \cdot e^{j\omega z} \ \mathrm dz$$

As you can see the difference is just in the argument $\omega$. I think now you can get the characteristic function of a Gaussian Random Variable. Just calculate the fourier transform of the PDF and replace $\omega$ by $-\omega$.

Fourier Transform of Gaussian PDF :

$$\mathscr F \big\{ f_Z(z) \big\} = \int^{\infty}_{-\infty} \frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{1}{2\sigma^2}(z-\mu)^2}\cdot e^{-i\omega z} \ \mathrm dz$$

Now, just consider the power of exponential and simplify as follows:
$$(\frac{-(z-\mu)^2}{2\sigma^2} - i\omega z )= -\frac{1}{2\sigma^2}(z^2+\mu^2-2\mu z + i2\sigma^2 \omega z)$$ $$= -\frac{1}{2\sigma^2}(z^2 -2(\mu - i\sigma^2 \omega)z + \mu^2 + (\mu - i\sigma^2 \omega)^2 - (\mu - i\sigma^2 \omega)^2)$$ $$ = -\frac{1}{2\sigma^2}((z - (\mu - i\sigma^2 \omega))^2 - (\mu - i\sigma^2 \omega)^2 + \mu^2)$$

Now, basically, the integral becomes the following:
$$\mathscr F \big\{ f_Z(z) \big\} = \int^{\infty}_{-\infty} \frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{1}{2\sigma^2}(z-(\mu - i\sigma^2 \omega))^2}\cdot e^{\frac{1}{2\sigma^2}((\mu - i\sigma^2 \omega)^2 - \mu^2)} \ \mathrm dz$$ The second exponential can come out of the integral since it does not contain $z$, then what you are left with is just area under a Gaussian PDF with mean $(\mu - i\sigma^2 \omega)$ and variance $\sigma^2$. Area under any PDF is equal to $1$. So, basically the fourier transform of the Gaussian PDF becomes:
$$\mathscr F \big\{ f_Z(z) \big\} = e^{\frac{1}{2\sigma^2}((\mu - i\sigma^2 \omega)^2 - \mu^2)} = e^{-\frac{\sigma^2 \omega^2}{2} - i\mu \omega}$$
Once you have fourier transform, you can get the characteristic function by substituting $\omega = -\omega$ as follows: $$ \varphi_Z(\omega) = e^{-\frac{\sigma^2 \omega^2}{2} + i\mu \omega}$$


Beyond "mid-side": decompose an audio stereo recording (coincident microphones) into Mid + 45° + 90° Context: I have done in the past stereo recordings in XY position (coincident microphones):

from a source far from at least 20 meters (example: piano in a big reverberant building).
Since the microphones are coincident, there are IID (interaural intensity differences), but, sadly, few ITD (interaural timing differences). Thus the recording is less lively, much less "outside your head" (when listening with headphones), that it could have been with a spaced pair of microphones. That's the final problem that I'm trying to solve, since I can't redo the recordings: respatialize the sound into something closer to a spaced-pair recording. (I have other recordings with a spaced pair non-coicident mics, and I confirm it would have been better). See also Algorithms to re-spatialize a stereo recording audio signal?.
(Any coding idea to achieve this general goal is welcome).

The option I'm now considering to achieve this goal in this question is to:

decompose the signal L[n], R[n] (left, right) into several layers
apply different ITD (Interaural Timing Differences) on each layer
mix the layers to get a new output signal

Example with Mid-Side:
Mid[n] = (L[n] + R[n]) / 2
Side[n] = (L[n] - R[n]) / 2

Out_L[n] = Mid[n] + Side[n + K1]    # K1 is a time-shifting parameter
Out_R[n] = Mid[n] - Side[n + K2]    # K2 is a time-shifting parameter

I tried this, but it does not really help to achieve the goal mentioned above.
Question: Is there a decomposition of a stereo signal that goes beyond than Mid (0°) + Side (90°) ?
Something like:
Layer1[n] =  Mid (0° to 30°)
Layer2[n] =  Diagonal (30° to 60°)
Layer3[n] =  Side (60° to 90°)

Note: The "Coincident Microphones" part of this answer of Algorithms to re-spatialize a stereo recording audio signal?
 could be useful, but I don't see exactly how to use this concretely in an applied algorithm.

There is a generic class of algorithms called "upmixers" that indeed de-compose a stereo signal into multiple direction and/or it's constituent directional channels. They typically work by chopping the signal into frames, evaluating correlation/phase/directional differences in multiple frequency bands, steering the origial signal on rules based on the analysis and then splicing them all together again.
It's tricky business and typically involves quite bit of trade off between the effectiveness of the spatial decomposition and spectral and temporal artifacts.
However, I'm not convinced that this would help in your case. You have a far field recording done in a reverberant space: the energy of your direct sound component (which is the only thing localizable) is probably tiny as compared to the energy of the early reflections and reverb. Trying to pull out directional components out of the reverberant field is hard and chances are they are decorrelated anyway, so applying time delays to them is unlikely to make much of a difference.
If you just want a little more space in the recording, you could try to boost the difference channel a bit with a high boost shelf filter. No boost for low frequencies, some for the highs. You can play around with the cutoff frequency, but I would vary the cutoff from 500Hz to 2kHz and do try boost of 2-3 dB or thereabouts.
If you are aiming for "more natural over headphones", your best shot is to match the inter-channel correlation as a function of frequency to that of a human (or dummy) head in a diffuse field. At low frequencies the two ear signals are highly correlated, at high frequencies they tend be very uncorrelated and they are partially correlated in between. You can analyze the correlation of your recording and compare to a binaural target curve and then increase/decrease correlation at specific frequencies as needed. There are various ways of doing that, but the easiest way is to adjust the sum/difference balance: more sum = more correlation, more difference = less correlation.


Can the overlap-add method be used with varying cut-off frequency filter? I'm writing a software analogue synth. I have been using the Overlap-add method to cheaply implement a low-pass filter using FFT.

My low-pass filter has had a constant impulse response; a constant cutoff frequency. Now I'd like to make it possible to vary the cut-off frequency while a note is playing (this is a common feature of analogue synths).

I believe this is impossible while still using the circular convolution theorem trick (that let's us use FFT to implement FIR-filters).

Am I right in thinking so?
If so, is there some other trick I can use to efficiently implement a band-pass filter whose cut-off frequency varies ("sweeps")?

The problem here is not the application of the FFT for filtering but the fact that you want to model a linear time-variant system using methods for linear time-invariant (LTI) systems.

If you used the linear convolution instead of the circular convolution + overlap-add you would still face the problem that you have to divide the audio signal into blocks of data during which the impulse response is constant.

One solution is to make the data blocks sufficiently short so that the impulse response can be considered quasi-invariant during one data block. For example, if your audio sample rate is 44.100 Hz and your FFT length is 1024 every impulse response change needs to be approximated by 43 discrete steps per second.


Finding Fourier transform of a discrete signal from its Z-transform Is it possible to find the Fourier transform of a discrete signal if you know its $\mathcal{Z}$-transform of?

DTFT is the Z-transform at the unit circle. So if $z=re^{j\omega}$
then for DTFT $r = 1$.

i.e If you have the Z-transform of a signal then plug-in $e^{j\omega}$ for every $z$


What's the difference between quantization and coding? From what I understand quantization takes an analog signal and assigns numerical values to it as convert it to a discrete signal. However coding converts a signal into binary. It seems like coding is a part of quantization--if you're assigning a value to an analog signal, and you're using a computer, you'd necessarily need to use binary. Nonetheless I see diagrams with coders as separate blocks compared to quantizers. What's the difference between them?

Quantization is simply the rounding associated with going from a higher precision numbering system to a lower precision, for example analog to digital, or floating point to fixed point, or from a 16 bit system where our min to max values have $2^{16}$ levels to a 14 bit system where those same values now only have $2^{14}$ levels (so the true value will be further rounded or truncated).

Coding can refer channel encoding which is the process of introducing redundancy or memory systematically (error correcting codes) so that we can detect and correct for errors. A simple example often given is sending the same data multiple times and then voting on which answer occurred the most; this is a repetition code.

Coding can also refer to source encoding which is the process of removing redundancy to whiten the data which provides many conveniences including data compression and using the channel efficiently, but also for the ability of the receiver to work optimally for channel estimation and timing recovery.

The two processes above are often done together which is interesting in that we first source code by removing redundancy (to compress data), and then we systematically introduce redundancy (to reduce errors)!

Encryption can be considered another form of "coding".


Detecting a specific pop in a real time audio signal I am listening for a very definite pop in a real time audio signal. So far I have managed to get the audio signal in at a sampling rate of 44100.0 Hz and 2048 frames per second. I have visualized the waveform and computed and visualized the FFT. I am currently looking for distinct features to recognize this pop sound (it is the sound of a ping pong ball on the table).

I found a few papers which suggested features such as the zero crossing rate of a percussive sound like this http://www.csl.sony.fr/downloads/papers/2000/gouyon-dafx2000.pdf. However I am struggling a bit to identify the percussive envelop within the signal. I would like to identify a peak in a certain frequency range within my FFT as well, then combine the existence of an envelope with a peak in the correct range to determine whether a sound was heard.

I have two questions, how would I identify an envelope in the time domain, or how would I approach this problem. Further, how could I identify a peak in the frequency domain when I have a very fine mesh FFT (lots of bins and noise). I also have a lot of pink noise as this is a real time audio recording, so I have louder frequencies toward the bottom end than the top end, but I haven't managed to get those frequencies out of my FFT (I don't need to worry about converting back to the time domain).

Identify envelopes beside natural noise
Identify FFT peaks beside natural noise
Any experience with noticing percussive transients?

Thanks

One of the simplest methods to optimally (minimising the error energy while restricted to LTI filters) detect the existance of a known signal inside another is to use a matched filter. A matched filter will have its impulse response equal to the signal being searched in time-reversed form. This resulting filter equivalently computes the correlation of the known signal with the measured signal. And the peaks of this correlator yields the position of the best possible matches.

Let the searched signal be $s[n]$ with $s= \{-2,-1,0,0,1,2\}$ with $s[0]=-2$ and $s[5]=2$ and length of signal is $L=6$.

Assume that this signal is embedded into various locations of another longer and uncorrelated signal x[n] with its length M > N. Then the resulting combined (by adding sifted versions of the test signal s[n] into x[n] at various locations) signal is called as y[n] of length M.

Now to optimally find the index of the locations into where the signal s[n] was added, we can use an LTI filter with impulse respose h[n] which is the time reversed version of s[n] as: $$h[n] = s[-n]$$

(you can refer to Simon Haykin Communication Systems Theory for a simple and accessible derivation of this optimal filter which is called as the matched filter. Note that the impulse response h[n] of this optimal filter is the same (except time reversal) as the signal being searched, therefore it's called as "matched": optimal detector's impulse response matches the signal being searched)

Following the example of $s[n] = \{-2,-1,0,0,1,2\}$ then $h[n] = \{2,1,0,0,-1,-2\}$
Now however $h[-5]=2$ and $h[0]=-2$. This is an non-causal signal. In practice (using computer system) it is wiser to shift it right, enough to make it causal. It is just a convenience and will not effect the operation performance. Then we will have the following same impulse response $h[n]= \{2,1,0,0,-1,-2\}$ with $h[0]=2$ and $h[5]=-2$ for mathematical convenience.

Now the output of the LTI filter is $$z[n] = h[n]*y[n]$$ This is the convolution operator, which enables the computation of the output of the any discrete time LTI filter and is expressed as $$z[n] = \sum_{k=0}^{k=M-1} {h[n-k]x[k]}$$ The dummy index k ranges within the intersection of the valid signal ranges which is k=0 to k=M-1 for this particular example. Also the output index n will also range from 0 to M-1, which is by convention...

So the output sequence z[n] is expected to inlcude some peaks, these peak locations after being thresholded indicates a possible detection of signal s[n] inside y[n]. (practically the exact location will be L samples back)

What if the signal s[n] is not only shifted but also scaled? Then the peak value will be reduced as much as the scaling. Hence the threshold must be selected very carefully. Or you must find a way to estimate the scaling a priorily.


Convex Optimization with $ {L}_{1, 2} $ Regularization Term I have an optimization problem such as follow:
$$\underset{X}{\operatorname{argmin}}\sum _s \left \| T_sX_{:,s} - Y_{:,s} \right \|^2_2 +\lambda\left \| GX \right \|_{2,1} \tag{1}$$
I have introduced a new variable:
$$\underset{X,Z}{\operatorname{argmin}}\sum _s\left \| T_sX_{:,s}-Y_{:,s} \right \|^2_2 +\gamma\left \| GX -Z\right \|^2_2+\lambda\left \| Z \right \|_{2,1} \tag{2}$$
And divided this problem into two sub-problems:
$$\underset{Z}{\operatorname{argmin}} \gamma\left \| GX -Z\right \|^2_2+\lambda\left \| Z \right \|_{2,1}\tag{3.1}$$
$$\underset{X}{\operatorname{argmin}}\sum _s\left \| T_s X_{:,s} - Y_{:,s} \right \|^2_2 +\gamma\left \| GX -Z\right \|^2_2 \tag{3.2}$$
I have found that (from a paper) for the eqn 3.1, solution is as follows:
$$Z_{n,:}=[GX]_{n,:} \operatorname{}max\left \{1-\frac{\lambda}{2\gamma \left \| [GX]_{n,:} \right \|_2} ,0\right\}\tag{4}$$
I am solving this problem iteratively. First, I solve for $Z$ using equation 4, then I am solving for $X$ using equation 3.2 and utilizing nonlinear conjugate gradient method. Something is clearly wrong since I cannot reach to a solution. And it is not like solutions are diverging but more like always 0. Do you see anything wrong so far? If not, I can upload my nonlinear conjugate gradient method. Thanx in advance and have a nice day!

The problem is given by:

$$\begin{equation}
\arg \min_{X} \frac{1}{2} \sum_{k} {\left\| {T}_{k} {X}_{:, k} - {Y}_{:, k} \right\|}_{2}^{2} + \lambda {\left\| G X \right\|}_{2, 1} \\ = \arg \min_{X} \frac{1}{2} \sum_{k} {\left\| {T}_{k} {X}_{:, k} - {Y}_{:, k} \right\|}_{2}^{2} + \lambda \sum_{l} {\left\| G {X}_{:, l} \right\|}_{2}
\end{equation}$$

In the above the MATLAB notation of : is used to select a column.
It is also assumed that the Mixed Norm $ {\left\| \cdot \right\|}_{2, 1} $ is operating on each column. In case working on each rows is needed one could easily transpose $ X $.

One could see above that the problem can be solved per column of $ X $ independently. Hence it can be solved in a separable manner column of $ X $ as a vector:

$$ {X}_{:, i} = \arg \min_{x} \frac{1}{2} {\left\| {T}_{i} x - {Y}_{:, i} \right\|}_{2}^{2} + \lambda {\left\| G x \right\|}_{2} $$

Now this is a simple problem which can be solved using Sub Gradient Descent and its accelerated variants.

The OP also mentions that $ {T}_{k} $ is the DFT matrix which is a Unitary Matrix namely it preserves the $ {L}_{2} $ norm. So the problem could be written:

$$\begin{aligned}
\arg \min_{x} \frac{1}{2} {\left\| T x - y \right\|}_{2}^{2} + \lambda {\left\| G x \right\|}_{2} & = \arg \min_{x} \frac{1}{2} {\left\| {T}^{H} T x - {T}^{H} y \right\|}_{2}^{2} + \lambda {\left\| G x \right\|}_{2} \\
& = \arg \min_{x} \frac{1}{2} {\left\| x - {T}^{H} y \right\|}_{2}^{2} + \lambda {\left\| G x \right\|}_{2}
\end{aligned}$$

In this form there is no need to calculate the DFT of $ x $ each iteration.

Some of the Math derivation is given in my answer at The Sub Gradient and the Prox Operator of the of $ {L}_{2, 1} $ Norm (Mixed Norm).


Beat envelope detection I am trying to write a python function that checks whether two lets say sine signals would generate a distinct beat envelope or not. I know that the beat frequency = f1-f2 where f1 and f2 are the frequencies of my input signal. Is there a mathematical measure on when two sine curves generate a beat envelope and when they dont?

If you add two cosines you simply get the product of the sum and difference frequencies.
$$cos(x) + cos(y) = 2\cdot cos \left( \frac{x+y}{2} \right) \cdot cos \left( \frac{x-y}{2} \right)$$
If the frequencies are very close together, than the difference frequency is close to zero and that's what creates the "beat".
The exact definition of what constitutes a beat and what doesn't depends on your application. For example: in audio I would call any frequency difference of less than maybe 4 Hz or so a beat since it creates the typical beat sound.


Using I/Q demodulation/LIA for a fixed-frequency application of several signals - do the phases of the input signals affect the result? My input is a single signal composed of the sum of four periodic signals, each 4 Hz apart, centered around some resonant f0 typically between 1200 - 1300 Hz; the actual frequencies will change occasionally, but not the distance between them. For example, if f0=1250, my frequencies will be [f1..f4] = 1244, 1248, 1252, and 1256 Hz.

A decent analogy would be four LEDs flashing at these frequencies, observed by a single detector.

This signal is sampled at 16384 Hz, and I perform IQ demodulation in software (described here) using the [f1..f4] as the reference frequencies to ultimately extract the amplitude of those four components from the original signal.

My question is about the phases of the four different components. I can control their phases, and I'm wondering if there are any best practices to make sure I get the cleanest possible separation of the four signals, or if the phase won't matter?

Under the assumption that your words "extract the amplitude of those four components" actually mean "extract the spectral magnitudes of those four components", then the initial phases of the original four periodic signals do not matter. It's straightforward to prove this behavior by modeling your down-conversion process using software.


basic hsb skin detection, neon illumination i hope this is the right place to ask. otherways sorry for my mistake and pleace advice me a better site.

i'm trying to implement a super simple skin detector using some range of hsb image. i'm using approach described here and here.

i'm try to use a video source from my webcam. if i use sun illumination it works quite well (not so good but quite good), but with neon light.. it is a mess. a lot of white region are detected and a lot of noise everywhere.

why?

i'm using the algorithm described in the second source:

convert image ho HSV color space
put white on the range 0 < H < 38
dilate filter
erode filter
blur filter

This actually might work better using a simple generative model in RGB rather than HSV.

Get a training image or several training images with some skin.
Manually select the skin pixels (e. g. by creating a binary mask)
Compute the mean and covariance of the skin tone in RGB (each should be 3-element vectors)
For an unknown pixel, compute its Mahalanobis distance from the mean, using the covariance.
Classify it as skin if the distance is less than a threshold.
Tune the threshold for best performance.

Edit:
I don't know if OpenCV has a function to compute covariance, but I can tell you how to do it yourself. Let's say you have $n$ RGB pixels.  You put them into an $n$ x 3 matrix, let's call it $P$.  Then compute $m$, which is the mean RGB vector by calculating the average of the columns.  $m$ will be a 1 x 3 vector. Subtract $m$ from every row of $P$ and call the resulting matrix $Q$.  Now to compute covariance, all you have to do is multiply $Q$ by the transpose of itself: $C = Q'Q$. Make sure that $C$ is 3 x 3.

Edit2:
The values you are getting seem to be too large.  To get the maximum covariance create the following matrix:

255 255 255
 0   0   0

and compute covariance of that.  You should get a matrix where every value is approximately
32513.  So make sure that your pixel values range from 0 to 255, and make sure you copy them into floats or doubles correctly.  The Mahalanobis distance is in the units of variance, so for that the numbers should be small.  Your threshold for skin classification should probably be less than 4.


What should be the dimensions of a Fourier Transformed image? I have applied Fourier Transform to the following image.

I have downloaded this image from the Internet. I re-sized (without maintaining the aspect ratio) it using MS Paint application of Win7 to make it 512x256.

Then I have used two applications to observe its Fourier Transformed appearance.

IPLab gives the following output:

ImageJ2-20160205 gives the following output:

As you can see, the first output is a 512x256 image. The second one is a 512x512 image.

Why are those outputs different?

How would they effect processing of an image?

Here is a third different output. So I have tried to perform an FFT with Matlab, on four versions of your image:

the grayscale version of the colored image,
the red, green, and blue channels.

I do get  the following, with $512\times256$ size, in a $\log$ scale.

If as suggested by @Peter K. you zero-pad below, with $512\times512$ size, one gets:

Looks like  the FFT of the Red component, padded.
If you have the opportunity to do other tests, we can settle this riddle


Fourier transform of cosine to the power of 3 How can I find the Fourier transform of

$$ f(x) = ( \cos(x) )^3$$

I know that for $ g(x) = \cos(x) $

$$\mathcal F \Big\{ g(x) \Big\} = \mathcal F \Big\{ \cos(x) \Big\} = \pi \Big [ \delta(w-\pi / 2) + \delta(w+\pi / 2) \Big ]$$

But using this pair of Fourier transform how to obtain the $ F \Big\{ f(x) \Big\} $ ?? Is there a direct/simple way to do that?

One way would be to use the power-reduction trigonometric identity:

$$
\cos^3(x) = \frac{3 \cos(x) + \cos(3x)}{4}
$$

Due to the linearity property of the Fourier transform, you can transform each term separately and take their weighted sum to get the transform of the entire expression. The relationship we will use (from line 304 here) is:

$$
\mathcal{F}\{\cos(ax)\} = \pi\left(\delta(\omega - a) + \delta(\omega + a)\right)
$$

Which assumes that you're using the non-unitary, angular frequency definition of the Fourier transform:

$$
\mathcal{F}\{x(t)\} = X(\omega) = \int_{-\infty}^{\infty}x(t) e^{-j\omega t}dt
$$

This would yield:

$$
\begin{align}
\mathcal{F}\{\cos^3(x)\} &= \frac 34 \mathcal{F}\{\cos(x)\} + \frac 14 \mathcal{F}\{\cos(3x)\} \\
&= \frac{\pi}{4}\left(3 \delta(\omega - 1) + 3\delta(\omega + 1) +  \delta(\omega - 3) + \delta(\omega + 3) \right)
\end{align}
$$


How do I find the difference equation from a transfer function? Most of the resources I found online go the other way. If I have the transfer function
 $H(z) = 1 - cos(\theta) \cdot z^{-1} + z^{-2} $
 how do I get the difference equation from it so that I can apply the transfer function to a set of data?

A transfer function $H(Z)$ can be written as $H(Z)=\frac{Y(Z)}{X(Z)}$. Then, your $H(Z)$ can be written as

$\frac{Y(Z)}{X(Z)}=1-\cos\theta~Z^{-1}+Z^{-2}$ or

$Y(Z)=X(Z)(1-\cos\theta~ Z^{-1}+Z^{-2})$

Now taking the inverse $Z$transform, we get the difference equation as

$y(n)=x(n)-\cos\theta~x(n-1)+x(n-2)$.

I hope this help you.


Is the modulation on my WiFi different to the modulation on cellular? Is the modulation method, or nature of modulation otherwise, between my phone and my WiFi modem, versus my phone and a cellular mast, different?

From google:
3G-5G appears to use OFDM —along with QPSK og QAM depending on channel conditions.
WiFi 802.11b used some CDMA-dialect, 802.11a and g uses OFDM along with BPSK, QPSK or QAM.
WiFi mainly use frequencies around 2.4Ghz, while cellphones have used lower frequencies around 900Mhz (due to larger distances). Both have recently added capability for using higher frequencies.


Weak sense stationarity sufficient for multitaper spectrum analysis? This is a general question on multitaper spectral analysis of signals and stationarity

I work with time series which are non stationary. In this regard I have two questions

Is weak sense stationarity of second order sufficient to apply multitaper spectrum to compute valid estimates from contiguous time series objects or strick sense stationarity is required.
Multitaper works with signals that are periodic, quasi periodic, nearly periodic, and semi periodic signals however I find hard to diferenciate between periodic and the rest as to be able to see if the series I am working with belong to any of the mentioned groups.

Your help will be welcomed

The answers to the questions are

Yes WSS is sufficient.

Quasiperiodic

Semiperiodic


Even more on: Kalman filter for position and velocity As a follow-on to this great question...

I'm trying to implement the model in the referenced article (and referenced Wikipedia article), but I'm not getting the results I think I should be getting.  Ultimate I would like to take this example and marry it with intermittent GPS data so as to simulate accelerometer failover in situations where GPS is lost.

I've build the following model in Simulink:

UPDATED MODEL:

and power it using the following UPDATED Matlab code:

% Configure simulation
sim_length=10;
Ts=0.1;

% State space equations
A = [1 Ts;0 1];
C = [1 0];

% Noise characteristics
sigma_a = .1;
process_noise_power = 10;
measurement_noise_power = 10;

% Noise matrices
G = [(Ts^2)/2; Ts];
H = [1];

% Process noise
Q = sigma_a^2;
% Q = G*G'*sigma_a^2;

% Measurement noise
R = measurement_noise_power;

% Steps (not needed)
N = 0;

% Run simulation
sim('full_discrete_2inputs');

% Show results
close all;
figure(1);
hold on;
plot(tout,x_actual,'b-');
plot(tout(1:length(y_measure)),y_measure,'m-');
plot(tout(1:length(x_est)),x_est(:,1),'r-');
legend('True Position','Measured Position','Filtered Position','Location','Best');
title('Actual vs Kalman Filtered Position');
xlabel('Sample time [s]');
ylabel('Distance from origin [m]');
grid on;
hold off;
shg;

And here are my UPDATED results:

So with that, I have the following questions:

Did I inject the noise at the correct location? Is it appropriate to inject it prior to the double integration?
What's a good real-world noise power level to use in this application?  I picked a number at random, but I imagine whatever I pick has to be considered in my Q or R matrices.  I know if I bump the noise up to something like 50, the filter output and measurement value go wild, but still track each other.
To get the simulation to work I had to change my H matrix to 1, but [1 0] was specified in the example.  Does anyone know why?

There are some problems:

  1-Did I inject the noise at the correct location? Is it appropriate to
  inject it prior to the double integration?

The noise you injected accounts for process noise at the acceleration state, which you seem to assume as constant plus random disturbance. But then your true state generation is wrong.

  2-What's a good real-world noise power level to use in this
  application? I picked a number at random, but I imagine whatever I
  pick has to be considered in my Q or R matrices. I know if I bump the
  noise up to something like 50, the filter output and measurement value
  go wild, but still track each other.

Which noise? Measurement or Process noise? Where is your measurement noise added? On the simulink model or on the script file?

  3-To get the simulation to work I had to change my H matrix to 1, but
  [1 0] was specified in the example. Does anyone know why?

H =[1 0] is correct if your state $X = [x ~~\dot{x}]^T$ has 2 elements such as position $x$ and velocity $\dot{x}$. Becasue your measuring just the first state $x$, the position.


Interpolation with an FIR filter How can I use interpolation with an FIR filter? I am more familiar with interpolation such as nearest distance interpolation, linear interpolation and so on.

Suppose a signal is given as the following series:

    time=    1 2 3 4 5 6 7 8
    signal=[ 1 1 2 2 3 2 7 6]

Then, if I ask what's the signal value at time 2, we can see signal[2] = 1; then how about the signal value at time 2.1? With a linear interpolation method, I can tell signal[2.1] = 0.9*1+0.1*2 = 1.1.

In this case, how could I know the corresponding FIR filter? How could I use the filter to obtain the signal value?

Let's view this from a slightly different angle:

If the signal is properly sampled, i.e. in accordance with the Shannon/Nyquist criterion, then the samples contain all information about the original signal. If not, all bets are off, so we'll skip this for now.
Interpolation is then equivalent to sampling the signal at a non-integer time. In your case you want to know what is x(t=2.1) by using the information x[1], x[2], x[3] ...
There is only one "correct" answer to "what is x(t=2.1). In your example it happens to be 1.05061339.... It can be calculated using the Whittaker Shannon interpolation (see Wikipedia) as $$x(t)=\sum_{n=-\infty }^{\infty }x[n]\cdot \mbox{sinc}(\frac{t-nT}{T})$$
This can be interpreted as the convolution with an impulse response. However, it's not "finite".
Unfortunately The Whitaker Shannon interpolation is not very practical: The impulse response is infinite, the number of zeros is infinite and it's falloff with time is slow, so it's not a great candidate for windowing to get a time limited impulse response.
The Whittaker Shannon interpolation is equivalent to convolution with the impulse response of an ideal low pass filter. All practical interpolation methods will also involve  a low pass filter. Even linear interpolation can be interpreted as a low pass filter, it's just a very bad one.
Interpolation filters are made "practical" by adjusting filter parameters to the requirements of the application: how much precision do you need, what's the spectral content of the signal, what latency can you tolerate, etc.
Filter parameters to adjust are filter order, cutoff frequency, rolloff, filter design method (equiripple, least squares, etc), linear or minimum phase, side lobe height, stop band attenuation, pass band ripple, etc.


How to decode speaker gender from spoken transcript I have 20 audio recordings (raw format) of people reciting the same ~50 word transcript. I would like to use audio signal analysis to decode speaker gender from audio recording.

Intuitively, it seems I need to align the audio signals and then cluster them. Assuming this is correct,

For alignment, would you recommend cross-correlation?
For clustering, would you recommend principal component analysis vs independent coordination analysis?

Thanks!

From your question, I take it that you are attempting to analyze, in a completely unsupervised way, unlabelled speech data so as to unveil two groups that should (hopefully) line up well with the "male" and "female" labels.

Some notes:

Cross-correlation is useless for speech application. The notion of speech similarity (be it at the higher level of the actual words pronounced, or at the lower level of the perceived timbre) does not line up at all with the raw waveform correlation. If you really want to align pairs of recordings, you'll rather have to match sequences of feature vectors (such as MFCC) rather than waveforms. Even then, the MFCC of two different speakers saying the same words do not match - at the very least you have to consider a linear transform from one set of vectors to the other.
You have not explained how you planned to use principal component analysis here (PCA is not a clustering method in itself). I assume that your intention was to represent each speaker with a humongous $N$-dimension vector representing the entire waveform of the recording for this speaker; run a PCA on this $20 \times N$ matrix to get a lower dimensional representation of say $20 \times 2$, and use this for clustering. It won't work because without further feature extraction, the raw waveform does not carry any meaningful information for speech analysis.

Here's how I would naively do it (and this wouldn't require any alignment):

For each recording, extract the pitch and MFCCs (and maybe their first and second derivatives) and train a GMM on them.
Compute the pairwise distance between the models. This gives you a $20 \times 20$ matrix with the distance between each model.
Once you have this matrix, it is ripe for use with spectral clustering techniques, or maybe something as simple as hierarchical agglomerative clustering.

Note that in this answer I have tried to follow your "clustering" approach - how to make difference in genders apparent in a collection of speech recordings, without any prior annotation.

If your problem is truly that of recognizing gender, and if you can afford a supervised approach (annotating data + training a model), then things would have to be done differently...


Fourier transform exercise I have this signal:
$$
X(f)= 2\delta(f)+ \delta\left(f-\frac 1{T_0}\right)+\delta\left(f+\frac 1{T_0}\right)+\textrm{rect}\left(\frac{f-\frac 4{T_0}}{\frac 2{T_0}}\right)+\textrm{rect}\left(\frac{f+\frac 4{T_0}}{\frac 2{T_0}}\right)$$

I must calculate the energy. How can I find
$$\int \lvert \delta(f-f_0)\rvert^2df\quad?$$

The signal, whose total energy you want to calculate, is periodic therefore it will have infinite energy...

To see that note, the following Fourier transform pair:
$$
x(t) = \cos(2\pi f_0 t) \longleftrightarrow X(f) = 0.5 \delta(f + f_0) + 0.5 \delta(f - f_0)
$$

And based on Parseval's relation,
$$
E_x = \int_{-\infty}^{\infty} |x(t)|^2 dt = \int_{-\infty}^{\infty} |X(f)|^2 df
$$

you can conlude that total energy is infinite


How should one select this filter coefficient for a pre-emphasis filter? My professor asked me why you choose $a=0.97$ in a pre-emphasis filter $H(z)=1-\frac{a}{z}$ for your pre-processing step.

What is the relation between $a$ and the filter?
How does this coefficient effect the performance?
Why do people use $a=0.9$ or $a=0.97$?

If you look at the spectrum of speech sounds, you will notice that less energy is present in the highest frequencies, with an overall decreasing slope. The goal of the pre-emphasis filter is to counter-balance this, and flatten the spectrum. This helps for one or several of the following reasons, depending on your application:

In an analog speech transmission/recording system, applying an emphasis filter on the emitter (recorder) and a de-emphasis filter on the receiver (player) improves the signal to noise ratio for the high frequencies - since they will be recorded/transmitted at a higher amplitude than without pre-emphasis. This is why vinyl records need a special kind of pre-amplifier (RIAA equalization).
In a digital speech processing system implemented with fixed point arithmetics, this avoids truncation errors, by allowing the FFT coefficients carrying the highest frequencies to have a larger magnitude, making better use of the available dynamic range.
In a speech analysis/modeling context (such as linear prediction or sinusoidal modeling), one seeks the model parameters which minimize the distance between the model and the signal. When applying these techniques on speech signals which have not been pre-emphasized, the model will spend a lot of its "parameter space budget" (sinusoids, poles...) trying to overfit the very first harmonics, while ignoring the highest order ones - while these harmonics might actually be relevant for speech intelligibility or phoneme discrimination. Applying pre-emphasis ensures that the model will fit the spectrum more evenly.
The low frequency band is occupied by sounds which are useless/harmful for speech recognition: DC offsets, microphone pops, AC 50Hz/60Hz buzz from ground loops or poor board design, etc. The pre-emphasis filter you mention, which is indeed a simple high-pass filter, helps getting rid of these.

As for the the actual value of $a$, you can try plotting the frequency response of the filter for these different values and/or see where the cutoff frequency is. You could also try to visualize its effect on a small corpus of speech signals - how balanced the spectrum becomes after pre-emphasis. It is largely application-dependent - for example, in some applications, poor microphone quality adds an extra 6dB roll-off of the higher frequencies which need more compensation. In speech recognition, this is the kind of parameter tuned a posteriori to try to get better recognition results - though my experience is that with modern models and recognition techniques it has little influence on performance.


Filter H(z) manipulation frequency response changes I'm trying to understand how altering the frequency response of a H(z) low pass filter, will visually alter it's frequency response plot.
For example: by doing H(z^2), would the frequency response represent the filter's cut out frequency happening much later than H(z)?
If I multiply H(z) for something like 0.5(1-z^(-1)) as I have zeros in 0 and 0.5 a sort of band-pass will show up between 0 and 0.5?
Thank you

It is just like ordinary functions manipulation.
You have to remember that z and the frequency are linked with $z = exp(j \omega \cdot T)$
$H(z^2)$ you are doing something equivalent to $H(2 \omega)$ on the continuous frequency domain, i.e. you are narrowing the filter band, and thus stretching its time response.
On the second example, you have to remember that $H_1(z) H_2(z)$ is the response of two cascading filter, that the product of the Z transforms is the convolution of the the time responses.
The filter you are suggesting $0.5 (1 - z^{-1})$ will be a sort of derivative estimate because a signal with response $x[n]$ will be transformed to $(x[n] - x[n-1])/2$ and this is a high-pass filter. You can also think of it as the inverse of the step response, you can spend some time contemplating a table of pairs of signals and their Z transform, for instance this and slowly but certainly you will build some intuition.


Low-Complexity Audio Compression? What options exist for discrete audio compression that is good quality, moderate compression (e.g. 3+:1), and low complexity (e.g. hundreds of processing cycles to encode each sample without the use of floating point), especially ones that are not patent-encumbered? The specific application is ambient sound - vocal and nonvocal audio.

I'm not very well-versed in this space of algorithms. I'm familiar with the high quality high compression algorithms, but even MP3 is too expensive for my purposes. The only low complexity moderate compression algorithm family I know is ADPCM, and I have no idea if any implementations of that could be considered good quality.

EDIT: Relevant question: what is the highest quality fixed-point ADPCM algorithm?

SBC was specifically created for this purpose. See http://en.wikipedia.org/wiki/SBC_(codec). I don't know what the current legal situation is.

Other alternatives depend on your application: what is your target bandwidth and what type of artifacts and errors can you tolerate and which ones not.


Improving Template Matching Algorithm for OpenCV Car detection using OpenCV
I am working on a small personal project where i have to know wether an image shown is a car or not. I am using the basic OpenCV python based template matching.
Positive Image / Template Image
Car Top View  :-

The simple template matching by using one of the positive image on the other is giving the required result.
Negative Image
But when we are using negative images like

the template finder is finding this as a positive match.
Methods tried and failed :-

Increase the threshold for template matching :- Increasing the threshold is causing many of the actual positive images to stop matching

Different types of template matching techniques that are there in OpenCV :- Not giving any better result

Using multiple templates to reinforce positive and negative matches:- Doesn't work well for all the cases

Is there  a better way to remove the template matching false positives that we are getting. Are there any feature matching or edge detection based techniques that can be used instead of template matching to improve my algorithm.

I assume that the image have same scales in terms of sizes.
With that out of the way, I would suggest try using the following:

Break up your image into parts. (the template image as well as the test image. may be 3 or 4 according to your choice.)
Then perform template matching on each part. Consider common threshold for every template part and then count it as positive only when all parts give a positive result.

(P.S. A similar solution was suggested on some other post. I don't exactly remember which post.)

Update: Partition your template and test images at exactly same points.


Difference between OFDM and FDM In Frequency Division Multiplexing we provide one channel to one user.

Example like BandWidth of user1's voice is 200Khz. I provide one channel of 200KHz in FDM for uplink as shown in the figure.
Now I will provide another 200KHz Bandwidth(channel 2) for user2's voice. Similarly for user3.
I will be having a guard band of 100 Hz for this.

Now What I understand about OFDM is that channel1 will be orthogonal to all other channels. Also, we divide channel 1 into orthogonal sub-carriers with bandwidth less than 200KHz.
Am I correct? If yes then will there be no ISI as my required Bandwidth is 200KHz.

Both are orthogonal in your example, in different ways. The FDM example provides orthogonality by properly spacing the signals in the frequency domain (their dot product is 0). OFDM is far more spectrally efficient where you have overlapping carriers that are orthogonal.

From an implementation standpoint OFDM can be achieved with a DFT. Each of the carriers is a much lower bit rate that turns a frequency selective fading channel into flat fading. This is not true for your FDM example. There are many other differences but these are a few.


Kronecker delta instead of dirac delta as correlation function of white noise From my understanding, if you have a sample $x_{t_1},\dots,x_{t_n}$ of $X_{t_1},\dots,X_{t_n}$ which are iid $N(0,1)$, then $x_{t_1},\dots,x_{t_n}$ is a sample path of Gaussian white noise.

However, it is stated that the correlation function $c(s,t)$ of white noise is $\delta(s-t)$ where $\delta$ is the Dirac delta function. I understand that if $s\neq t$ then $c(s,t) = 0$. However, if $s=t$, you want $c(s-t)=1$. Why isn't the kronecker delta function used instead?

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


Comparing Image FFT with Sine wave FFT So, When we take the FFT of any Composite or Single Sine Wave we get exact frequency on x-axis of the plot.
The point where we get the original frequency in FFT is depend upon the Sampling Rate and the Length of the actual signal, i.e

If there is Sine wave of 100KHz and the Sampling Frequency equals to 1MHz, and If we are plotting the 10 Oscillations i.e t = 10/f then the length of the Signal would be 100.
and If we plot the FFT of that signal than the First point in FFT represent 10KHz i.e 1MHz/100 = 10KHz. So the actual signal frequency will be plotted on the 10th point of FFT graph, ( 10KHz*10 = 100KHz ).

This is how we can see the Frequency of our actual signal in Frequency domain.

but what about an Image ?

Consider a binary 8x8 image ,

image = [zeros(8,4) ones(8,4)]

image =

     0     0     0     0     1     1     1     1
     0     0     0     0     1     1     1     1
     0     0     0     0     1     1     1     1
     0     0     0     0     1     1     1     1
     0     0     0     0     1     1     1     1
     0     0     0     0     1     1     1     1
     0     0     0     0     1     1     1     1
     0     0     0     0     1     1     1     1

The FFT will be calculated for each row and for each column which mean the number of samples to calculate each FFT will be 8.

I know that the FFT values have no frequencies, it has only complex amplitudes, then how does the FFT is plotted for image, How can identify the actual frequency like i did in case of Sine wave by dividing the Fs with length of Signal ?

Here is 1D FFT of the Signal,

x = [0 0 0 0 1 1 1 1]; %what are the frequencies here

    % Samples = 8

    f = fft(x)
    plot (abs(f))

Can somebody explain me please what are the frequencies in the FFT plot ! ???
I just need one good example that explain the graph of FFt of this signal like i did for a sine wave.

x = [0 0 0 0 1 1 1 1] is one cycle of a (not bandlimited) square wave, with a DC offset, with a period equal to the period of the window.

....''''....''''....''''

The DC component is the FFT value at index 1 (remember that Matlab's array indexing starts at 1 instead of the more sane 0-based indexing), and to get the true value you have to divide by the number of samples (4/8 = 1/2 DC offset).  The square wave is the components at index 2 and 4.  Because the signal is real, the components at 6, 7, 8 are just mirror images of 4, 3, 2.

A square wave is made up of odd harmonics with 1/n amplitude, so with fundamental f1

frequency   f1, amplitude   1
frequency 3*f1, amplitude 1/3
frequency 5*f1, amplitude 1/5
frequency 7*f1, amplitude 1/7
...

Since your number of samples is so low, you only have room for the first 2 harmonics before you hit the Nyquist frequency.  So the value at index 4 is the 3rd harmonic, which should be approximately 1/3 of the value at index 2, the first harmonic (though not exactly 1/3 because of aliasing).

If your sampling frequency is 1 MHz, and the value at index 2 represents the fundamental frequency, then you know the square wave fundamental frequency is (index-1)*fs/N = (2-1)*1MHz/8 = 125 kHz, or 1/8th of the sampling frequency, as you would expect.

Try experimenting with longer signals and it will be easier to see what's happening.


Where can I get a database of speech signals? For a classification sort of experiment using Matlab, I am in need of a database of male and female voices. So where can I get links of such open databases.

you can also find speech corpus from here

English:

American National Corpus Open to All

Huge Corpus from Vox Forge Repository

Indian Local:

Hindi Language Corpus

Bengali Corpus

Speeches from Indian Language Consortium


Determining optimal binary decision rule threshold from observations with unknown priors? Given only observations of a binary signal perturbed by Gaussian noise with unknown prior information, how can I estimate the optimal decision threshold?

( No, this is not a homework question)

Specifically, I am think about the following model: $Y$ is a two-state $(H_0,H_1)$ Random variable :

$P(Y|H_0) \sim \mathcal N(\mu_0,\sigma)$
$P(Y|H_1) \sim \mathcal N(\mu_1,\sigma),\quad \mu_0 < \mu_1$
$P(H_0) = \pi_0$
$P(H_1) = 1-\pi_0$

with unknown parameters: $\mu_0, \mu_1, \sigma, \pi_0$.

The Maximum a Posteriori Log-likelihood threshold could be computed from those parameters if I knew them. I was originally thinking about how to estimate the parameters first in order to get to the threshold $Y_t$. But I'm thinking it may be more robust to directly estimate $Y_t$.

Thoughts:
Normalizing the observations (subtracting the sample mean and dividing by standard deviation) reduces the parameter space into 2 dimensions: $\pi_0$ and $\frac \sigma{\mu_1-\mu_0}$.

My intuition is that it would be difficult to derive the right decision threshold you expect to find:

$$\tau = \frac{1}{2}\left(\mu_0 + \mu_1\right) - \frac{\sigma^2}{\lVert\mu_0 - \mu_1\rVert^2} \log \frac{\pi}{1 - \pi}\left(\mu_0 - \mu_1\right)$$

From the global statistics you are considering (sample mean: $\pi \mu_0 + (1 - \pi) \mu_1$ ; standard deviation: more complex expression but I doubt it would involve a log).

I would approach the problem this way:

If the assumption that $\sigma$ is small can be made

I'm mentioning that, because keep in mind that the decision threshold is affected by $\pi$ only if $\sigma$ is sufficiently high to allow both classes to overlap. If the $\mu$s are distant by more than a few $\sigma$, class prior probabilities have nothing to say in the decision process!

Run k-means on your observations ($\sigma$ is small and is shared by both classes, so k-means is in this case EM for the mixture model). If you just want to binarize these observations and no other data, you can stop here.
If you have new observations to binarize, and you know they are generated by the same process, you can use the class centroids found by k-means on your training data as estimates of $\mu$, and use the middle as a decision threshold.

If no assumption about $\sigma$ can be made

Run the EM algorithm (with pooled, diagonal covariance) on your training data. Use the inferred "soft class membership" variables to binarize your observations.
Compute the decision threshold $\tau$ from the parameters given by EM to binarize new data generated by the same process.


What are the practical constraints on designing Sensing matrix in compressed Sensing? In a typical compressed sensing scenario,  $y=Ax$, where $x$ is a sparse signal and $A$ is the sensing matrix.
To reconstruct the sparse signal $x$ from $y$,  $A$ should posses the Restricted Isometry Property (RIP). Therefore, if a matrix does not possess RIP the reconstruction of $x$ will not be possible. How can we check a matrix for RIP ?
Secondly what are the practical constraints on designing the sensing matrix $A$ ? Like should the matrix consists of real values or complex? What distribution to use for generating the matrix $A$ ? etc,.

Checking for RIP of a matrix is an NP-Hard problem which means it is not computationally feasible to accomplish. RIP is used in matrix design mostly in theoretical aspects. Stealing @David 's comments, RIP is a tight condition.
To design sensing matrices, another weaker condition named Coherence is used, which its computational complexity is tractable ($O(n^2)$ where $n$ is length of matrix). No matter what type of elements form the matrix, you can always check for its coherence.
A similar question has been asked here:
Compressive sensing: numerical generation of RIP matrices . Check it out for more details.


Phase of a signal with multi delayed windowed (rectangular window) signal with different amplitudes I have a signal $x(t)$ which is composed of multi delayed linear chirps with different amplitudes and phases but they share the same $w$ and $\alpha$, the signal can be expressed as follows:

$$\DeclareMathOperator{\rect}{\mathrm{rect}}
 x(t) =\sum_{i} a_i \cos\left(\omega t+\frac{\alpha}{2} t^2 +\phi_i\right) \times \rect\left(\frac{t-T_i}{\Delta T_i}\right)
$$

where $a_i$ can be any real number. The signal $x(t)$

differentiated w.r.t   time $t$

My goal is not to estimate the different parameters $a_i$,$\phi_i$,  $T_i$ or $\Delta T_i$ but to estimate the time varying phase term $\omega t+\frac{\alpha}{2} t^2 $.

Any ideas/tips?

Assuming they all have the same $\omega$ and $\alpha$, meaning they are all following the same frequency ramp, then you can multiply by a ramp in the opposite direction of a starting guess while measuring the derivative of the phase of each carrier (a simple approach to see this would be to observe the result with an Short-Time Fourier Transform). Point is once you have found (and removed) the frequency slope, you will be left with distinct and constant frequencies, offset in frequency based on the time difference of their starting position and the frequency ramp rate (which also gives you information as you converge toward the final frequency ramp).

Since you do not care about the amplitude, then a metric that simply looks at the frequency components independent of amplitude would work well for this.

For the frequency ramp correction, it would be simplest to do that by multiplying by a complex frequency ramp ($e^{-j(\omega t + \alpha /2 t^2)} $ ) rather than a cosine and then observe just the positive frequencies in the complex output (What you would observe from bins 0 to N/2 in the STFT).


How to remove the bump from the signal using filters? I want to have only the slope of the signal. How can I do? I have tried to apply some filters but so far I have not been able to hit the nail on the head. I use Matlab.

A bump like this one is likely to be wide-band, especially with the sharp onset. Plus, the line may be hard to deal with in the Fourier domain.
 Hence, the combination is complicated to remove with a classical linear filter. The problem is very akin to baseline, background or trend removal, answered elsewhere here.

Several options are possible, for instance:

use a non-linear filter, based on a median, or a minimum/maximum statistics,
use morphological operators: a rolling ball, lot of straight segments, etc.
use a knowledge on the data model, like a linear equation: $y=ax+b$, or the fact that the bump is "above",
combine the above in a variational formulation, using appropriate data fidelity and penalty.

In your example, I suspect that a classical linear fit with robust distance (like a least-absolute distortion) could do the job. I will call all the above filters, in the wide sense that you will remplace a value with respect to some sort of combination of the others.

You can also call the following robust regression, LAD fitting. An example at work:

% Standard and Robust fit of a degree 1 polynomial w/ a bump
nSample = 1000;
% Create a similar composite signal
time = linspace(0,5,nSample)';
polyCoef = [0.2  0];
dataLine = polyval(polyCoef,time);
dataParabola = -8*(time-2).*(time-3);
dataParabola(dataParabola < 0) = 0;
data = dataLine+dataParabola;

% Use Matlab curve fitting toolbox
optsRobust = fitoptions('Method','LinearLeastSquares','Robust','LAR');
[fitObject,gof] = fit(time,data,'poly1',optsRobust);
h1=plot(fitObject,time,data);
grid on


How to estimate radial distortion from lens characteristics? Is there a way to estimate the radial distortion coefficients without calibration, given the lens characteristics, such as its diameter and its focal length?

Typically for camera calibration a pinhole camera model is used, and then radial distortion is added on top of that. Radial distortion is a 4th or 6th degree polynomial with 2 or 3 coefficients. During calibration, these coefficients are not derived from anything. They are initially assumed to be 0, and then estimated from data using numeric optimization.

What I am asking is whether or not there is an equation relating lens characteristics that you can read off the manufacturer's spec (focal length, aperture, field of view) to the radial distortion coefficients from the typical polynomial distortion model? In other words, is it possible to get some approximation of the distortion coefficients from the camera spec, before doing the calibration.

There is a lot of variability because distortion can be compensated for optically and lens designs differ. Some lenses are marketed as rectilinear, most not. You would be making the distortion worse for some lenses by doing any correction advised by focal length alone. For focal lengths shorter than about 30 mm, the distortion seems dominantly barrel, according to the data of:

Neale, W., Hessel, D., and Terpstra, T., "Photogrammetric Measurement Error Associated with Lens Distortion," SAE Technical Paper 2011-01-0286, 2011, doi:10.4271/2011-01-0286.

They only tested zoom lengths though.


why LTE doesn't use DC subcarrier in the downlink? I know LTE doesn't use the DC subcarrier for OFDM/OFDMA modulation for the downlink, but no author gives a detailed explanation to this stuff.

Can you tell me the reason behind it?

Thank you!

The DC subcarrier is normally omitted because in a zero-IF receiver it overlaps with DC offsets in the hardware which are very difficult to remove.


What is the coherence bandwidth of the channel? I read the following sentence.  "In wideband systems, the transmission bandwidth of a single channel is much larger than the coherence bandwidth of the channel."  What is  'coherence bandwidth of the channel'?

If you have a multipath channel with impulse response

$$h(t)=\sum_{k=0}^{m-1}a_k\delta(t-\tau_k)\tag{1}$$

the corresponding frequency response is

$$H(j\omega)=\sum_{k=0}^{m-1}a_ke^{-j\omega\tau_k}=e^{-j\omega\tau_0}\sum_{k=0}^{m-1}a_ke^{-j\omega(\tau_k-\tau_0)}\tag{2}$$

Assume that the delays are ordered such that $\tau_0<\tau_1<\ldots <\tau_{m-1}$. Note from $(2)$ that the periods of the complex exponentials are given by $1/(\tau_k-\tau_0)$, and the smallest period is $1/(\tau_{m-1}-\tau_0)$, which is simply the inverse of the largest delay difference. This largest delay difference is called the channel's delay spread. Its inverse, i.e., the smallest period of the complex exponentials in the frequency domain is referred to as the coherence bandwidth.

The importance of the coherence bandwidth is the following: in a frequency range which is considerably smaller than the coherence bandwidth, the channel's frequency response can be considered flat (i.e., it doesn't change much). So if your signal's bandwidth is small compared to the channel's coherence bandwidth, it will not be distorted very much. If the signal bandwidth is larger than the coherence bandwidth, the signal will experience frequency-selective fading.


What does it mean to combine signals coherently? What does the expression 'combine signals coherently' imply?

It means that they are in phase---think of the signals as waves. Adding noncoherent signals results in cancellation, or fading.


How do I resample an audio signal without altering its pitch? I'm working on concatenating videos again, this time with different playback speeds. The video stream part is easy enough, but the audio part is what's stumping me at the moment.

Currently, I've managed to do it by just multiplying the timestamp by a factor inversely derived from the desired playback speed, so that if the playback is to be half speed, the timestamps are multiplied by two.

Not only that, I've also altered the audio samples (in the form of a byte array) so that when the playback is half speed, every sample (of four bytes) is repeated, and when the playback is double speed, every second sample is discarded.

My problem is that this method also alters the pitch of the resulting audio track - the half-speed playback lowers the pitch, the double-speed playback raises it.

How can I resample these samples without altering their pitch?

Simple resampling would not work in your case. Resampling, while stretches/compresses the signal in the time domain, also does an inverse operation in the frequency domain.

Some basic time stretching algorithms you might want to look into (which do not result in pitch change) are:

1) OLA - overlap and add

2) PSOLA - pitch synchronous overlap and add

These work well for monophonic signals if that's what you have.

3) Phase vocoder

This works for polyphonic signals as well. But you may have to be careful about transients while using a phase vocoder approach.

There may be time stretcher plugins you could look at which use these. You could also find some open source code and time shift the signals in MATLAB/Python, etc.

Hope this helps. Let us know how it works out!


Remove audio signal component under certain frequency in python I want to remove all signal components under 10KHz in one audio file. Below is my code
from scipy.signal import butter, filtfilt

def butter_highpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='high', analog=False)
    return b, a

def butter_highpass_filter(data, cutoff, fs, order=5):
    b, a = butter_highpass(cutoff, fs, order=order)
    y = filtfilt(b, a, data)
    return y

(fs,rate)=wav.read('46.wav')
y=butter_highpass_filter(rate,10000, fs, order=5)

Below are the two outputs before and after applying the HPF.

Why should the be a mistake?
A linear system (like your filter) can never introduce frequency components that weren't there. So, the tone can't move anywhere. And because you don't show the mapping of colors to power, you can't notice any change in amplitude.
So, I bet your signal is properly attenuated, you just forgot to actually check whether it is.
Also, filfilt is almost certainly not the function you're looking for, read the documentation to that function before using it. You're even importing lfilter!
Furthermore, you label this as "BB (baseband?) signal" in your plot. Further guess: you really don't want a Butterworth filter, but a FIR with a linear phase.
(Using an FIR with a linear phase would most definitely make filtfilt obsolete, as you can just shift the signal by the constant group delay.)


Initialize Hi Pass Filter I want to remove the DC from a ADC, so I designed a two-pole, 1 Hz high-pass Butterworth filter for a signal sampled at 10kHz. The filter coefficients are:

a =  [ 0.99955581 -1.99911162  0.99955581] (a0, a1, a2)
b =  [ 0.          1.99911142 -0.99911182] (don't use, b1, b2)

Once the filter gets going (after 1 second) it work great, but it has a wicked response in that first second. I tried initializing the first two y's with both 0's and the future 1 second average. Testing the filter with white noise, here is an example of the initial transient. The initial variance of the white noise is 1 and the transient wanders up to 370 before settling down!

Question: What, if anything, can I do to calm this down?

BTW, my one-pole filter doesn't do this! Is it a stability issue with the two pole filter with a large ratio of sample frequency to cutoff frequency?

See related question: Is initializing a digital filter's output with no "momentum" a non-trivial task?

The easy way is to subtract, prior to filtering, from each input sample the first input sample. This will eliminate the step in the beginning, avoiding the wavy step response. That's roughly what the other original poster said they did. Or because you have noise, you can take the average of the first few input samples as a better estimate of the DC bias and subtract that from each input sample, prior to filtering.

If you want to save the work of subtracting then changing the filter's state variables in the beginning is the way to go, and my answer to the linked question shows how to do that (with $x$ as the estimated DC bias and $y=0$ as you have a high-pass filter with zero DC response but this can be validated). Not sure if you can easily touch the state variables in your prototyping environment without rewriting the filter code.


What Do Skewness and Kurtosis Represent? I understand that the question could mean a lot of things but I am thinking specifically to image processing. For example, I know that the mean can be a basic texture feature that represents the average pixel value of the image. That is, it is useful for determining an image background. Then I know that the standard deviation can be an indicator of how spread is the histogram and is useful for indicating what other pixel values also belong to the background.

In similar terms to this, what do Kurtosis and Skewness represent? I have read about them in various posts but I can't find anything related to what they represent in terms of what I have explained above.

I am not sure if I have came to the right place but if I haven't feel free to point me in the right direction.

This paper has a nice write up on those and many other measures.

For skewness:

  In terms of digital image processing, Darker and
  glossier surfaces tend to be more positively skewed than
  lighter and matte surfaces. Hence we can use skewness in
  making judgements about image surfaces.

This is because skewness measures how "lopsided" the distribution of pixel values are.

For kurtosis:

  In digital image processing kurtosis values are
  interpreted in combination with noise and resolution
  measurement. High kurtosis values should go hand in
  hand with low noise and low resolution.

I'm not sure I agree with this.  Images with moderate amounts of salt and pepper noise are likely to have a high kurtosis value.


Windowed N-point FFT In the previous question,

Smoothing in frequency domain and impulse reponse of the filter in time-domain

I guess I have made a mistake in-terms for applying windowing to FFT.

After reading "Understanding Digital Signal Processing - Richard G. Lyons" few things are rather very clear now. Nevertheless to be certain I would like to ask whether my understanding is correct or not.

Frequency Domain Windowing is simply given by,

$$X_{\mathrm{three-term}}(m)=\alpha X(m) - \frac{\beta}{2}X(m-1)-\frac{\beta}{2}X(m+1)$$
where, $\alpha$ and $\beta$ are the co-efficients of Hamming window or Hann window.

Further
$$X_{\mathrm{three-term}}(0)=\alpha X(0)-\frac{\beta}{2}X(N-1)-\frac{\beta}{2}X(1)$$
and
$$X_{\mathrm{three-term}}(N-1)=\alpha X(N-1)-\frac{\beta}{2}X(N-2)-\frac{\beta}{2}X(0)$$

So after taking fft(X(m)), the windowed N-point FFT is simply the above three equations?

All these formulas are from Understanding Digital Signal Processing - Richard G. Lyons and all the credits goes to the author of the book.

Some questions

Is Windowed FFT called smoothing?
Do I need to do anything else apart from ifft(X_3term(m)) to perform inverse of $X_{\mathrm{three-term}}(m)$?

"Is Windowed FFT called smoothing?" No.
"Do I need to do anything else apart from ifft(X_3term(m)) to perform inverse of Xthree−term(m)" No.
Here is what is going on:
Suppose
$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} nk } $$
This is the conventional definition of the DFT which is unnormalized.
Now do the following with it:
$$
\begin{aligned}
Y[k] &=  \alpha X[k] - \frac{\beta}{2}X[k-1]-\frac{\beta}{2}X[k+1] \\
&= \alpha \sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} nk }
- \frac{\beta}{2} \sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} n(k-1) }
- \frac{\beta}{2} \sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} n(k+1) } \\
&= \alpha \sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} nk }
- \frac{\beta}{2}  \sum_{n=0}^{N-1} e^{ i \frac{2\pi}{N}n} x[n] e^{-i \frac{2\pi}{N} nk }
- \frac{\beta}{2}  \sum_{n=0}^{N-1} e^{-i \frac{2\pi}{N}n} x[n] e^{-i \frac{2\pi}{N} nk } \\
&=  \sum_{n=0}^{N-1} \left[ \alpha - \frac{\beta}{2} e^{ i \frac{2\pi}{N}n} - \frac{\beta}{2} e^{-i \frac{2\pi}{N}n} \right] x[n] e^{-i \frac{2\pi}{N} nk } \\
&=  \sum_{n=0}^{N-1} w[n] x[n] e^{-i \frac{2\pi}{N} nk } \\
\end{aligned}
$$
Where
$$
\begin{aligned}
w[n]  &=\alpha - \frac{\beta}{2} e^{ i \frac{2\pi}{N}n} - \frac{\beta}{2} e^{-i \frac{2\pi}{N}n}  \\
&= \alpha - \beta \cos \left( \frac{2\pi}{N}n \right) \\
\end{aligned}
$$
So
$$ y[n] = \text{DFT}^{-1} ( Y[k] ) =  w[n] \cdot x[n] $$
You might also enjoy:

The Zeroing Sine Family of Window Functions


How to find out amplitude of harmonic after diff type of FT? I am new in DSP, so I apologies if it's dump questions but I didn't find needed information on the Internet.
So the questions is:
Assume that we have simple harmonic signal described by next function $$x(t) = A \cos(2 \pi f_0 t)$$ After Fourier transformation amplitude of initial harmonic in spectrum changed as next:

For FT: $A/2$
For DFT and FFT: $(A \times N)/2$ where $N$ is the number of samples

Am I right about the last one for discrete transformations? And what about DTFT what is amplitude at frequency $f_0$ equal to? Is it equal to $A/(2T)$ where $T$ is sampling period of the signal.

Lets find the spectrum of the signals $x_c(t) = A \cos(\Omega_0 t)$ and $x_d[n] = x_c(nT) = A\cos(\Omega T n)$

I-First consider the CTFT (continuous-time Fourier transform) $X(\Omega)$ of $x_c(t)$:
$$x_c(t) = A \cos(\Omega_0 t) = A\frac{e^{j\Omega_0 t} +e^{-j\Omega_0 t} }{2} \longleftrightarrow X_c(\Omega) = A \pi \delta(\Omega-\Omega_0) + A \pi \delta(\Omega+\Omega_0) $$ The CTFT spectrum involves two impulses at $\pm \Omega_0$. Their magnitudes are infinite whereas their weights are $A \pi$

II-Second consider the DTFT (discrete-time Fourier transform) $X_d(e^{j \omega})$ of $x_d[n] = A\cos(\Omega_0 n T)$

Approach-1: Considering directly the signal $x_d[n] = A\cos(\Omega_0 n T)$ where $\Omega_0 T = \omega_0$ and by using the DTFT properties:

$$x_d[n] = A \cos(\Omega_0 T n) = A\frac{e^{j\omega_0 t} +e^{-j\omega_0 n} }{2} \longleftrightarrow X_d(e^{j \omega}) = A \pi \delta(\omega-\omega_0) + A \pi \delta(\omega+\omega_0 ) $$ The DTFT spectrum involves two impulses (in $-\pi < \omega < \pi$) at $\pm \omega_0 = \Omega_0 T$. Their magnitudes are infinite and their weights are $A \pi$

Approach-2: Considering the relation between DTFT $X_c(\Omega)$ and CTFT $X_d(e^{j\omega})$ due to sampling $$\boxed{ X_d(e^{j\omega}) = \frac{1}{T} X_c(\frac{\omega}{T})}$$

Given that the CTFT $X_c(\Omega) = A \pi [\delta(\Omega-\Omega_0) + \delta(\Omega+\Omega_0)]$ we conclude on the DTFT as:
$$X_d(e^{j\omega}) =  \frac{1}{T} A \pi [\delta(\frac{\omega}{T}-\Omega_0) + \delta(\frac{\omega}{T}+\Omega_0)]$$

Recalling the property of the impulse as:
$$ \delta(\frac{\omega}{T}-\Omega_0) = \delta( \frac{\omega- \Omega_0 T }{T} ) = T \delta( \omega- \Omega_0 T )=T \delta( \omega- \omega_0 )$$

Yields for DTFT $X_d(e^{j\omega})$:
$$X_d(e^{j\omega}) =  \frac{1}{T} A \pi [\delta(\frac{\omega}{T}-\Omega_0) + \delta(\frac{\omega}{T}+\Omega_0)] = A \pi [\delta(\omega-\omega_0) + \delta(\omega + \omega_0)]$$ The same result of approach-1 is observed: The DTFT spectrum involves two impulses (in $-\pi < \omega < \pi$) at $\pm \omega_0 = \Omega_0 T$. Their magnitudes are infinite and their weights are $A \pi$

III-Third consider the $N$-point DFT of the signal $x_d[n] = A \cos(\Omega T n) = A\cos(\omega_0 n)$

Approach-1: Considering directly the $N$-point DFT sum of the discrete-time signal $x[n] = A\cos(\omega_0 n)$ for $n=0$ to $n=N-1$.

(NOTE for the simple analysis we shall assume there's an integer $k_0$ such that $\Omega_0 T = \omega_0 = \frac{2 \pi}{N} k_0$ is satisfied; i.e., the frequency $\omega_0$ resides exactly on a DFT bin frequency for the index $k = k_0$, for $0 \leq k < N$.)

$$X[k] = \sum_{k=0}^{N-1} A \cos(\omega_0 n) e^{-j\frac{2 \pi}{N} k n } = \frac{A}{2} \sum_{k=0}^{N-1}[e^{j\omega_0 n}+e^{-j\omega_0 n}]  e^{-j\frac{2 \pi}{N} k n }  $$ Since we assume that $\omega_0 = \frac{2\pi}{N} k_0$ then we get:

$$X[k] = \frac{A}{2} \sum_{k=0}^{N-1}[e^{j\frac{2\pi}{N} k_0 n}+e^{-j\frac{2\pi}{N} k_0 n}]  e^{-j\frac{2 \pi}{N} k n }   = \frac{A}{2} \sum_{k=0}^{N-1}e^{j\frac{2\pi}{N} (k_0 -k) n} + \frac{A}{2} \sum_{k=0}^{N-1}e^{-j\frac{2\pi}{N} (k_0 +k) n}$$

We can evaluate the sums using the relation
$$\sum_{k=0}^{N-1}e^{j\frac{2\pi}{N} (k-k_0) n} =
\begin{cases}
N  &, \text{ for } k = k_0 \\
0 &, k \neq k_0 \\
\end{cases}$$
Which yields for DFT $X[k]$ as:

$$X[k] = \begin{cases}
0 &, \text{ for } 0 \leq k < k_0 \\
\frac{A}{2} N &, \text{ for } k = k_0 \\
0 &, \text{ for } k_0 < k < N-k_0 \\
\frac{A}{2} N &, \text{ for } k = N-k_0 \\
0 &, \text{ for } N-k_0 < k < N-1 \\
\end{cases}$$

Approach-2 considering that the finite length signal $x_v[n]$ is obtained from the infinite length signal $x_d[n] = A\cos(\omega_0 n)$ by multiplying it with a window $w[n]$ of length $N$. And considering its effect on the DTFT of the windowed signal $x_v[n] = x_d[n] w[n]$ and considering that DFT $X_v[k]$ is obtained as the uniform samples of the DTFT $X_v(e^{j\omega})$ of $x_v[n]$:

Since $x_v[n] = x_d[n] w[n]$ then;
$$X_v(e^{j\omega}) = \frac{1}{2\pi} X_d(e^{j\omega}) \star W(e^{j \omega})$$

We previously computed $X_d(e^{j \omega})$ as:
$$ X_d(e^{j \omega}) = A \pi \delta(\omega \pm \omega_0)$$

And it can be shown that the DTFT of the window is:
$$W(e^{j \omega})= e^{-j \omega ( \frac{N-1}{2} ) } \frac{ \sin(\omega N/2) }{\sin(\omega /2) }$$

And the DTFT $X_v(e^{j \omega})$ becomes:
$$ X_v(e^{j \omega}) = \frac{1}{2\pi} X_d(e^{j \omega}) \star W(e^{j \omega}) = \frac{A}{2} W(e^{j (\omega-\omega_0)}) + \frac{A}{2} W(e^{j (\omega+\omega_0)}) $$

Finally the DFT $X_v[k]$ of $x_v[n]$ will be the samples of this $V(e^{j \omega})$ at the frequencies $w_k = 2\pi k /N$. Under this special sampling we have the following:
$$X_v[k] = X_v(e^{j \omega})|_{w_k = \frac{2 \pi}{N} k} = \frac{A}{2} W(e^{j (\frac{2 \pi}{N} k-\omega_0)}) + \frac{A}{2} W(e^{j (\frac{2 \pi}{N} k + \omega_0)}) $$

and replacing again $\omega_0$ with $\frac{2 \pi}{N} k_0$ we have:
$$X_v[k] = \frac{A}{2} W(e^{j (\frac{2 \pi}{N} (k - k_0))}) + \frac{A}{2} W(e^{j (\frac{2 \pi}{N} (k + k_0)}) $$

Similar to approach-1, and considering $W(0) = N$ we have
$$
W[k-k_0] = \begin{cases}
0 &, \text{ for } 0 \leq k < k_0 \\
N &, \text{ for } k = k_0 \\
0 &, \text{ for } k_0 < k < N-k_0 \\
N &, \text{ for } k = N-k_0 \\
0 &, \text{ for } N-k_0 < k < N-1 \\
\end{cases}
$$

And therefore we conclude for DFT $X_v[k]$ as
$$
X_v[k] = \begin{cases}
0 &, \text{ for } 0 \leq k < k_0 \\
\frac{A}{2}N &, \text{ for } k = k_0 \\
0 &, \text{ for } k_0 < k < N-k_0 \\
\frac{A}{2}N &, \text{ for } k = N-k_0 \\
0 &, \text{ for } N-k_0 < k < N-1 \\
\end{cases}
$$


understanding discrete-time convolution in LTI systems I'm trying to understand the discrete-time convolution for LTIs and its graphical representation. standard explanations (like: this one) start with the idea of decomposing an input signal $x[t]$ into a sum of unit impulses. We have unit impulse function $\delta$:
$$
\delta[n]=
\begin{cases}
    1,& \text{if}\ n = 0\\
    0,& \text{if}\ n\neq0
\end{cases}
$$

then $x[t]$ can be represented as:

$$x[t] = \sum_{k = -\infty}^{+\infty}x[k]\cdot\delta[t-k]$$

so the input signal at $t$, $x[t]$, can be represented this way because $\delta[t-k]$ is going to be 1 only if $t = k$, and 0 elsewhere. this makes great sense.

the confusion for me begins in Transparency 4.2 of above link. the left hand side makes perfect sense, but what does the right represent? how does one unit impulse generate a sequence of responses? what does $x[0]h[n]$ represent here? How do you go from this to Transparency 4.9 where $h$ is used instead of $\delta$?

put differently: can someone explain the intuition behind $h$ (the response of a an LTI system) being reflected and time shifted? in the above trick for representing simply the input signal $x[t]$, there's no time shifting really; $\delta[t-k]$ seems like just a mathematical trick for getting a $1$ in the right place and zero elsewhere; it's not intuitive for me to think of it as time shifted.

if we use the fact that in LTI, the response of a system to input $n$, $y[n]$, is just the superposition of the responses to unit impulses of input signal $x[n]$, why can't we write convolution this way instead?

$$y[n] = \sum_{k=-\infty}^{+\infty}x[k]h[n]$$

You know that $h(n)$ is the system's response to the input $x(n)=\delta(n)$. Because of time-invariance, the response to $x(n)=\delta(n-k)$ is given by $h(n-k)$. If you write the input signal as

$$x(n)=\sum_kx(k)\delta(n-k)\tag{1}$$

then, due to linearity, you can write the output signal as the sum of the responses to each individual term $x(k)\delta(n-k)$ in (1), which is $x(k)h(n-k)$. So you get the convolution sum

$$y(n)=\sum_k x(k)h(n-k)\tag{2}$$

To compute (2) for a given value of $n$ you need to multiply $x(k)$ with $h(n-k)$ for all relevant values of $k$. Since you have $h(n\color{red}{-k})$ in (2) (i.e. with a minus sign before $k$) you need to flip the impulse response (shifted by the current value of $n$), multiply it with $x(k)$ and compute the sum.


Fourier series - time and frequency domain confusion I am computing the fourier series of the following function between $[-0.5, 0.5]$
$$\displaystyle f(t) = \frac{1}{2} - |t|$$
According to the definition of Fourier Series the coefficients are given by
$$\displaystyle \hat f_n = \frac{-(-1)^n}{2\pi ^2 n^2} + \frac{1}{2\pi ^2 n^2}$$
That's easy to calculate. So far so good.

With the coefficients I know I can do the approximation of the function $f(t)$. But then how to plot, the frequency spectrum??

For example, I defined the function $f(t)$ in matlab and plot it from $-0.5$ to $0.5$ and got the triangle function. Now, I have the coefficients that depend on $n$ but what's the range of $n$? Plotting the coefficients will give me the frequency spectrum?

Please help me clarify this I'm really confused, seems I can't get the relation.

The coefficients of the Fourier series that you have computed are, in effect, the spectrum of the periodic signal consisting of the sum of signals $f(t)$ delayed in time or advanced in time by integer multiples of $1$ second.  Mathematically, the Fourier transform of a periodic function has impulses in it with the impulse amplitudes being the Fourier series coefficients. The impulses occur at intervals of the period, in this instance $1$ second,
and a "plot" of the spectrum typically shows a bunch of vertical arrows (a common representation of impulses) with the amplitude marked just above the arrowhead. In your case, assuming that you have computed the complex form of the Fourier series,
you would get a series of such arrows extending from $n=0$ out to $\pm \infty$.
Note that your formula
$$
 \hat f_n = \frac{-(-1)^n}{2\pi ^2 n^2} + \frac{1}{2\pi ^2 n^2}
$$
is saying that $\hat{f}_n = 0$ when $n$ is even. Also, putting $n=0$ in the formula
makes $\hat{f}_0$ undefined which is incorrect since $\hat{f}_0$ is the
"DC value" which is $\frac{1}{4}$ by inspection.. So please do check your work.
In particular, note that while integrating by parts as you are likely
to be doing, the antiderivative
$$\int e^{-j2\pi nt}\,\mathrm dt = \frac{e^{-j2\pi nt}}{-j2\pi n}$$
cannot be used when $n=0$: $\int 1 \,\mathrm dt = t$.

Note: if you computed the sine/cosine form of the Fourier series, then you
should have found that the coefficients of the $\sin(2\pi nt)$ terms are all
$0$ and your spectrum will extend from $n=0$ to $n=\infty$ only.


changing angular data to remove the sharp changes I have a 200 by 200 matrice named theta which shows an angle in a 200 by 200 box. The angle is in the interval [-90,90]. When I plot theta in a 2d plane, there are some positions in which the angle jumps from -90 to 90. This is shown for a particular choice case in the below picture:

As it is clear from the picture, there are positions where theta jumps from -80 to 80. I want to make theta even by adding (or subtracting) 180 to it in the positions with sharp change. The reason that I can only add (or subtract) 180 is that in the problem that I am trying to solve, has a symmetry which allows adding 180 degrees to the data. I appreciate it if someone could solve this problem.

The script that I use to plot it is:

[a, b] = contourf(theta, 88);
    set(b,'EdgeColor','none');
    colormap(summer);
    shading flat;
    colorbar;

Using the proposed answer I get the following:

Could you do something like: (1/2)*unwrap(2*theta)?

edit:
More concretely, I think about this. I would think that the numerical example and asserts are sufficient to see if it can be adapted to your situation.

It is standard usage of the unwrap() function to replace large jumps with the "phase alias" that cause the smallest jump in phase, applied to a 2-d input. As unwrap expects radians rather than degrees, I have to convert back and forth, and as your phase "aliases" at 180 degrees rather than 360, I am scaling by 2 and 1/2.

th = repmat([70 80 -90 -80], 2, 1)
th =

70    80   -90   -80
70    80   -90   -80

th2 = 2*th;
deg_to_rad = 2*pi/360;
th_unwrapped = (1/deg_to_rad)*((1/2)*unwrap(th2*deg_to_rad, [], 2))
th_unwrapped =

70    80   90   100
70    80   90   100
assert(all(rem(th-th_unwrapped, 180)==0, [1 2]))
assert(all(diff(th_unwrapped, [], 2)==10, [1 2]))


Extract voice from magnetic tape There is a recording which contains very noisy voice. (voice
  with a low Signal-to-Noise ratio).

It is a female voice of a stationary speaker (right in front of the microphone) and there are many other recordings of the same
  speaker with much better quality.

What can be done to denoise the audio?

From the other recordings, it might be possible to extract statistical values that allow some model to filter the speaker's voice.

Existing questions only address snoring, have non-random noise, want to identify a specific artifact or a moving speaker. The question is-it-possible-to-filter-out-a-persons-voice-out-of-100-of-other-voices is similar, but the answer there lacks details and is more about distinguishing, less about extraction.

Some of the tapes are direct copies, some are copies of copies, which have also been recorded using a lower speed setting (4.75 cm/s). This affected the quality prior to digitalization.

I am afraid there's not much you can do.  The voice part seems to have gone through the equivalent of a low pass filter with a cut off of around 1000 Hz.  Basically, all of the speech components above 1000Hz are gone.

The filtering action may not have been an intentional filter, but may have been due to the improper biasing of the tape during recording.  If it is an old tape, it may simply have deteriorated over time.  Also, the playback head may have needed degaussing.

Running it through a high pass flattens the frequency response, but doing so pushes the level down so far it drowns the signal in the noise - no help there.

The best result I got was from using a very steep low pass with a cutoff of 1000Hz together with a very steep high pass with a cutoff of 160 Hz.  That gets rid of the noise by only passsing what is left of the actual speech, but obviously it can't recover what was lost.

Your real problem is not the noise, it is the lost frequency range.

This is the spectrum of the bad recording:

This is the spectrum of the good recording:

As you can see, there's a lot of stuff missing from the bad recording.  So, it isn't simply a problem of removing noise.  The problem is that there's stuff that's just GONE.

Look at the range from 1000Hz to 7000Hz.  There's lots of stuff there in the good recording, but in the bad one it is just a flat spectrum a good 30dB below the voice peaks aroudn 400Hz.

Some of what's missing might be buried in the noise, but recovering it would cause artifacts that are worse than the noise and muffled sound.

Looking at just the noise, it doesn't seem like there would be much to recover out of it.  It looks just like the portions with speech (except for between 160Hz to 1000Hz,) so anything that is in there is going to be buried really deep.


Short term energy signal explanation I'm trying to calculate the short term energy of a digital signal
in python, using a rectangular window, but don't understand
the basics of how this works, a pointer to any good explanation would be great.

Specific things I'm confused about:

Where do I initialise the window? (at position 0,window_length,
or -window_length,0, and then work through the signal? ...)
If I sum up the squares of all samples at every position inside the window,
don't I end up with a 'square signal', where all the values are the
same for the length of the window step?

Thanks.

For initializing window you need the start point would depend on your overlap. If the overlap is 50% then you would start at window_length/2
If you are using a rectangular window and with 50% overlap in that case you want to use a square root rectangular window on the input and output so that there are no scaling issues. so basically you window will be scaled also dependent on your overlap.

Also it would be worthwhile to consider using other windows dependent on your application. As in a lot of other windows the scaling is taking into account.