What Properties Make Certain Wavelets "Better" Than Others In Image Compression?

Overview

The short answer is that they have the maximum number of vanishing moments for a given support (i.e number of filter coefficients). That's the "extremal" property which distinguishes Daubechies wavelets in general. Loosely speaking, more vanishing moments implies better compression, and smaller support implies less computation. In fact, the tradeoff between vanishing moments and filter size is so important that it dominates the way that wavelets are named.  For example, you'll often see the D4 wavelet referred to either as D4 or db2.  The 4 refers to the number of coefficients, and the 2 refers to the number of vanishing moments. Both refer to the same mathematical object. Below, I'll explain more about what moments are (and why we want to make them disappear), but for now, just understand that it relates to how well we can "fold up" most of the information in the signal into a smaller number of values. Lossy compression is achieved by keeping those values, and throwing away the others.

Now, you may have noticed that CDF 9/7, which is used in JPEG 2000, has two numbers in the name, rather than one. In fact, it's also referred to as bior 4.4. That's because it's not a "standard" discrete wavelet at all. In fact, it doesn't even technically preserve the energy in the signal, and that property is the entire reason people got so excited about the DWT in the first place! The numbers, 9/7 and 4.4, still refer to the supports and vanishing moments respectively, but now there are two sets of coefficients that define the wavelet. The technical term is that rather than being orthogonal, they are biorthogonal. Rather than getting too deep into what that means mathematically, I'll just review the factors which led to using non-energy-preserving biorthogonal wavelets in the first place.

JPEG 2000

A much more detailed discussion of the design decisions surrounding the CDF 9/7 wavelet can be found in the following paper:

  Usevitch, Bryan E. A Tutorial on Modern Lossy Wavelet Image
  Compression: Foundations of JPEG 2000.

I'll just review the main points here.

Quite often, the orthogonal Daubechies wavelets can actually result in increasing the number of values required to represent the signal. The effect is called coefficient expansion. If we're doing lossy compression that may or may not matter (since we're throwing away values at the end anyway), but it definitely seems counterproductive in the context of compression. One way to solve the problem is to treat the input signal as periodic.
Just treating the input as periodic results in discontinuities at the edges, which are harder to compress, and are just artifacts of the transform. For example, consider the jumps from 3 to 0 in the following periodic extension: $[0,1,2,3] \rightarrow  [...0,1,2,3,0,1,2,3,...]$.  To solve that problem, we can use a symmetric periodic extension of the signal, as follows: $[0,1,2,3] \rightarrow [...,0,1,2,3,3,2,1,0,0,1...]$. Eliminating jumps at the edges is one of the reasons the Discrete Cosine Transform (DCT) is used instead of the DFT in JPEG. Representing a signal with cosines implicitly assumes "front to back looping" of the input signal, so we want wavelets which have the same symmetry property.
Unfortunately, the only orthogonal wavelet which has the required characteristics is the Haar (or D2, db1) wavelet, which only as one vanishing moment. Ugh. That leads us to biorthogonal wavelets, which are actually redundant representations, and therefore don't preserve energy. The reason CDF 9/7 wavelets are used in practice is because they were designed to come very close to being energy preserving. They have also tested well in practice.

There are other ways to solve the various problems (mentioned briefly in the paper), but these are the broad strokes of the factors involved.

Vanishing Moments

So what are moments, and why do we care about them? Smooth signals can be well approximated by polynomials, i.e. functions of the form:

$$a + bx + cx^2 + dx^3 + ...$$

The moments of a function (i.e. signal) are a measure of how similar it is to a given power of x.  Mathematically, this is expressed as an inner product between the function and the power of x.  A vanishing moment means the inner product is zero, and therefore the function doesn't "resemble" that power of x, as follows (for the continuous case):

$$\int{x^n f(x) dx = 0 }$$

Now each discrete, orthogonal wavelet has two FIR filters associated with it, which are used in the DWT. One is a lowpass (or scaling) filter $\phi$, and the other is a highpass (or wavelet) filter $\psi$. That terminology seems to vary somewhat, but it's what I'll use here. At each stage of the DWT, the highpass filter is used to "peel off" a layer of detail, and the lowpass filter yields a smoothed version of the signal without that detail. If the highpass filter has vanishing moments, those moments (i.e. low order polynomial features) will get stuffed into the complementary smoothed signal, rather than the detail signal. In the case of lossy compression, hopefully the detail signal won't have much information in it, and therefore we can throw most of it away.

Here's a simple example using the Haar (D2) wavelet. There's typically a scaling factor of $1/\sqrt{2}$ involved, but I'm omitting it here to illustrate the concept. The two filters are as follows:
$$ \phi = [1,1] \\  \psi = [1,-1]$$

The highpass filter vanishes for the zero'th moment, i.e. $x^0 = 1$, therefore it has one vanishing moment. To see this, consider this constant signal: $[2,2,2,2]$.  Now intuitively, it should be obvious that there's not much information there (or in any constant signal).  We could describe the same thing by saying "four twos". The DWT gives us a way to describe that intuition explicitly. Here's what happens during a single pass of the DWT using the Haar wavelet:

$$
[2,2,2,2] \rightarrow_{\psi}^{\phi} \left\{ \begin{array}{rr}
\left[2 + 2, 2 + 2\right] = \left[4,4\right] \\
\left[2-2,2-2\right] = \left[0,0\right]
\end{array}\right.
$$

And what happens on the second pass, which operates on just the smoothed signal:

$$
[4,4] \rightarrow_{\psi}^{\phi} \left\{ \begin{array}{rr}
\left[4 + 4\right] = \left[8\right] \\
\left[4-4\right] = \left[0\right]
\end{array}\right.
$$

Notice how the constant signal is completely invisible to the detail passes (which all come out to be 0).  Also notice how four values of $2$ have been reduced to a single value of $8$.  Now if we wanted to transmit the original signal, we could just send the $8$, and the Inverse DWT could reconstruct the original signal by assuming that all the detail coefficients are zero. Wavelets with higher-order vanishing moments allow similar results with signals that are well approximated by lines, parabolas, cubics, etc.

Further Reading

I'm glossing over a LOT of detail to keep the above treatment accessible. The following paper has a much deeper analysis:

  M. Unser, and T. Blu, Mathematical properties of the JPEG2000 wavelet
  filters, IEEE Trans. Image Proc., vol. 12, no. 9, Sept. 2003,
  pg.1080-1090.

Footnote

The above paper seems to suggest that the JPEG2000 wavelet is called Daubechies 9/7, and is different from the CDF 9/7 wavelet.

  We have derived the exact form of the JPEG2000 Daubechies 9/7 scaling
  filters... These filters result from the factorization of the same
  polynomial as $Daubechies_{8}$ [10]. The main difference is that the
  9/7 filters are symmetric. Moreover, unlike the biorthogonal splines
  of Cohen-Daubechies-Feauveau [11], the nonregular part of the
  polynomial has been divided among both sides, and as evenly as
  possible.

  [11] A. Cohen, I. Daubechies, and J. C. Feauveau, “Biorthogonal bases
  of compactly supported wavelets,” Comm. Pure Appl. Math., vol. 45, no.
  5, pp. 485–560, 1992.

The draft of the JPEG2000 standard (pdf link) that I've browsed also calls the official filter Daubechies 9/7. It references this paper:

  M. Antonini, M. Barlaud, P. Mathieu, and I. Daubechies, “Image coding
  using the wavelet transform,” IEEE Trans. Image Proc. 1, pp. 205-220,
  April 1992.

I haven't read either of those sources, so I can't say for sure why Wikipedia calls the JPEG2000 wavelet CDF 9/7. It seems like there may be a difference between the two, but people call the official JPEG2000 wavelet CDF 9/7 anyway (because it's based on the same foundation?).  Regardless of the name, the paper by Usevitch describes the one that's used in the standard.


Are the DFT/FFT multiplications in the fourier domain complex multiplications?

Yes, this is correct - the DFT/FFT of a real signal (for example an image) is not necessarily real; so if you want to compute a convolution through multiplications in the Fourier domain, the coefficients to multiply will be complex numbers; and will thus cost 4 real multiplications and 2 real additions if a naive implementation is used (unless the architecture you are writing code for has hardware support for complex multiplications).


how to find sample rate of sum of 2 cos signals

Your signal has maximum frequency component of 10Hz which means a sampling rate of 20Hz according to Nyquist theorem should suffice. Then the cut-off frequency of 5Hz in the ideal filter should be enough to retain s(t) and eliminate v(t). However in the digital domain, typically a much higher sample rate is needed to process and display continuous time signals. I used a sampling rate of 100Hz in a simple python script below. Cutoff frequency of the selected FIR filter is still 5Hz to match s(t). It is customary to normalize sampling frequency to \pi in which case cut-off frequency becomes \pi/20 instead of \pi/4 as you showed for the theoretical case.
Below python script:
import numpy as np
from scipy.signal import firwin, lfilter
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt

# sampling frequency in Hz
Fs = 100

# center frequencies of a linearly combined signal
fc1 = 5
fc2 = 10

# Observe duration in seconds, forumulate a time-axis
T = 5
tvals = np.arange(0, T, 1/Fs)

fig=plt.figure(figsize=(10,10))
s_t = np.cos(2*np.pi*fc1*tvals)
v_t = np.cos(2*np.pi*fc2*tvals)

plt.subplot(3,1,1)
plt.plot(tvals,s_t)
plt.ylabel('low-pass signal')
plt.grid()
x_t = s_t + v_t

plt.subplot(3,1,2)

plt.plot(tvals, x_t)
plt.ylabel('combined signal')
plt.grid()

# Filter Length
N = 31
taps = firwin(N, fc1/Fs)

# Simple FIR filtering
filtered_signal = lfilter(taps, 1.0, x_t)
ll = len(filtered_signal)
plt.subplot(3,1,3)

# Plot final result
plt.plot(tvals[int((N -1)/2):ll], filtered_signal[int((N - 1)/2):ll])
plt.grid()
plt.xlabel('Time')
plt.ylabel('filtered signal')
plt.tight_layout()
plt.show()

will generate a figure like:
Observe the recovered signal after filtering is the same (minus the filter delay effect) as the s(t)


Extracting feature from audio track using DWT

Four features is a very short description for an audio signal. The hint: "Write a program to transform an audio track to time-frequency domain" is important. Since audio signals are in general not fully stationary, windowed-Fourier or Short-term Fourier transforms can be important. Such representations are generally redundant, yet one can expect that the most energetic features in those domains represent the main energy in the time domain.

One first bet is to use a time-frequency spectrogram, and pick the most important features, or the highest local maxima.

Then, if needed, you can think about reducing the redundancy with filter-banks or discrete wavelet transforms. You can further  incorporate auditory system weighting, or classification priors.


How to check if a signal is power signal or energy signal?

let's assume the signal, $x(t)$ is not identically zero for all $t$.
an "energy signal" (what i would prefer to call a "finite energy signal") is such a signal, $x(t)$ with a finite energy:
$$ 0 \ < \ \int_{-\infty}^{\infty} |x(t)|^2 \ dt \ < \ +\infty $$
BTW, sometimes for mathematical ease, we require a stricter sense of finite "energy":
$$ 0 \ < \ \int_{-\infty}^{\infty} |x(t)| \ dt \ < \ +\infty $$
and a "power signal" (what i would prefer to call a "finite power signal") is such a signal, $x(t)$ with finite power:
$$ 0 \ < \ \lim_{T \to +\infty} \frac{1}{T}\int_{-\frac{T}2}^{\frac{T}2} |x(t)|^2 \ dt \ < \ +\infty $$
i think that is the most fundamental definitions of the two classes of continuous-time signals.  you can do a very similar definitions for discrete-time signals, $x[n]$.


Why is CDMA so complicated, wouldn't this simpler solution work?

Your solution works, and it is known by the name "time division multiple access" (TDMA). Basically, you assign time slots to each transmitter.

CDMA has several advantages compared to TDMA, but it also has its own disadvantages. Both schemes are inferior to more modern techniques such as OFDMA.

Main disadvantages of TDMA:

With data rate constant, the time slices become narrower as the number of transmitters increases. This means that the bandwidth also increases. With CDMA, you can increase the number of transmitters without increasing the bandwidth, by adding orthogonal code words.
TDMA requires very precise synchronization between transmitters. It is very hard to keep a large number of mobile units synchronized so that they won't transmit during the time slot of another unit. This may be solved using guardbands, but that decreases efficiency.

A couple of advantages of CDMA:

More resistance to fading and multipath. CDMA reflections are orthogonal to each other, so they cancel out.
If the codes are private, then there's a measure of inherent privacy.


This line seems suspicious to me:

f = fs/2*linspace(0,1,nfft/2+1);

Your frequency bins should be:

$$f_k = \frac{f_s}{\tt nfft} k$$

so

$$ \begin{array}{rcl}
f_{{\tt nfft}/2 + 1} &=& \frac{f_s}{\tt nfft} ({\tt nfft}/2 + 1)\\
&=& \frac{f_s}{2} + \frac{f_s}{\tt nfft}
\end{array}$$

whereas your linspace setting has it mapping to just $\frac{f_s}{2}$.


Show that $S[n]= \frac1M \sum e^{\frac{2\pi l n}M}$ for $l = 0$ to $M - 1$

HINT:

First of all, Eq. $(70)$ is wrong, there's an imaginary unit $j$ missing in the exponent. If you add that missing $j$, note that for $n=kM$, $k\in\mathbb Z$, you have $e^{j2\pi ln/M}=1$, so the result of the scaled sum equals $1$. In order to see that the sum equals zero for all other values of $n$, use the formula of the geometric series.

This answer is also relevant to your question.


What is the Z-transform of $0.8^{n+2}u(n-1)$?

HINT:
If you rewrite $y[n]$ as
$$y[n]=(0.8)^3(0.8)^{n-1}u[n-1]\tag{1}$$
does it become easier to find the $\mathcal{Z}$-transform?


Is there an adjective describing a filter with kernel that has zero mean?

"Zero-Mean" is the word that's commonly used to describe signals and signals with a zero average. "This is a zero-mean filter."
If you really mean a filter that is specifically meant to cancel the DC component, a "DC blocker" is a name for that.


Random (Over) Sampling signal and perfect Reconstruction in Nyquist form?

For the moment restrict the discussion to 1D signals. Sampling operation can be performed as uniformly or non-uniformly. The uniform sampling is the most obvious, simplest and the preferred form, unless otherwise stated.

Nonuniform sampling can be performed in a number of ways, periodically or aperiodically, sampling instant given by a formula or being completely random.

Depending on the nonuniform sampling form, the reconstruction criterions may change. For example when you have a formula based nonuniform sampling instant $t_n$, you can (if possible) find a kernel $\phi_n(t)$ for the reconstruction integral that provides exact (perfect) recovery of a continuous time signal $x(t)$ from its nonuniform samples $x_n$. In this context the $sinc()$ kernel is a special case of the generalized sampling theorem WKS (Whittaker-Kotelnikov-Shannon) with sampling instants given by $t = n T_s$. This was historically called as the Nyquist sampling, who most successfully emplyed it into the earliest digital pulse communication systems.

The primary requirement that Nyquist (uniform) sampling theorem would provide exact reconstruction of a bandlimited baseband signal $x(t)$ from its uniform samples is that the sampling rate satisfy the Nyquist criterion; $F_s \geq 2 F_c$ where $F_c$ is the signal's bandwidth.

Now for the non-uniform sampling strategies, the most common restriction for perfect signal recovery is that on the average the Nyquist rate associated with the bandlimited signal is maintained. So when there is a deficiency of samples at some short interval $I_1$ then there must be an excess of samples at another closeby short interval $I_2$ so that the average Nyquist rate preserved.

The above condition can be stated as the local Nyquist rate criterion.

Note that unlike uniform sampling case, for nonuniform sampling strategies and exlcusively true for the random sampling strategy, the sampling times $t_n$ in addition to the samples $x_n$ should be known in order to reconstruct the signal perfectly.

You can look for Whittaker's, Yen's, Shannon's papers for classical view on the nonuniform sampling of signals.


What is the unit of x-axis of a sound signal?

The below signal has a duration of 3 second and sampling frequency 44100. So should i write Time(s) for the unit of x-axis.

No. In your graph, the X axis seems to be simply the sample number. I mean, it's a bit strange to ask us, considering you generated the graph, but if you benefit from having it spelled out:

In DSP, signals are represented by a sequence of samples. And if you plot these samples, there's no inherent "units" on the X axis. Your signal has been reduced to a sequence of digital numbers, and the fact that there was a sampling rate that implies some period of time between each sample is simply not there in DSP.


it's not "modern", but if done right, correlation methods are quite robust.  i personally like Average Squared Difference Function (ASDF)
$$ Q_x[n_0,k] = \frac{1}{2L+1} \sum\limits_{n=n_0 - \lfloor k/2 \rfloor - L}^{n_0 - \lfloor k/2 \rfloor + L} \left( x[n]-x[n+k] \right)^2 $$

but the more familiar Average Magnitude Difference Function (AMDF)
$$ Q_x[n_0,k] = \frac{1}{2L+1} \sum\limits_{n=n_0 - \lfloor k/2 \rfloor - L}^{n_0 - \lfloor k/2 \rfloor + L} \Big| x[n]-x[n+k] \Big| $$
might be more to your liking.

$x[n_0]$ is the middle of the neighborhood of samples where you are determining the pitch.  $2L+1$ is the number of terms being averaged.  you probably do not need to divide by $2L+1$ if you don't want to, since the values of $Q_x[n_0,k]$ are relative to each other. $\lfloor k/2 \rfloor = k/2$ for even $k$ and $\lfloor k/2 \rfloor = (k-1)/2$ for odd $k$.

look for a value of $k$ that minimizes $Q_x[n_0,k]$.  you will likely want to interpolate around that integer value of $k$ to find a fractional value that is the "true" minimum.  i usually just use simple quadratic interpolation. if $Q_x[n_0,k_m]<Q_x[n_0,k_m-1]$ and $Q_x[n_0,k_m] \le Q_x[n_0,k_m+1]$ then the interpolated minimum (and a candidate for the period) is

$$ P_m = k_m + \frac{1}{2} \frac{Q_x[n_0,k_m+1] - Q_x[n_0,k_m-1]}{2Q_x[n_0,k_m] - Q_x[n_0,k_m+1] - Q_x[n_0,k_m-1]} $$

if the input $x[n]$ is extremely periodic, there might be many values of $k$ that locally minimize $Q_x[n_0,k]$ because

if $x[n+P]=x[n] \quad \forall n$

then $x[n+2P]=x[n]$ or $x[n+3P]=x[n]$

and $Q_x[n_0,P]$ or $Q_x[n_0,2P]$ or $Q_x[n_0,3P]$ are all equally minimum (and close to zero). so then usually you want to pick the minimum that has the smallest $k$.


How can I synchronize signal from an LED and TTL pulses better?

What you are probably missing here is that the frame rate of a webcam can be all over the place, especially if you are capturing video by getting a handle to the imaging device and then requesting frames by it using a timer.

In any case, wherever you are capturing the frame, capture its timestamp as well. This will allow you to cross reference the frames with the voltage time series better.

One of the first things that you have to do is synchronise at least the first pulse. So, independently of what you are going to do in subsequent processing steps, the first two pulses (one from the voltage time series and one from the video time series) must be synchronised. To do this, you simply shift one (or the other) time series so that the rising edges of the pulses coincide.

Now, extract the same amount of waveform in terms of seconds from both time series. So, extract for example 5 seconds from each. One is going to be 150 samples long the other 150000 samples long. Always subsample the more detailed one.

The two waveforms now will be as synchronised as your equipment allows them to be. If you plot them side by side you might find some mismatches because your camera's frame rate has some jitter.

The only way to escape this is by "counteracting" the FPS jitter by dropping voltage time series samples so that the voltage time series (which has higher resolution) "catches up" (or actually, lags behind) with the video.

You can do this either by matching each individual pulse's width or by Dynamic Time Warping (DTW).

The first method is straightforward: Extract a pulse period from both waveforms, measure the pulse width (in either waveform it has to be 150ms, 4.5 frames in the video waveform, 4500 samples in the voltage waveform), downsample the voltage waveform by a factor equal to the ratio between the pulse widths, move on to the next pulse period.

The second method basically does the same thing and can even return a signal which is exactly this "timing advance" signal by which you have to transform one of the waveforms so that it catches up with the other. In this case, things are even more straightforward: Downsample the voltage to the video sampling frequency and this will give you an approximate matching. Then, with both waveforms in the same sampling frequency, run them past DTW and apply the transformation.

Hope this helps.


FIR filter design and implementation - sample rate and number of taps

I'm unable to visualize the real-valued inverse-DFT FIR filter coefficients you describe in your Step 5. But no matter. My answer for A: When you use a large number of DFT bins to yield fine granularity of your frequency-domain representation of your desired frequency response, the greater will be the number of nonzero coefficients remaining after you window the output of your inverse DFT. And the larger the number of nonzero coefficients, the greater will be the number of filter computations needed per filter input sample. Only you can decide how many multiplications and additions you're willing to perform per filter input sample. By the way, I'm wondering how many nonzero coefficients you're contemplating at this time.

My answer for B: I'm not sure what the asterisk symbol means in your y(n) = SUM_over_m[h(m)*x(n-m)] expression. In any case, you typically implement your FIR filter by convolving your filter's input sequence with your filter's nonzero coefficients. If your have 100 nonzero coefficients then you must perform 100 multiplies and 99 additions per filter input sample. (Depending on your computing hardware, there may be a way to reduce the number of multiplies per input sample to 50. Ask me how.) Now if you have many hundreds of nonzero coefficients then in that case, yes, you'll use what's called the "fast convolution" filtering method which involves both forward- and inverse-DFTs. But CA Charlie, all of this discussion depends on whether your filtering is "real-time" of "offline." As for the 48 kHz data rate, it affects two things: (1) it means your input signal's spectrum can contain no spectral components greater than 48/2 = 24 kHz; and (2) in a real-time implementation you must perform all your "per sample" filtering computations in less than 1/48,000th of a second. (Sorry CA Charlie, I don't understand the question in your final paragraph.)


real impulse response

An impulse response that is real-valued.

effect of such a response on ROC

nothing specific, since realness doesn't imply any kind of stability/convergence, far as I can tell. It does tell you something about poles and zeros always appearing in conjugate pairs, but that's it.

restrict the possible area of the poles and zeros

no


The "dimensions" of the spectrogram are not chosen based on where will the spectrogram be fed to but rather depend on your application. Therefore, it is key to understand the spectrogram itself first, as a means of generating features for one or more signals and to an extent, understand the Discrete Fourier Transform (DFT) as well, which is the key operation the spectrogram is based on. In fact, the "dimensions" you are referring to stem entirely from the way the DFT works and determine its spectral "resolution".

If you were to run DFT on some discrete time signal ($x[n]$), of maximum length $N$, you would get back a set of coefficients $X[k]$ that decompose the signal to a weighted sum of sinusoids at different frequencies.

The range of physical frequencies in $x[n]$ is $\left[ 0 \ldots \frac{F_s}{2} \right]$. However, the DFT decomposes the signal into $N_{FFT}$ frequency bins (indexed by $k$) with $N_{FFT} \le N  $. Therefore, the physical frequency that is associated with the $k^{th}$ frequency bin is $f_k = F_s \cdot \frac{k}{N_{FFT}}$.

So what?! So, if you are creating a classifier that operates on some $x_m[n]$ sample obtained from the $m^{th}$ speaker and is trying to figure out the gender of the speaker via the range of frequencies audible, then you would have to tune your $N_{FFT}$ in such a way as to be able to assign different frequencies between 85Hz to 225Hz to different frequency bins. If you did not do that, then your spectral features (the set of $k$ frequency bins) would be insensitive to your class $C$ (here $C \in \left\{ male,female \right\}$ and the classification would report poor performance.

Alright, so, we keep this consideration of spectral resolution in mind. This actually considers the vertical dimension of the spectrogram. Let's have a look at the spectrogram itself now.

The spectrogram is the repetitive application of Short Time Discrete Fourier Transforms (STFT) to a signal in overlapping sections. Let's denote the overlap here with some factor $q \in \left\{ 0 \ldots 1\right)$. When $q=0$ there is no overlap and the spectrogram is evaluated independently on a set of windows (for example $n_{w1} \in \left\{0 \ldots 10 \right\}, n_{w2} \in \left\{11 \ldots 20 \right\}$ and so on). As $q$ approaches $1$, the STFT starts running on overlapping windows (for example $n_{w1} \in \left\{0 \ldots 10 \right\}, n_{w2} \in \left\{5 \ldots 25 \right\}$ and so on).

This overlap is desired in the STFT for a number of reasons that are related to the maths of the FT but also the signal itself and the desired analysis on it. The bottom line is that now we have to worry about one more parameter, $q$.

So what?!? So, to extend the previous example: Before, we were trying to understand what is the gender of the speaker just by examining a sample of their voice. Let's now assume that we have the recording of a conversation between males and females and we want to produce a list of time instances and their gender classification to examine (for example) if males or females tend to dominate a given discussion. Obviously, we have to tune $N_{FFT}$ from above but now we also have to tune $q$. Too much overlap means a high temporal resolution but also large amount of data (and processing). Too little overlap and we might miss some transition as some speaker was interrupted by some other speaker and then they resumed. If this phenomenon happens within the duration of a window, the STFT cannot resolve different males or females. It's a gray area that contains a mixture of both fundamental frequencies and it will produce some decision that might not be accurate.

As you can see, our discussion so far has not involved the classifier itself. The classifier is indeed some $C = f(x_m,\Theta)$ that takes the $m^{th}$ sample and maps it to one of the elements of the set of classes $C$ based on some parameter set $\Theta$. But the only thing we have been talking about so far is how do you construct this set of features that are derived from $x_m$ so that then the $f$ can do its mapping.

Let's take a look at that now: Remember how playing around with $q$ produced overlapping windows? Let's take this to the extreme. Let's think about what happens when $q=0.95$. Then, two neighbouring windows actually look very much alike. They differ only in 5% of the samples! Consequently, their spectra (adjacent vertical lines of the spectrogram) are very much alike. Therefore, it is as if we are feeding the classifier redundant data. We are showing the classifier lots of similar examples. Therefore, continuing with the above example, the classifier starts forming an opinion that MALE means all of these similar data that we are feeding into it when MALE, in our world, in the physical world, is conceptually all the different MALE voices across all recordings. A large STFT overlap could lead to overfitting!

To brake this similarity of adjacent samples, we introduce whitening. What does whitening do? It shuffles the input matrix, to reduce the covariance between adjacent samples, this is all that the whitening matrix $W$ does (but does it in a mathematical way that takes into account the covariance of the input. It produces an optimal shuffling, so optimal that any covariance between samples has been reduced to zero.)

Of course, when you bring Principal Component Analysis (PCA) into the mix, the game is changing slightly because PCA is still on the feature generation side. Applied to our previous example, PCA would enable us to get an indication of how many distinct speakers there might be in a given conversation and express the rest of the samples as coefficients to these distinct speakers (or waveforms).

So:

  As a multi-dimensional array what dimensions should an audio sample's spectrogram have before PCA whitening;

It depends on your application and specifically the spectral ($N_{FFT}$) and temporal ($q$) resolutions necessary to resolve unique examples for the classifier.

  What dimensions should the pre-preocessed object have so that it may be fed into a DBN?

It doesn't matter. You can choose to feed the whole thing and let the network "choose" the best input vectors by adjusting its weights OR you could run feature selection and reduce the dimensionality of your input (and the effort to train the network and the time it takes to train it, etc).

Hope this helps.

EDIT:
The paper states "For the application of CDBNs to audio data, we first convert time-domain signals into spectrograms.
However, the dimensionality of the spectrograms is large (e.g., 160 channels). We apply PCA whitening to the spectrograms and create lower dimensional representations. Thus, the data we feed into the CDBN consists of $n_c$ channels of one-dimensional vectors of length $n_V$ , where $n_c$ is
the number of PCA components in our representation.". Therefore, the input can be considered as "2D".


Resampling scattered data with MATLAB's $\tt interp$ and $\tt resample$

In your question, you state that your ultimate goal is to compute the discrete Fourier transform (DFT) from nonuniformly-sampled input data. It is not necessary to resample your data at all.  In fact, there are several fast algorithms, known variously as the unequally spaced fast Fourier transform (USFFT) and the nonuniform fast Fourier transform (NUFFT) that do exactly what you want: take nonuniform samples as input and produce equally spaced samples of the DFT as output. The remainder of your question, dealing with issues related to MATLAB's resampling routines, are a red herring.

Here are some publications describing the USFFT algorithms:

G. Beylkin, On the Fast Fourier Transform of Functions With
Singularities,  Applied and Computational Harmonic Analysis, 2, pp.
363-381, 1995. Link to preprint.
A. Dutt and V. Rokhlin, Fast Fourier Transforms for Nonequispaced Data, SIAM J. Sci. Comput., 14, pp. 1368, 1993.

There are many commercial and free implementations of these algorithms. You can access some of them here.

EDIT: Here is a "quick and dirty" approach that may help you.

Based on your problem description and comments, you are just trying to get a sense of the spectral content of your signal, which does not require high accuracy. Here is one way to formulate your problem so that the spectrum can be easily estimated using any of the USFFT implementations in the above link.

You have a signal $f(t)$ defined on the $T$-duration interval $t\in[-T/2,T/2]$, and you wish to estimate its Fourier transform
$$\hat{f}(\omega)=\int_{-T/2}^{T/2}f(t)e^{-i 2 \pi \omega t}\,dt$$
at a number of values of $\omega$. You have collected $N$ samples of your function
$$f_n=f(t_n),\quad n=1,\ldots,N$$
at the unequally spaced time locations $t_n$, $n=1,\ldots,N$.

In general, you would solve this by formulating a least squares problem for the Fourier transform $\hat{f}$, and using the USFFT to rapidly solve the system of equations, but since you do not require high accuracy, here is a "quick and dirty" solution.

Use the trapezoid rule from Calculus 101 to discretize the above integral, giving the approximation
$$\hat{f}(\omega) \approx \frac{\Delta_1}{2}f_1 e^{-i2\pi\omega t_1}+\left(\sum_{n=2}^{N-1}\frac{\Delta_{n-1}+\Delta_{n}}{2}f_n e^{-i2\pi\omega t_n}\right)+\frac{\Delta_{N-1}}{2}f_N e^{-i2\pi\omega t_N}$$ where
$$\Delta_{n}=t_{n+1}-t_n$$
is the temporal spacing between samples $n+1$ and $n$.

Define the new sequence
$$g_n=\begin{cases}\frac{\Delta_1}{2}f_1, & n=1 \\ \frac{\Delta_{n-1}+\Delta_{n}}{2}f_n, & n=2,\ldots,N-1 \\ \frac{\Delta_{N-1}}{2}f_N, & n=N\end{cases}$$
Now pick some frequency grid you care about. For example, say you want $2000$ equally spaced frequencies between $0$ Hz and $25$ Hz, and define
$$\omega_m = \frac{m-1}{2000-1}25,\quad m=1,\ldots,2000.$$

Now your problem amounts to evaluating the sum
$$\hat{f}(\omega_m) = \sum_{n=1}^N g_n e^{-i2\pi\omega_m t_n},\quad m=1,\ldots,M.$$

This looks like an FFT, except the $t$-values are unequally spaced and the $\omega$-values are unusual. If you do not care about speed, then just evaluate the sum directly using your favorite programming language.  It will take a minute or two of computer time, but so what? If you are just doing this a few times to get a sense of what your signal's spectrum looks like, then it does not matter how fast it is. If at some point in the future you decide that you do care about speed, then the sum can be evaluated extremely rapidly using a single call to a USFFT library, which you can download from the link I posted (or you can use one of the many commercial implementations).  Making the switch to the USFFT will only require you to modify a handful of lines of code.


Detecting blobs that connects to any other blob, maybe with OpenCV

Flood fill  the background, then hough transform to pick out rough circles II and III


Filtering artefacts and filtering short signals

All real-time filtering (as opposed to post processing) wirh FIR and IIR filters will have start up transitions based on the state of the filter at start up. For optimum rejection of AC noise , Instead of a Butterworth Filter consider using an 2nd order IIR notch filter with the notch set at your AC frequency (such as 50 Hz). A design for this is further explained on this site and is simple and effective for this purpose, and happens to be demonstrated with a notch at 50 Hz. (Transfer function of second order notch filter).
To have shorter segments, partition into segments at the output of the filter so that the filter’s memory can be maximized to meet the filtering requirement as designed (a filter’s transition bandwidth between passing a frequency of interest and rejection is dependent on the span of its memory—- for tight rejection long memory is required.)
If the application is for very short durations at completely different instances in time, that can make no use of the prior immediate time domain signal, then such traditional FIR and IIR filtering techniques will not be suitable beyond what they are capable of for the given time span. That said, the IIR nulling approach with its rejection optimized at the interference frequency may be sufficient and simple.


MFCC extraction for speaker verification

You get a vector of 39 coefficients for each short frame of the signal; not for the entire signal. A speech signal contains segments of voice and silence, and of course when voice is present, there are different phonemes spoken in sequence - so the signal is not stationary and the MFCCs will fluctuate. Extracting the coefficients on a single frame spanning the entire signal would blur all this.

You might also be tempted to compute the average of the MFCC vectors. There are indeed a handful of applications in which this might work (for example categorizing sounds, or very coarse music genre classification), but I don't advise you to do this for speaker verification.

In speaker verification or recognition applications, you perform the classification or model scoring for each frame. For example, let us say you have trained a GMM for a speaker, and then are given another signal and want to know if this was spoken by the same person. Extracting the MFCC on the signal will give you a sequence of 100 MFCC vectors. You score each vector with the GMM and make a yes/no decision for each vector. Finally, you use voting to aggregate the decisions over the 100 frames. Or alternatively, you average the GMM scores for the 100 frames and use the result to make the yes/no classification.

When people report that they have used 39 coefficients, they mean "per frame".


How do I get the Time information after STFT?

How, can I find the change of frequency with Time (900000 sample) or 15 min?

OK, so you need 900000 frequencies.  Well, one way to do that is to take each element of your 637 length vector and expand it to (up to) $H$ samples.

  Moreover, how does changing the window length and overlap affect the outcome of STFT?

The longer the window length, the better the frequency resolution is. The shorter the length, the worse the frequency resolution will be.

Before we can answer your question, we will need to know what form your signal takes.

If it's a single tone in noise, then you may be able to use one of these or one of these or this approach.


Let me give you a quick answer, focused on ways you can improve and research for a deeper knowledge by yourself.

1 - Kolmogorov and Shannon  Entropies are related to measures of information and strongly related to one another. (See the wikipedia page for Kolmogorov complexity for further info.)

2 - Boltzmann entropy, which is an older concept, is a completely differenty and unrelated concept from the theory of Thermodynamics. (Again, wikiepdia...)

3 - The formal definition of kolomogorov entropy uses the Shannon entropy in the formula.

If you are interested in the information entropy concept, i suggest you read the Book by Cover &Thomas, which has an entire chapter on Kolmogorov Complexity.


Is this HW question even valid? (linear systems)

there is an SE specifically for signal processing.  sorta surprized to see tags for this at the Physics SE.

what you are suspecting is precisely true.  LTI is a property of the system, not of the input.  both 1. and 2. are obviously linear and time-invariant.

one slight misconception.  in DSP we usually use square brackets, as in "$x[n]$", to depict discrete-time signals (as an alternative to "$x_n$").  a discrete-time signal is just a sequence, but we want to leave the subscript for the same purpose we might with continuous-time signals.  i.e. given a vector of continuous-time signals: $\{x_1(t), x_2(t), ... x_N(t) \}$, if we were to sample those signals, the corresponding discrete-time vector would be $\{x_1[n], x_2[n], ... x_N[n] \}$.  the argument (inside of square brackets) to a discrete-time signal is always an integer.

that said, "$\delta[n]$" is the Kronecker delta, not the Dirac delta.  (the latter the math guys like to tell us is not even a "function".)

$$ \delta[n] \ \triangleq
\begin{cases}
1, & \text{if }n = 0 \\
0, & \text{if }n \ne 0
\end{cases} $$


Can an infinite bandwidth signal be an infinite-time signal?

Let's consider the continuous time case.

Take a Gaussian $f(t) = \frac{1}{\sqrt{2 \pi}} e^{-t^2/2}$, which is positive for all time.

The Fourier transform of a Gaussian is another Gaussian: $F(\omega)= c e^{- d \omega^2 }$ for some constants $c,d$  that you can work out on your own, so a Gaussian occupies all time and all frequency. That is, you can be both infinite in time and infinite in frequency.

Your trade-off between resolution in time and frequency is given by the uncertainty principle. The Gaussian attains equality in the uncertainty principle.

It seems like you're somewhat confused as to what Nyquist's sampling theorem says -- see this page for details (or your favorite signal processing text). You're assuming the signal is band-limited (i.e. is non-zero in some finite interval of frequencies) in order to be able to sample the signal at some prescribed rate in order to be able to perfectly reconstruct the signal.


Converting transfer function that is a sum of unusual rational polynomials to finite difference equation


In block codes (channel coding), in practice how is the message length ($k$) fixed?

There are a few rules of thumb:

For a fixed code rate $r = \frac{k}{n}$, the amount of gain you observe from applying coding typically increases as you increase the block size $n$. Accordingly, the complexity of decoding the larger blocks is often larger.
If you decrease the number of information symbols in the block $k$ and keep $n$ fixed, the code rate decreases. This will result in higher coding gain, but lower information throughput for a given channel symbol rate.
If you increase the number of information symbols in the block $k$ and keep $n$ fixed, the code rate increases. This will result in lower coding gain, but higher information throughput for a given channel symbol rate.


I am not really sure what is meant by "fake face" but for the purposes of this response, let's consider that you have two classes "REAL" and "FAKE" and representative example data points from each class are positioned in some $n$ dimensional space, with $n$ denoting the number of features.

In the described situation, the Minimum Distance Classifier (MDC) is used to discover the decision boundary between the two classes (the MDC can work with more than two classes if required but let's just constrain it to two here).

The MDC discovers the decision boundary by first estimating the centroid of each class by averaging the feature values. Once this is done, the classifier is essentially "trained". The knowledge of what it means to belong to one or the other class is basically depending on the position of the centroids in multidimensional space.

The next step is to find the line that passes through the two centroids and its vertical that is further at the half-way point between the two centroids. This vertical line at the half-way point is the actual decision boundary. Anything landing to the right of this line is classified as belonging to one class and anything landing on the left of this line (or on the line itself) is classified as belonging to the the other class.

For the purposes of illustration, I will be re-using one of the figures from a previous link here, let's call the "O"s -> REAL and the "X"s -> FAKE:

So, here is a simple "exercise". Forget about the vertical line being at the HALF-WAY point between the two classes and let's try and draw it a bit to the LEFT, that is, closer to the REAL class.

Obviously, this will classify some of the examples of Class REAL as Class REAL (True Positives), some of the examples of Class REAL as Class FAKE (False Negatives), some of the examples of Class FAKE as Class REAL (False Positives) and some of the examples of class FAKE as Class FAKE (True Negatives).

This "exercise" just gave us the 4 numbers we need to calculate the two metrics of False Positive Rate (FPR), True Positive Rate (TPR) that we need to derive one point of the ROC curve at some threshold value.

So, in the case of the MD classifier, your threshold is equivalent to the position of the decision boundary between the two classes. By "sweeping" that decision boundary all the way from the centroid of one class to the centroid of the other class and deriving FPR, TPR for each value, you will be able to derive the complete ROC curve.

Another way to "visualise" this process with respect to the textbook examples of 'sweeping the threshold value', is to imagine that you rotate your feature space clockwise, until the decision boundary is vertical to the x-axis. If you were to do this, you would be practically eliminating one of the features and you would be "back" to the one dimensional case that is depicted here (for example), with each "gaussian" curve representing the density of the points above that x-axis value. (This is suggested here just as an aid, in case it helps you to make the connection between the position of the decision boundary and the threshold. You don't have to actually rotate the feature space).

Hope this helps.


What data should I use to test an FFT implementation, and what accuracy should I expect?

If you want to verify an FFT algorithm for correctness, in the sense that it performs the desired function that has the known properties of the discrete Fourier transform, then you can use the approach proposed in:

  Ergün, Funda. (1995, June). Testing multivariate linear functions: Overcoming the generator bottleneck. In Proc. Twenty-Seventh Ann. ACM Symp. Theory of Computing. (p. 407–416).

The above paper is referenced by the makers of FFTW as their method of choice for verifying that a particular FFT implementation does what it should. The proposed technique distills the function into three main components that are verified with separate tests:

Linearity: The DFT (along with its other cousin transforms in the Fourier family) is a linear operator, so for all values of $a_1, a_2, x_1[n], x_2[n]$, the following equation must hold:

$$
\mathrm{FFT}(a_1 x_1[n] + a_2 x_2[n]) = a_1\mathrm{FFT}(x_1[n]) + a_2\mathrm{FFT}(x_2[n])
$$

DFT of the unit impulse: A time-domain signal equal to the Kronecker delta function is applied to the input of the FFT algorithm and the output is checked against the known DFT of the unit impulse function (it transforms to a constant value in all output bins). If the FFT algorithm provides an IFFT, it can be tested in reverse to show that it yields the unit impulse function again.
Time shift: Two sets of data are applied to the input of the FFT algorithm; the only difference between the two in the time domain is a constant time shift. Based on the known properties of the DFT, this should effect a known linear phase shift between the two signals' frequency domain representations, where the slope of the phase shift is proportional to the time shift.

The authors of the paper assert that these tests are sufficient to validate the correctness of an FFT implementation. I've not used this technique in the past, but it does seem to make sense, and I would trust FFTW's authors (who have produced a great piece of free software) as credible authorities on good approaches to the validation problem.


You switched the $a$ and $e$ in your third line of code. If you write it like this, then it will be behave as expected:

H = 1.0./(1.0-a*exp(-j*W));


This is how my professor is finding the frequency response of an LTI system when given the impulse response. Is this wrong?

Your professor is right, and you're almost right too. The filter is clearly an FIR filter, but because its frequency response can be expressed as a geometric series, a recursive implementation is possible.
If you write the transfer function as a rational function you get
$$H(z)=2\frac{1-z^{-12}}{1+z^{-2}}\tag{1}$$
which is almost the same as you got, apart from a constant factor.
Note that there is a pole-zero cancellation in $(1)$, so the filter is FIR with impulse response
h = [2  0 -2  0  2  0 -2  0  2  0 -2];
So apart from rounding errors both

b = 2*[1,zeros(1,11),-1];
a = [1,0,1];
[H1,w] = freqz(b,a,1024);

and

[H2,w] = freqz(h,1,1024);

must give the same result, which turns out to be true. On my machine I get

max(abs(H1-H2))
ans =    7.9623e-14

Note that the phase is ambiguous, so the phase values are only identical up to an additive multiple of $2\pi$.


Follow up question regarding: "Complex sampling can break Nyquist?"

To represent a complex baseband signal in the range $[-B/2,B/2]$ you need either $2B$ real samples or $B$ complex samples per second.

To take real samples, you write the signal as $I(t)+jQ(t)$, with $I(t)$ and $Q(t)$ real, and sample each one at $B$ real samples per second, for a total of $2B$ real samples per second. It should be clear that you can't downsample without aliasing.

You can also take complex samples at $B$ samples per second (which is just a different way of thinking about the real sampling described above). This allows you to represent the entire band $[-B/2,B/2]$. Downsampling in any way will also result in aliasing.

To be more specific regarding your application: You have a bandpass signal with bandwidth $B$ and take real samples at rate $2B$. When you calculate the complex envelope, the I and Q components can be sampled at a rate of $B$ real samples per second each. The complex envelope can also be sampled at a rate of $B$ complex samples per second. In any case, you end up with a total of $2B$ real samples per second.


It should be clear from the definition of the STFT that the result is invariant under global phase changes. With the restriction to real input signals that still provides for a sign ambiguity, or a phase difference of $\pi$. That means, even if your reconstruction is perfect, the result may differ by that phase from the original.

Another reason for the failure to converge may be the existence of nearly disconnected components in your signal. Say you have two sound components that are temporally separated so that there is no STFT frame that overlaps with both. Each component may be reconstructed with a different sign, giving you four different reconstructions with identical STFT. The same works for components that are separated well enough in frequency direction. Each component doubles the ambiguity of the reconstructed signal and will increase the probability of a deviation from the original signal.

The final reason I can see lies in the nature of the iterative algorithm. Fixed point iterations, of which kind this algorithm is, can fail to converge to a proper fixed point, keep oscillating or even diverge. Unless the maps involved are well behaved in a certain sense, predicting convergence is hard. It seems that your iteration doesn't settle, and it's unclear if adding iterations will help. So this may apply here, together with the other two problems.


Yes, you can always use pixel units for $\sigma$. In 2-D you will have $\sigma_X$ and $\sigma_Y$.


Does only filter coefficients can represent it's own filter characteristic?

You can compute the filter's frequency response to check the behavior of the filter:

$$H(\omega)=\sum_{n=0}^{N-1}h[n]e^{-jn\omega}\tag{1}$$

where $h[n]$ are the filter coefficients, and $N$ is the filter length (i.e., the number of coefficients). The frequency response $H(\omega)$ is a complex function. You can compute its magnitude and phase to evaluate the filter's magnitude and phase characteristics.

You can compute (and plot) the frequency response using Matlab's freqz command.


Yes your computation is correct. Indeed, you could have seen this even easier, had you considered that the output computation sum was

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

And this is clearly an LTI (linear time-invariant), FIR (finite impulse response) system with associated impulse response

$$h[n] = \delta[n-1] + \delta[n] + \delta[n+1] + \delta[n+2] $$


Smoothing 3D Data for the Second Derivative

I have always found that a "Median Filter", https://en.wikipedia.org/wiki/Median_filter   to be effective on streaming data.
A brief intuitive description is below; but you must take steps to insure you restrict the bandwidth to the real data source you expect.  While you might not know the sources/spectrum/timing of noise and spikes (in which case the median filter is actually pretty good) you hopefully know the signals and characteristics of the signals you are looking for.  For extensive analysis, still iffy IMHO, along the lines of PDF (probability distribution function) and spectral analysis: there are some IEEE articles; but my analysis of most of the articles were that they weren't really useful in the face of unknown interference.  Aside from accommodating only the real signal bandwidth I wouldn't try z-s domain filtering without a really good reason.

------------------------ Median filter abbreviated description
Basically you pass the data into a shift register and pick the middle value.  Examples always show a sort process; but since you lose only one value and gain only one value during a shift, you can bypass doing full sorts each step
Having said that: if the spikes are in the input data then that is where this should be applied; but if the spikes are real "jerks" (i.e. occurring in the second derivative) then that is where it should be applied.
Since the derivative process amplifies higher order noises one has to careful.  If I was doing what you describe I would impose the median filter on the source data, first differences, and then on second differences (acceleration).  And then examine the results; look inside the algorithm and understand what data is tossed!
There is an alternate form: you examine a shift register of the data for min's and max's and just toss the extremes.
Both techniques incur some data loss that has to be taken into account.  You can interpolate, which is a little expensive, or just accommodate a lower sampling rate; say one in ten samples will always be tossed.


You can try to improve one both sides. One issue is that smooth separable kernels $w$ or integration scheme like trapz tend to reduce rank, and pose conditioning problems. Did you use pinv or X = (A)\(b) already?

You can invest on non-separable or higher order integration rules, see for instance n-dimensional simplex quadrature;
You can invest in regularization, for instance with the simple diagonal update: $B_c = B + c*\operatorname{eye}(\operatorname{size}(B))$. If this is not sufficient, there are quantities of constrained or penalized methods (LASSO for instance), implemented for instance in the CVX toolbox (e.g. Premade solvers for specific problems (vector variables)).


Your understanding is correct. There is no difference between a phase of $\pi$ and a phase of $-\pi$. You can always add or subtract integer multiples of $2\pi$ to the phase without changing anything because $e^{j2\pi n}=1$ for $n\in\mathbb{Z}$. Clearly you have $e^{j\pi}=e^{-j\pi}=-1$.

The probable reason why they used opposite signs in the phase plot for positive and negative frequencies is because then you can directly see from the figure that the phase (of a real-valued signal) is an  odd function, but that's a purely cosmetic matter.


You need the right audio features, and the modern way to do this is with deep neural networks, which learn them automatically. I would use one built for speech recognition (it's basically the same problem), then model the users by the output of the top layers. Accents should manifest as clusters of the features, and languages as superclusters. I don't have much experience here but I'm making an educated guess.

Here is a relevant paper (Deep Neural Networks for Acoustic Modeling in Speech Recognition), and here is a relevant student project (Accent Recognition with Neural Network).


The resulting up chirp traverses across BW and arrives at the same place it started.

not quite exactly right, but very close: it ends up in the DFT bin before the bin it started from. This makes the chirps (frequency domain) cyclically shifted versions of the straight-up "prototype chirp". (my name invention)
Maybe this way of looking at it helps intuition:
Let's assume the system is synchronized in time¹ and let's look what happens when you multiply the unshifted downchirp with the upchirp.
You're multiplying a complex sinusoid with a frequency slope that is the exact inverse of the slope of the other sinusoid. In other words, you're multiplying two sinusoids which only differ in the sign of their exponents. Now, multiplying two exponential functions with the same base ($e^\cdot$) will lead to a new exponential function with the sum of the exponents of the factors as exponent. Yay! That means up- and down-chirp cancel out, you get a constant $e^{j0}=1$.
Now, if the chirp is shifted, the sum of the exponents doesn't cancel - but it becomes a linear function – as argument of the complex exponential: a tone!
A tone that on top of that falls in the raster of the DFT. So, a single peak after FFT.

As per DSP principles, the multiplication of received and locally generated base band chirps is the addition of the two instantaneous frequencies of the chirps in red and purple

yep.

Then, I have drawn a line which represents their summation in yellow.

And here you simply forgot to realize that your positive frequency is outside Nyquist bandwidth. You need to subtract the sampling rate, and you'll see that it aliases to the same frequency as the rest of the yellow line.
Don't forget that your in discrete time and frequency domain: the spectrum must be imagined to repeat infinitely every multiple of the sampling rate.

¹ Unsolicited words of advice: if you find literature that proposes cool signalling schemes and uses the sentence Assume the system is synchronized, don't invest too much time into that scheme until you find information on how that can be synchronized... synchronization is one of the hard parts, and it has broken many a system's promises.


Study materials for Decision feedback equalizer

Not a comprehensive source of details, but for me this paper/tutorial has been the most helpful in terms of principles explanation. Also here you can find a basic working example of DFE implementation in python, that you may find helpful.


Does a filter need to completely attenuate high frequencies to be considered low-pass?

Ideal low pass filter is supposed to attenuate/remove all the frequencies above the cutoff. but in practice the filters have non zero magnitude even for frequencies in their attenuation region, still they are called low pass.

if the magnitude of filter is quite high for frequencies above cutoff it might be called bad lowpass filter, but is still a low pass because it attenuates the frequencies above cutoff.


Minimum FFT sample time - guitar frequencies and harmonics

Your highest frequency is 831 Hz; since undersampling to ignore the sub-82 Hz would only save you 10% sampling rate, let's ignore that possibility.

So, Nyquist says your sampling rate should be at least 1.662 kHz; now, your problem is, as @fibonatic mentioned, that a the freq. resolution of any DFT is alwas

$$\Delta f = \frac{f_\text{Nyquist}}{N_\text{DFT}}\\
\implies\\
N_\text{DFT} = \frac{f_\text{Nyquist}}{\Delta f}$$

and the time taken by capturing of samples you'll need for that DFT is

$$T = \frac{N_\text{DFT}}{f_\text{sample}} = \frac{N_\text{DFT}}{2f_\text{Nyquist}} = \frac{\frac{f_\text{Nyquist}}{\Delta f}}{2f_\text{Nyquist}}= \frac{1}{2\Delta f}$$

So no matter at which rate you'd sample, there's nothing oyu can do about the time period needed for a full-resolution DFT.

You could just pad a shorter observation with zeros, and that will increase your resolution. Point is that you trade resolution for quality, and the variance is going to grow significantly.

The DFT is simply not the right choice here; if you need 2Hz resolution at the lower frequency, you get the same at the higher frequencies without any need.

Honestly, think about it: you're planning to use a DFT that gives you 4096 frequency bins, where you'd need to distinguish maybe 100 tones, if that many.

I'm pretty certain that the cheap guitar tuners you get don't work that way. I'd assume they simply employ dedicated hardware (not firmware like you'd need to write) to implement a simple filterbank – possibly even a IIR filterbank, and if I had a guess, at very low sample bit width, with simple biquads. Or, even simpler, they just have couple of analog filter bandpasses, and just undersample the resulting sample streams and do the same digital filters on all the octaves they preselect that way.

If I had to do this, I'd atually do the same on your microcontroller (by the way, 8192-FFT sounds gigantic for a Cortex-M4, but you know your device better than I do), but with different architecture. Do a rough calculation: If you keep the sampling rate low, how many FIR filters of the different lengths you'll need can you let run over the data? Does implementing your filterbank as a PFB filterbank help? Does making a hierarchic decision (e.g. just first an "octave selector", then depending on the result of that, throwing away half, $\frac38$,$\frac78$,… of the samples and then doing a tone-detection filterbank.


Best way to do accurate amplitude measurements using FFT

You can interpolate using the transform of your chosen window function as the interpolation kernel (similar to Sinc reconstruction), and then use successive approximation to find the peak of the interpolated spectrum.
Pre-constructing a polyphase table using an extremely oversampled transform of your chosen window might or might not assist in this computation.


There's a good overview article which appeared in 1996 in the IEEE Signal Processing Magazine: Splitting the unit delay: tools for fractional delay filter design. The nice thing about it is that there's also a set of related Matlab files available. These routines will allow you to design such a system.

As for shifting by an integer number of samples, what you did is obviously correct.


Are you looking for a specific function to perform this task? Otherwise this could probably be done with a simple for loop.

Determine your unique rows using the function

  bins = unique(A, 'rows');

Then just loop through each row of A and count each instance of each bin you find in A.

This is probably good enough if you dont care about how long the program runs (which still isnt long at all if your A matrix isn't ridiculously large), but keep in mind MATLAB doesn't "like" running loops.


When you look on a signal's DFT transform what you see is given in normalized units in the range $ [ -\pi, \pi ] $ where $ \pi $ equals the sampling rate.

Since the DFT is periodic and the way FFT work when one calculate the DFT in modern software (MATLAB / Python, etc...) the DFT is given in the $ [0, 2 \pi ] $ range.

In order to view it in the more intuitive way one could use fftshift() (MATLAB).

When you display it then, probably you won't see it gets to zero anywhere.
What you'll probably see (If the signal is band limited with bandwidth lower than the sampling rate) is a big drop in the magnitude.


How do I find a system's impulse response from its state-space repersentation using the state transition matrix?

You can approach the problem using the state transition matrix by solving the standard non-homogeneous ODE in the first equation. The solution to $\dot{x}(t)=A x(t) + B u(t)$ is

$$x(t)=x_0 e^{At}+\int_{0}^te^{A(t-t')}Bu(t')dt'$$

where $x_0=x(0)$. The quantity $e^{At}$ is called the state transition matrix (also the solution to the homogeneous ODE), which I'll refer to as $\Xi(t)$ (I don't recall the standard notation for this). Taking $x_0=0$, the equation for $y(t)$ becomes

$$y(t)=C\int_0^t\Xi(t-t')Bu(t')dt'+Du(t)$$

The above equation gives you the output as the input convolved with the system impulse response and indeed, you can take the Laplace transform of the above equation to verify. Noting that the Laplace transform of $\Xi(t)=e^{At}$ is $(sI-A)^{-1}$ and convolutions in the time domain become products in the s-domain, we get

$$Y=C(sI-A)^{-1}BU+DU$$

which gives you the same transfer function as in your question.

Regarding your comment on the fully Laplace transform approach being long, I wouldn't necessarily say it is so. However, state transition matrix approach might be simpler to implement, because several operations involving it can be computed with simple matrix multiplications and nothing more.


How to use 'randi' in MATLAB to generate equally probable samples which are either 0 or 1?

Generate random integer uniformly distributed between 1 or 2. Then you can shift it to your liking.
X=randi(2,size,1)-1


Is this a robust method for analysing correctness of a filtered signal? Is there a more robust method?

Given you can have a "reference waveform" representing the characteristics of your expected signals, to test the "correctness" of the signal extracted from noise use the Normalized Correlation Coefficient. I reference other posts below where I give the details of computing and using the Normalized Correlation Coefficient:

Noise detection

How can I find SNR, PEAQ, and ODG values by comparing two audios?

This will be an interesting approach to compare your different algorithms as the result of removing portions of your signal of interest or keeping noise contributions will both reduce the correlation coefficient as expected; so this will be an excellent approach to accurately indicate which algorithm you use is optimum for your waveform.


It is not clear if you want to visualize the signal or process/transform it ("graphic equalizer" is an effect/transformation adjusting the signal level independently in various frequency bands, not a visualization - a more correct term for the representation you posted an image of is "spectrum analyzer").

It is not clear if you are working in the analog domain or the digital domain, and what kind of technology you want to employ (analog electronics? computer code running on a DSP chip or microcontroller? on a desktop computer?).

Are your signals audio signals? In which case you have to be aware that a meaningful representation should display the 20 Hz .. 20kHz range (range of human hearing), and that it should be subdivided logarithmically.

To get this kind of visualization, there are two solutions:

A bunch of band-pass filters (and a low-pass + a high-pass for the lower and higher channels). You then measure the (log-) energy in each band to get the "height" of each bar. This can be done in the analog domain (a bunch of op-amps for the filter, a rectifier to measure the signal level, and a log-amplifier to get the energy on a log scale; or with dedicated ICs like the MSGEQ7); or in the digital domain with digital filters.
A FFT + some grouping of adjacent bins (which is more or less equivalent to computing the filter bank in the frequency domain). Note that this can only be used for digital signals.


Calculating PSD, Cross PSD by using FFT outputs

Cross PSD (fourier transform of cross-correlation) is $X(f)^*Y(f)$, so it will $(a-ib)(c+id)=(ac+bd)+i(ad-bc)$. I think you missed to take conjugate like you previous question (Calculating Cross Power Spectral density between two complex signals).

Cross PSD tells a lot about correlation between two different signals. If Cross PSD is white (meaning near flat at all frequencies) with very low value implies both signals are uncorrelated. If $x$ and $y$ are input and output of a system, you can measure the coherence between $x$ and $y$. The relationship between PSD of $x$ and $y$ for the LTI system $h$ is $G_{yy}(f)=|H(f)|^2G_{xx}(f)$. You can read more at https://en.wikipedia.org/wiki/Coherence_(signal_processing)


How does mic combined with loudspeaker work?

By using Acoustic Echo Cancellation. Since the signal that goes to the speaker is known you can subtract it out at the microphone, provided you know the exact transfer function from speaker to microphone. The determination of the transfer function is typically done through an adaptive filter. Similar techniques can be used to reduce reflections, reverberation and echoes from the room.

In addition, there are often multiple microphones using beam forming. The beams can be steered so that the pickup from the speaker is minimized (basically steering a null in the polar pattern to the speaker) and it can also be used to maximize pickup of the desired sound.

All of these things are really difficult, complicated and hard to get right. Cancellation and adaption often create artifacts that significantly reduce speech intelligibility. The quality of commercial implementation varies all over the place. In my name brand laptop, it's so bad that it's worse than useless and I simply turned it off. I have tried multiple USB speakerphones and also found a huge variability in quality. Finally I found one that actually works :-) but product recommendation are out of scope of this site. Let's just say that the one pictured is fairly decent.


A frequency bin 'only' covers the value of its center frequency (+ potential spectral leaks to that center frequency). Is that correct?Is that correct?

No. All frequencies show up in all frequency bins EXCEPT frequencies that are an integer multiple of $f_{\Delta}=f_s/N_{FFT}$ sample rate divided by FFT length.

In other words a frequency bin is affected by ALL frequencies except the center frequencies of the other bins.

You need to take into account that the FFT assumes that the signal is periodic with N. For a periodic signal contains only integer multiples of the fundamental frequency.


Machine Vision Techniques for Sign Language Recognition and Interpretation(using Matlab)

I am afraid that the question is trying to probe a very wide range of issues to receive a truly meaningful answer.

In general, this would be a classification problem. As such, you would need to decide on a number of features which would reflect the key changing characteristics of your images and a classifier which, after being trained, would be able to discriminate between your different gestures. For more information please see: http://en.wikipedia.org/wiki/Statistical_classification

The issues of how many training images you need per gesture and which features and classifier, all depend on your environment and specifics of your application.

The simplest thing you can do is template matching, where you would simply compare an acquired image of a gesture with stored images of gestures. For more information please see: http://en.wikipedia.org/wiki/Template_matching

But of course, this approach is a tiny step forward. What about changing backgrounds, different pairs of hands, scaling of the objects because of different distances to the camera, slight rotations, occlusions... ? And then, once the static recognition is done, you will have to consider the version that works on uninterrupted video.

There are things you can do to deal with these issues but for robust recognition, it is likely that you will need to move to 3D data representing the movement of hands in space. (For example: http://www.sciencedirect.com/science/article/pii/S1077314210001748 ).

You can acquire 3D data in a number of ways, e.g. virtual reality gloves (please see: http://dx.doi.org/10.1109/IJCNN.1999.832699), devices such as kinect (please see: http://dx.doi.org/10.1109/ICSESS.2012.6269439) and plain simple accelerometers (please see: http://dx.doi.org/10.1109/JSEN.2011.2166953).

The simplest method to start with 3D data is the one based on accelerometers. In this case you would be working with time-series data (e.g. timestamp, acceleration vector).

Happy to amend the response if you provide more specific data about your application.


Phase difference measurement of a 100 kHz tone

The achievable rms accuracy is approximately the reciprocal of the SNR as an rms quantity, specifically  $10^{-35/20}$ in your case assuming the small angle criteria, if your noise is all phase noise. It is likely your SNR is equally phase and amplitude noise, so the phase result as limited by your noise would be 3 dB better after hard limiting the signal to eliminate the AM component.

I linked other posts that include multiple variants for phase detector implementations, depending on if you have a complex or real representation of the signal. The simple approach is to hard limit the signals (as most phase detectors are sensitive to AM) and then use a multiplier. The output of a linear multiplier is proportional to the cosine of the phase angle  between the two inputs.

Phase Detectors: Phase synchronization in BPSK

Phase Detection for complex signals: How to correct the phase offset for QPSK I-Q data

Phase Detection for Carrier Recovery: High modulation index PSK - carrier recovery


Removing known sequence from an audio sample

From the Description, I think Adaptive filters will work out for this application.
As you have the reference signal("I have an audio recording of that same movie playing on that same television with no conversation in the background") and the signal from which you need to remove this reference.


Relationship between sampling frequency and number of coefficients

For digital filters, all frequencies (and bandwidths) are relative to the sampling frequency. So a fixed bandwidth (in Hz) becomes relatively more narrow for a higher sampling frequency. You would get the same number of coefficients (for the same performance) if the ratio $B/f_s$ were constant, where $B$ is the filter's bandwidth and $f_s$ is the sampling frequency.

I think that is basically what you already suspected, even though I think that both reasons you gave are actually equivalent.


How to implement an improper discrete-time transfer function into Simulink

I think you missed the point in the link you posted. Implement the improper function as a differentiator. You cannot get rid of the proper part but only replace it with differentiators. Hope this helps!


What is the difference between phase noise and frequency noise?

Phase Noise and Frequency Noise are not two different noise sources, they are artifacts of the same noise, it is just a matter of what units you want to use.  Frequency and Phase are directly related as frequency is phase changing with time, so if you have one you will always have the other; frequency and phase are related by derivatives and integrals: the derivative of phase equals frequency  (as a change in phase versus time is frequency).

For a waveform in the time domain that represents the frequency changing with time as $f(t)$, the equivalent function as phase changing with time would be:

$$\phi(t) = \int f(t) dt$$

And likewise to go from phase to frequency:

$$f(t)= \frac{d\phi(t)}{dt}$$

Note this can be very confusing as we are talking about a waveform whose value represents frequency changing with time; so to be very clear this waveform I am describing is indeed in the time domain (as the independent variable is time), and we are observing how the frequency changes versus time. We can take the Fourier Transform of this waveform, and then observe the rate of change of frequency (which is the frequency of our frequency variable....even more confusing!).

However getting that straight is important to understanding the difference between Phase Noise and Frequency Noise:

Integration in the time domain is a low pass filtering response in the frequency domain, as indicated by the Laplace Transform property for integration:

$$L\{ \int x(t)dt\}=\frac{1}{s}X(s)$$

And similarly the derivative in time results in a multiplication by s in Laplace, which represents a high pass function (as you can get the frequency response (or Fourier Transform) for both by setting $s=j2\pi f$)

Therefore a Phase power spectral density and Frequency power spectral density are related as follows:

$$S_f(f) = (2\pi f)^2 S_{\phi}(f) $$

Where $S_f(f)$ is the power spectral density due to frequency fluctuations, and $S_\phi(f)$ is the power spectral density due to phase fluctuations.

For example, if you had a white frequency noise process (meaning the FREQUENCY power spectral density was flat for all frequencies), the PHASE power spectral density would be going down at the rate of $1/f^2$ or 20 dB/decade.

Here is a plot I have showing the relationship between Phase Noise and Frequency Noise as time domain plots. Notice how the phase is smoother but wanders off relatively slowly with varying random offsets (low pass with random walk), while the frequency plot has higher frequency content but does not deviate from 0 for long durations (high pass filtered).


What is correlation kernel and compare with gaussian kernel

Referencing the paper, the convolution kernel can be written as:

$d(|y-x|)$ is the Euclidean distance between the centre pixel $y$ and a location in the kernel $x$.  $\alpha$ is used to ensure that the entire area under the kernel is 1.  However, you did not specify how big this kernel is.  As such, we will specify the rows and columns of this kernel to be M and N respectively.  Let's also assume that the size of the kernel for each dimension is odd.  The reason why is because the shape of the kernel will be an even square and makes implementation easier.  As such, here are the steps that I would perform to do this:

Define a grid of X and Y co-ordinates, and ensure that the centre pixel is at 0.
Compute each term in the convolution kernel without the $\alpha$ term.
Sum up all of the terms in the kernel, then divide every value in this kernel by this term so that the entire area of the kernel is 1.

Let's do this step by step:

Step #1

We can do this by using meshgrid.  meshgrid (in this case) creates a 2D grid of (X,Y) co-ordinates.  X defines the horizontal co-ordinate for each location in X, while Y defines this vertically. By calling meshgrid(1:m, 1:n), I am creating a n x m grid for both X and Y, where each row of X progresses from 1 to m, while each column of Y progresses from 1 to n.  Therefore, these will both be n x m matrices.  Calling the above with m = 4 and n = 4 computes:

m = 4;
n = 4;
[X,Y] = meshgrid(1:m, 1:n)

X =

 1     2     3     4
 1     2     3     4
 1     2     3     4
 1     2     3     4

Y =

 1     1     1     1
 2     2     2     2
 3     3     3     3
 4     4     4     4

As such, we simply have to modify this, but we ensure that the middle is at (0,0), and also ensure that the size of X and Y are odd.  Let's say that M = 5 and N = 5.  We can then define our X and Y co-ordinates like so:

M = 5;
N = 5;
[X,Y] = meshgrid(-floor(N/2):floor(N/2), -floor(M/2):floor(M/2))

X =

-2    -1     0     1     2
-2    -1     0     1     2
-2    -1     0     1     2
-2    -1     0     1     2
-2    -1     0     1     2

Y =

-2    -2    -2    -2    -2
-1    -1    -1    -1    -1
 0     0     0     0     0
 1     1     1     1     1
 2     2     2     2     2

As you can see here, the centre pixel for both X and Y are defined as (0,0).  Everywhere else has its (X,Y) co-ordinates defined with respect to the centre.

Step #2

We simply have to compute the Euclidean distance between the centre pixel to every point in the kernel.  This can be done by:

dis = sqrt(X.^2 + Y.^2)

dis =

2.8284    2.2361    2.0000    2.2361    2.8284
2.2361    1.4142    1.0000    1.4142    2.2361
2.0000    1.0000         0    1.0000    2.0000
2.2361    1.4142    1.0000    1.4142    2.2361
2.8284    2.2361    2.0000    2.2361    2.8284

Doing some quick calculation checks, you can see that this agrees with our understanding of Euclidean distance.  Moving to the left by 1 from the centre is a distance of 1.  Moving to the left by 1 then up by 1 gives us a Euclidean distance of $\sqrt{1^2 + 1^2} = \sqrt{2} = 1.4142$.  Doing similar checks with each element will demonstrate that this is indeed a Euclidean distance field from the centre pixel.  After we do this, let's compute the kernel terms without the $\alpha$ term.

kern = 1 ./ (1 + dis)

kern =

0.2612    0.3090    0.3333    0.3090    0.2612
0.3090    0.4142    0.5000    0.4142    0.3090
0.3333    0.5000    1.0000    0.5000    0.3333
0.3090    0.4142    0.5000    0.4142    0.3090
0.2612    0.3090    0.3333    0.3090    0.2612

Step #3

The last step we need is to normalize the mask so that the total sum of the kernel is 1.  This can simply be done by:

kernFinal = kern / sum(kern(:))

kernFinal =

0.0275    0.0325    0.0351    0.0325    0.0275
0.0325    0.0436    0.0526    0.0436    0.0325
0.0351    0.0526    0.1052    0.0526    0.0351
0.0325    0.0436    0.0526    0.0436    0.0325
0.0275    0.0325    0.0351    0.0325    0.0275

This should finally give you the correlation kernel that you are seeking.  You can now use this in convolution (i.e. using imfilter or conv2).

Hopefully I have answered your question adequately.  Good luck!


Kalman filtering in image processing, resources?

Here are a couple of examples from the computer vision system toolbox:
Using kalman filter for object tracking

For a more in-depth explanation the best book is Multiple Target Tracking with Radar Applications by Samuel Blackman.


From the documentation of the CMSIS-DSP FFT function:

  The FFT functions operate in-place. That is, the array holding the
  input data will also be used to hold the corresponding result. The
  input data is complex and contains 2*fftLen interleaved values as
  shown below.

It is possible to implement the FFT in such a way, that it does not require additional buffers but uses the memory of the input array for its output. This method is usually referred to as "in-place" calculation. From the above citation I think that the CMSIS-DSP library does exactly that (though we can not know for sure without the source code).

So in this case, yes, a memory size of 16 KB is sufficient.

That an in-place calculation is possible can be understood from the butterfly diagram, see for example here.


I would approach this by first thinking of this as an upsampled version of a different signal, x(m) where m=2n, which has values at every integer.

$$x(m) = ma^{|m|/2} = 2na^{|n|}$$

So x(m=2) has the same value as x(n=4) since m=2n and x(m) is non-alternating

Then all you have to do is find the z transform of the non-alternating sequence x(m), G(z) and use the upsampling property of the z-transform(http://web.stanford.edu/class/archive/ee/ee264/ee264.1052/upsampling.pdf) to know that if:

$$x(m) \rightarrow G(z)$$
$$x(n) \rightarrow F(z)=G(z^{2})$$

If my explanation isn't sufficient let me know and I'll try to help more.


Let's proceed step-by-step:

  I need some more explanation on what is gradient? How we calculate gradient? please explain me in a lucid manner I have zero knowledge about gradient.please explain me with easy explanation

a gradient relate to a measure of the derivate, that is, the coefficient of variation. If you think of a measure of position, the derivate is the speed. Now, we do not have a single point, but a field of pixel with values. Now think of it as a landscape, the luminance as the height and the gradient as the slope. It has both a magnitude (how steep it is) and a direction.

In the formulation that you cite, the slope is computed thanks to the divergence operator. In this case it equals to the sum of slope in the directions of the reference frame (that is, in x and y). It is used to allow a more general formulation in non-euclidean spaces.

  This gradient magnitude equation is what I am unable to understand please see the following equation and explain me in detail. I don't Know why they use divergence and gradient

the basic idea behind this equation is to control diffusion by the contours:

diffusion helps to smooth out data by averaging locally. one drawback is that if you do it too long, you lose all information.
contours is the abstract concept of a boundary. Think of it for instance as the contour of a object in an image or the transient in the square wave that you show. To detect the position of these contours, you may use the gradient, and specifically the amplitude of the gradient (how steep the slope is)
now comes the heuristic of anisotropic diffusion: if you diffuse, but that you control the speed of this diffusion by the borders, then, you obtain a better result.

These equations just lay down these concepts quantitatively.


If you know Linear Algebra then the process can be explained using Linear Algebra in a very simple way.

Let's say you have Orthonormal Basis for $ \mathbb{R}^{k} $ given by the Columns of $ A \in \mathbb{R}^{k \times k} $.

Now imagine having a vector $ x \in \mathbb{C}^{k} $.
How would represent the vector $ x $ in the set of coordinate of matrix $ A $, namely using the Orthonormal Basis?

Well, you'd calculate the projection (Inner Product) between any column of $ A $ to the vector $ x $ by $ {x}_{A} = {A}^{H} x $.
How would you recreate $ x $ from the columns of $ A $ using $ {x}_{A} $?
Well, you can see that $ x = A {x}_{A} $ since $ A $ is Unitary Matrix.
Basically you sum the columns of $ A $ using the weights of $ {x}_{A} $.

So, Now just replace $ A $ with the DFT Matrix and you have the whole process of DFT and Inverse DFT.

In practice, the Inverse DFT and the DFT are very similar, multiplication and summation.
The only difference is on one direction (Forward) we multiply by $ {e}^{- \frac{2 \pi j}{N} k n} $ the signal and on the backward we multiply the DFT coefficients by $ {e}^{\frac{2 \pi j}{N} k n} $.


Programming language for sound manipulation

For the first few experiments I would recommend using a scripting language like Matlab or Python. They're much easier to understand and much quicker to write than "lower level programming languages" like C++.

Matlab has a signal processing toolbox and can read and write audio files, do windowing, FFTs etc. as well as a very simple playback mechanism. Basic 2D and 3D-plotting tools are also built right into the language itself.

For Python there is the SciPy library that does all kinds of scientific computing, including wavefile handling, windowing and FFT as well. Plus, there is the Matplotlib that does all of the plotting for you.

In C++, you would have to find a UI toolkit that works for your OS, write a drawing widget to plot stuff in and then fill it with data. Plus, obviously all the memory management; wich will probably become the biggest annoyance and time-sink if you're not experienced in doing it.

To implement live playback and recording I suggest using portaudio, a library written in C. There exists a Matlab .MEX utiliy called PlayRec and there are Python bindings called PyAudio.

From my experience, if you want to use portaudio I'd suggest using Python. Compiling the .mex utility is a huge pain while Python seemed to "just work".

If you manage to get it to work but still want to do it in C++ you can then port the codebase (just a couple of lines in both Matlab or Python) while still using portaudio for the audio input/output.

For a superficial comparison of the different code-complexities you can take a look at some work I wrote last year. Each implementation does the same thing, just in a different language:

FTM filter written in C/CUDA (plus some benchmarking stuff)
FTM filter written in Matlab (plus some plotting stuff)
FTM filter written in Python (just the filter, no extras)


what symmetry does system function H(z) have if h[n] is real?

it's because $e^{j\omega}$ is such a signal that has the property:

$$ (e^{j\omega})^* = e^{-j\omega} = \frac{1}{e^{j\omega}}$$

then going back from the DTFT to Z-transform is not obvious. So it's better to start from definition of Z-transform to define the symmetry:

$$ X(z) = \sum_{n} x[n] z^{-n} $$

and apply this to a signal which is real; $x[n] = x^*[n]$

$$ x[n] \xrightarrow{ Z-transform} X(z) \implies x^*[n] \xrightarrow{ Z-transform} X^*(z^*) $$

and conclude for real $x[n]$ :
$$ X(z)  = X^*(z^*) $$

Derivation:

$$\begin{align}
\mathcal{Z} \{x^*[n] \} &= \sum_n x^*[n] z^{-n} \\ \\
& = (\sum_n x[n] (z^*)^{-n})^* \\ \\
& = X^*(z^*) \\
\end{align}
$$

As you can see, on the unit circle where $z = e^{j \omega}$, this becomes:

$$ X(e^{j\omega}) = X(e^{-j\omega})^* $$
the conjugate-symmetry property...


I don't understand what exactly you are doing, and your Nyquist frequency is only 50KHz if you're sampling at 100KHz (so you won't get 1MHz resolution).  Please clarify if you need a better answer.

Regardless, Python has several ways to compute a PSD.  One is with matplotlib.  For example:

matplotlib.pyplot as pl
# ...load in your data to dataArr...
sampFreq = 100000
Pxx, freqs = pl.psd(dataArr, NFFT=256,Fs=samplFreq,detrend=detrend_mean,window=window_hanning,noverlap=128,sides='onesided',scale_by_freq=True)

and you get a nice plot, along with the PSD data in Pxx and the frequencies in freqs.


Neural Network learning project based on 8 wave signals over 1 second at 1 sample every 10 ms ( hence 100Hz )

...a neural network that can decide wether a pattern produced by the movement of a hand near capacitive sensors is as expected, or random.
The neural network is supposed to learn himself how the different channels react, in wich order, so i don't have to tell anything to the programm concerning the physical distance between two electrodes or whatever.

At this point, I am not sure if you wanted to create a neural network that can learn how to interpret the readings from the 8 capacitive sensors to infer the position of the hand over the sensor overall, or, a classifier that can understand a specific gesture.
I think that what you are trying to do is the latter. That is, you want to move your hand in a particular gesture above the sensor "bar" and have a classifier recognising if this was one of the "known" gestures or someone just passing nearby.
If this is indeed what you are trying to achieve then you would need to look more closely into types of neural networks that can take into account time such as recurrent neural networks.
In this case, the input layer of the network would have 8 nodes, the hidden layer would have a number of neurons / layers that would depend on how long is the gesture to be learned and the output layer could well be a single node with a 0-1 output.

And another question i had, once trained, will th algorithm be exportable on an arduino ?

You can train this "offline" and then try to run the network on an Arduino, but I am not sure if you will get the same performance as you would get on a desktop computer. For two reasons:

You might run out of memory quickly on the Arduino depending on the number of nodes your network will end up having; and

The Arduino does not do floating point arithmetic. Yes, float does exist as a data type but in most boards it is emulated in software. If you have a board with a CPU that has a Floating Point Unit and the compiler can "see" it, then you might be in luck. But, without an FPU, without 32 (or at least 16) bit architecture and without much memory, then, the Arduino does not do floating point arithmetic.

There is however another way by which you can do gesture recognition with less computational complexity so that it is easier for an Arduino to run it.
Sum the outputs of all 8 capacitive sensors to produce one output and then run a "traditional" pattern recognition technique which would basically try to recognise a single waveform. In this way, you still do not have to make any additional specifications regarding the geometry of the bar.
Suppose that you configure your sensor "bar" in this way. If you run your hand from left to right, you will get an output with 8 "bumps" as your hand passes closely from each one of the sensing elements. If you move your hand from left to right and also do a "wave" in mid-air, then you will get an output with 8 "bumps" whose amplitude will vary depending on the distance of your hand from the sensors.
Obviously, the problem that this method has is that it cannot discriminate symmetrical gestures. For example, moving horizontally left to right and right to left. But, it can still discriminate between non-symmetric gestures (for example, moving left to right horizontally and moving left to right in an increasing "line" (or ramp). These two gestures would produce 2 distinct waveforms that it is fairly easy to classify.
You can do this with a number of ways from very simple to very complex with varying degrees of accuracy (and of course computational complexity).
One of the easiest (for Arduino too) would be to use a k-NN classifier.
By the way, a resource you might find useful is the Gesture Recognition Toolkit. In addition to the software itself, the wiki contains a lot of information on materials and methods it uses.
Hope this helps.


You might try padding the shorter (and even the longer) speech samples with a wide range of different (but realistic) background noise environments to create a larger training set (all to an equal training input size).  This might allow the DNN to be trained to be more robust at inference in noise.

You might even want to add variations of noise and a range of slight time-pitch variations to all your dataset samples to enlarge your training set.


Sampling signal after root raised cosine

I can visibly see from your figure that you have done the two cascaded RRC's properly. I circled the samples I would expect to see for your 00111 data sequence after RC filtering. Clearly you are using 8 samples per symbol, and a longer data sequence would continue in similar fashion after the initial delay that you have as shown.

To demodulate this, you typically would have a timing recovery loop to establish the ideal sampling locations. There are several approaches to establish timing, but one recently discussed on this forum is the use of the Gardner TED (Timing Error Detector) which would create an error signal that you could use to control a timing sample location (increment up or down a position pointer based on the accumulated error, forming a timing recovery loop. This would be done using a (discardable) training sequence with lots of transitions, typically as part of the header for the data packet that would be significantly longer than your short data set.

Another synchronizer I like is the M&M synchronizer (Mueller & Meuller) as well as excellent approaches I learned from fred harris using polyphase filters for timing synchronization (thank you fred!).

Below are links to recent discussions on the Gardner TED, but in particular I want to copy here the figure for the eye diagram of the transmitted symbols after they have gone through a root raised cosine (RRC) filter ("Before RRC Filter") the received signal after it has been passed through a second RRC filter as you are doing. Conveniently this was done with same alpha= 0.25 roll-off factor you have used, and what it shows is the continuous time representation of every possible trajectory between symbols based on previous symbol patterns. In your case, you have eight samples per symbol, landing on one of these trajectories.

Links to recent discussions on the Gardner TED:

Location of Matched Filter

Gardner Timing Recovery for Repeated Sybmols


Let's say the audio is given by $ x $.
Let's define Low Pass Filter by $ h $.

Whet you need to do is:

$$ y = x + \alpha (h \ast x) $$

Where $ \alpha $ is the amplification level of the Low Pass data.
Once you set which frequencies you want to amplify it is easy to design.


MATLAB, error while trying to implement a Designed Filter object

The issue has been resolved. This is due to the absence of the filter object fields:

                    FilterFunction:
                  FiltfiltFunction:
GetCoefficientsForAnalysisFunction:

In MATLAB 2015a these fields are absent from the filters objects. But in MATLAB 2016a (which is I am using, these fields seem to be required.) As the files were imported after creation in MATLAB 2015a, that seems to be the problem. Recreating, the fields in MATLAB 2016a solved the issue. The fields have also been created.

The code ran to check:

filt = D_sig50(1);
sf = struct(filt) ;
sf.FilterFunction

Created in MATLAB 2015a:
output:

                    FilterFunction: []
                  FiltfiltFunction: []
GetCoefficientsForAnalysisFunction: []

These fields in the filter object are empty (after loading these filters in MATLAB 2016a)

Created in MATLAB 2016a:

                    FilterFunction: @filterIIR
                  FiltfiltFunction: @filtfiltIIR
GetCoefficientsForAnalysisFunction: @getCoeffsForAnalysisIIR

The fields are now present in the filter object, and the code works now.


Explain the Adaptive Part of Adaptive Algorithms - Kalman Filter and Least Mean Square / Constant Modulus

Adaptive Filters are called "Adaptive" when they can adapt to changes in data.
In the filters you mentioned above, which are part of the Linear Filters family the property means their coefficients are changing over time.

Linear Filters are basically weighing and summing the data.
For instance, given no prior information on data you may want to have exact weight for any data given.

Yet in most Signal Processing use cases we'd like to (Or might want to) give higher weight for each data sample according to its properties such as SNR, How good it fit the model, etc...

In the classic Kalman Filter the coefficients are set according to a model matrix - $ P $.
For Kalman filter this matrix stands for the Covariance of the estimation.
It changes according to 2 other matrices which are properties of our model and data - $ Q $ - The model confidence level matrix, $ R $ the data confidence model (For Gaussian data it is basically the SNR).

Since the algorithm support the cases those are changing over time it is called adaptive.

For instance, in the case of tracking a target using RADAR the algorithm can change the matrix $ R $ according to the SNR of the measurement.

The Recursive Least Squares have similar properties (Though it depends only on one factor, something similar to the matrix $ R $ of the Kalman Filter).
One could even see the RLS as a private case of the Kalman Filter (Or the Kalman Filter as a generalization).

There are many extension for the Kalman Filter for many cases.
The UKF (Unscented Kalman Filter) you mentioned is built to handle Non Linear cases of the Model.
But adaptive concept stays the same, the weights are changing according to properties of the Data and the Model.


2 mm period between samples is 500 samples per meter.  That's your sampling frequency.

Matlab's digital filter design functions (is that what you're using?) take a normalized frequency from 0 to 1, where 1 is the Nyquist frequency (1/2 the sampling frequency; the highest spatial frequency you can measure at that sample rate).

So if you wanted to filter at a spatial frequency of 100 cycles per meter (features that cycle every 1 cm), for instance, you would note that the Nyquist frequency is 250 cycles per meter, and enter a value of 100/250 = 0.4 for the normalized frequency.