Half of fundamental frequency in Fourier transform I understand why there are peaks in integer multiples of the fundamental in the Fourier transform. But why is there often a peak in half the fundamental?

I am testing for the guitar specifically. I am doing an STFT on an guitar audio monophonic signal and then get the 4 strongest frequencies in each note. I always get the fundamental in those 4 and of course some multiples of it but also half of it.

When I play A3 (220Hz) in my guitar, the fifth string which is A2 (110Hz) also vibrates a bit: it is what is called Sympathetic resonance. Besides other non-linear effects, this could be the case.
This phenomenon also occurs sometimes when recording bat calls, you can  observe resonance at half the fundamental frequency.

Here, the fundamental frequency is at 65 kHz (upper orange tag) and there is resonance, probably generated by the microphone, at half that frequency (lower orange tag).


Plotting the Phase Response I would appreciate it very much if someone would be able to provide some clarity on plotting phase responses.

For instance, given that the frequency response of a filter can be written as $$H(\exp(j*\Omega))=\cos(\Omega)$$ How would I go about determining the phase response?

The answer key demonstrates a piecewise solution:

This is what I understand so far:

$\cos(\Omega)$ is strictly real, so the phase carries values of $0$, or $\pm\pi$
Between the values of $-\pi/2$ and $\pi/2$, $\cos(\Omega)$ is positive, so the phase is $0$

I understand since the function $\cos(\Omega)$ is strictly real, phases where the function is not positive should either be $\pi$ or $-\pi$ (as opposed to $-\pi/2$ or $\pi/2$, for when it's imaginary), but how do you know for which range it's $\pi$, and for which range it's $-\pi$?

This might be a stupid question, but I've spent a lot of time thinking it through, so I'm thinking I'm missing something, and would appreciate any clarity.

Thank you so much!

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.


Find smoothed first derivative from signal with noisy slope How can I get filtered first derivative from a noisy signal that has slowly changing slope in form of y=kx+b? k can slowly change in time, and I want to estimate its value.

I have tried 3 different approaches:

Take derivative as dx(i) = (x(i)-x(i-99))/100
Smooth with sliding mean window = 100, then take derivative as dx(i) = (x(i)-x(i-99))/100
Simple IIF filters (e.g. y(i) = 0.99*y(i-1) + 0.01*x(i), then take derivative as dx(i) = y(i)-y(i-1) and again filter with similar IIR, e.g. dy(i) = 0.95*dx(i-1) + 0.05*dx(i)

Problems:

Least-squares, regression and FIR filters (except rectangular window) have high computational cost, since I have to translate it to micro-controller with no DSP. That is why I can use only rectangular windows and IIR filters (they have low order).
If I find first derivative first, then smooth, it will be very noisy. So, I should smooth the original signal first, then find derivative from smoothed signal (and perhaps smooth the derivative again!).
I should play with filter parameters manually and it is hard to understand the frequency response of the whole system.

Question:

Maybe there is a single special (optimal?) IIR filter for this specific problem - finding smoothed first derivative from signal with noisy slope?

I think least squares is going to be the best approach, and that's not going to be that computationally expensive (I think! Please correct me if I'm wrong).

The gradient can be estimated from a sliding window of your data using:
$$
\hat{k} = \frac{\sum (x_n - \bar{x})(y_n - \bar{y})}{\sum (x_n - \bar{x})^2 }
$$
where the sum over $n$ is taken over the window.

If I implement an example in R (code below) then it seems to do the right thing. The picture shows the noiseless data (red) and noisy data (black) with a change in $k$ at half-way.  The window length I've chosen is 20 points, and the output shifts at about the right point.

This estimate may still be too noisy for you; in that case, increase the window length. This picture shows what happens when window_length is increased to 100.

R Code Only Below

#42364

T <- 1000
k1 <- 0.006998
k2 <- 0.0188728493
k_true <- c(k1*rep(1,T/2), k2*rep(1,T/2))
c_true <- c(2.90238432*rep(1,T/2), 2.90238432*rep(1,T/2) + (k1-k2)*(T/2+1))
t <- seq(1,T)

y_true <- k_true*t + c_true

y_noisy <- y_true + rnorm(T)*0.1

k_hat <- rep(0,T)

window_length <- 20
for (t_idx in t)
{
  oldest_idx <- max(t_idx-window_length, 1)
  window_idx <- seq(oldest_idx, t_idx)
  x_bar <- sum(window_idx)/length(window_idx)
  y_bar <- sum(y_noisy[window_idx])/length(window_idx)

  xy <- sum((window_idx - x_bar)*(y_noisy[window_idx] - y_bar))
  xx <- sum((window_idx - x_bar)*(window_idx - x_bar))

  k_hat[t_idx] <- xy /xx
}

par(mfrow=c(2,1))
plot(t, y_noisy, type='l')
lines(t,y_true, col='red') #, lwd=3

plot(t,k_hat)
lines(t,c(k1*rep(1,T/2), k2*rep(1,T/2)), col='red')


GPS timing and computational complexity I have been using Garmin GPS device for years. Just curious which algorithm is used in the device:

Navigation equations that can calculate both GPS device's position and its clock bias. In this way, the time display should be accurate with the clock in the satellite. However, the computational complexity could be very high. See https://en.wikipedia.org/wiki/Global_Positioning_System.
Use TDOA method similar to what is explained at https://en.wikipedia.org/wiki/Multilateration. The computational complexity could be much lower, but you will not be able to obtain the receiver's time bias.

Note: This question is based on the assumption that the satellites' positions and TOFs have already been extracted from the received satellite signals.

Computational complexity lies in the correlator banks you need to despread the GPS signal, not in the position and time calculation (these are not separate things. You can't have one without the other). In fact, the TDOA system is very likely far more computationally intense, since it needs to do more cross-correlations.

You need to use the first method; 2. can't work, you'd need a clock that is about as stable and thus as large, power-hungry and expensive as that on the GPS satellites. Even with that clock, it wouldn't work well, because you still need to account for a lot of effects, including special relativity, that you can't solve via simple TDOA.

Computational complexity is still rather high, indeed, but really, not a problem since the late 1990s, where highly integrated chips for GPS became the norm. What it does do is drain the battery of your smartphone!


Signal bandwidth and channel bandwidth I have a channel bandwidth $B_c$ and a signal bandwidth $B_s$ with $B_s<<B_c$. What is the problem when the signal bandwidth is smaller than the channel bandwidth? For example, my channel bandwidth is $B_c=20$ kHz and my signal bandwidth is $B_s = 10$ Hz.

No problem at all if $B_s \ll B_c$.
However, $B_s = 10\text{Hz}$ is considerably small that makes synchronization very difficult especially for low-cost receivers.
For example, Sigfox Ultra Narrow Band uses $10\text{Hz}$ bandwidth (are you studying the standard?) that requires precise but expensive oscillators like TCXO. Even with that, there must be random access enhancement to deal with the frequency imprecision.
There is also some effort to design transceivers without TCXO for such narrow bandwidth. Check e.g. Foxy.


Audacity Spectogram view: how to see exact frequency and intensity under the cursor? In the Spectogram view, simply clicking anywhere with the selection tool gives a vertical line to the top ruler, and shows at the bottom the exact sample, H:MM:SS, or other time coordinate you'd want.

I'm looking for the exact same functionality for frequency (Hz, or say musical pitch such as C7+14 cents) and intensity (-25dB).

Should Audacity not be able to do this, is there another tool that can read a WAV and do it?  I've looked at Spec and it can't but it should be easy enough to modify (I'm a software engineer).

What you are asking is possible using heuristics but impossible in general because of the uncertainty principle. A single sample is just a value, it does not even make sense to ask about its frequency or frequencies. Roughly speaking (though not strictly true), the concept of "frequency" necessarily implies some form of repetition, thus it only applies to a range of samples, not a single sample. You can, of course, look at a fixed-length window around the sample you are interested in, but the result (both the number of frequencies and their accuracy) will depend a lot on the window length.

I think Sonic Visualiser has pretty decent heuristics for instantaneous frequencies.


3dB-Cut off frequency of moving average I've read the thread: Cut off frequency of moving average

And I use the second answer in my algorithm to calculate the 3dB cut off frequency of my filter, which works great, as my filter length is usually above 300. I verified it with the step response.

But I would like to have a source or derivation for this formula.

I tried by hand with taylor series stopping after the second and third term. I come close but not exactly to the formula and mapple gives me a valid but extremly complicated result.

Hope you guys can help.

kind regards

Slev1n

Consider a zero-phase moving average of length $N$:

$$\text{y}[n] = \begin{cases} \displaystyle\frac{\text{x}[n] + \displaystyle\sum_{k=1}^{\frac{N-1}{2}}\left(\text{x}[n+k] + \text{x}[n-k]\right)}{N},&n\in\mathbb{Z}&\text{for }N\text{ odd}\\
\displaystyle\frac{\displaystyle\sum_{k=1}^{\frac{N}{2}}\left(\text{x}[n+(k-\frac{1}{2})] + \text{x}[n-(k-\frac{1}{2})]\right)}{N},&n+\frac{1}{2}\in\mathbb{Z}&\text{for }N\text{ even}
\end{cases}$$

Even-length filters operating on discrete sequences with integer time indices cannot be zero-phase. We have circumvented this by enabling the output time indices to always have a fractional part of $\frac{1}{2}$, in case of $N$ even. As a real-world example, if the input was sampled every midnight, the zero-phase moving average of even length would be calculated for every noon. This unusual indexing conveniently gives the same zero-phase form of frequency response $F_N(\omega)$ for both $N$ odd and $N$ even:

$$F_N(\omega) = \begin{cases}
\displaystyle\frac{\displaystyle1 + \displaystyle\sum_{k=1}^{\frac{N-1}{2}}\left(e^{ik\omega}+e^{-ik\omega}\right)}{N}&\text{for }N\text{ odd}\\
\displaystyle\frac{\displaystyle\sum_{k=1}^{\frac{N}{2}}\left(e^{i(k-\frac{1}{2})\omega}+e^{-i(k-\frac{1}{2})\omega}\right)}{N}&\text{for }N\text{ even}
\end{cases}\\
= \frac{\sin(\frac{N\omega}{2})}{N\sin(\frac{\omega}{2})}$$

Unfortunately the frequency response has no symbolic solution for the -3 dB cutoff frequency $\omega_c$, such that:

$$F_N(\omega_c)=\sqrt{\frac{1}{2}}.$$

Strictly speaking $\sqrt{\frac{1}{2}}$ is about -3.01 dB, but I think that's what people mean when they say -3 dB, because otherwise it is just an arbitrary number. An approximate frequency response $\hat{F}_N(\omega)$ uses an integral instead of a sum:

$$\hat{F}_N(\omega)=\frac{1}{N}\int_{-\frac{N}{2}}^\frac{N}{2}e^{ik\omega} = \frac{2\sin(\frac{N\omega}{2})}{N\omega}$$

The main lobes of the true (sum) and the approximate (integral) frequency responses converge at large $N$:

We can prove the convergence by introducing functions $G_N(\chi) = F_N(\omega)$ and $\hat{G}_N(\chi) = \hat{F}_N(\omega)$ with the argument normalized such that $\omega = \frac{2\pi\chi}{N}$, bringing the first zero of both functions to $\chi = 1$:

$$G_N(\chi) = \frac{\sin(\pi\chi)}{N\sin\left(\frac{\pi\chi}{N}\right)}\\
\hat{G}_N(\chi) = \frac{\sin(\pi\chi)}{\pi\chi}\\
\lim_{N\rightarrow \infty}G_N(\chi) = \frac{\sin(\pi\chi)}{\pi\chi}$$

$G_N(\chi)$ is known as the $N$-periodic band-limited impulse train. Its limit at large $N$ and the function $\hat{G}_N(\chi)$ are both the $\text{sinc}$ function. Unfortunately the -3 dB cutoff frequency has no symbolic solution in the approximation $\hat{F}_N(\omega)$ either. For different $N$, the approximation only differs from the $N = 1$ approximation by a mapping $\omega \rightarrow \omega N$, so it suffices to solve the approximate -3 dB cutoff frequency $\hat\omega_c(N)$ numerically for $N = 1$:

$$\frac{2\sin(\frac{\hat\omega_c(1)}{2})}{\hat\omega_c(1)} = \sqrt{\frac{1}{2}}\\
\Rightarrow\hat\omega_c(1) = 2.78311475650302030063992,$$

giving the approximate cutoff frequency for arbitrary $N$:

$$\hat\omega_c(N) = \frac{\hat\omega_c(1)}{N} $$

That appears to be another, simpler approximation than Massimo's. For your $N > 300$ there should be no problem using it. Massimo's and this answer's constants are related by:

$$\frac{\hat\omega_c(1)}{2\pi} = 0.442946470689452340308369.$$

I looked a bit further and found that Massimo approximates $F_N(\omega)$ with $\hat{F}_M(\omega)$, choosing $M$ such that the (limits of the) second derivatives of the frequency response and the approximation match at $\omega = 0$:

$$F''_N(\omega)=\frac{\sin(\frac{N \omega}{2}) \left(2 \sin(\omega) \cos(\frac{\omega}{2}) + (N^2 - 1) \sin(\frac{\omega}{2}) (\cos(\omega) - 1)\right)}{2 N (\cos(\omega) - 1)^2}  - \frac{2 N \sin(\omega) \sin(\omega/2) \cos(\frac{N \omega}{2})}{2 N (\cos(\omega) - 1)^2}\\
\hat{F}''_M(\omega)=\left(\frac{4}{M \omega^3} - \frac{M}{2 \omega}\right) \sin\left(\frac{M \omega}{2}\right) - \frac{2 \cos(\frac{M \omega}{2})}{\omega^2}\\
\lim_{\omega\rightarrow 0}F''_N(\omega) = \frac{1-N^2}{12}\\\lim_{\omega\rightarrow 0}\hat{F}''_M(\omega) = -\frac{M}{12}\\
\lim_{\omega\rightarrow 0}F''_N(\omega) = \lim_{\omega\rightarrow 0}\hat{F}''_M(\omega)\\
\Rightarrow M = \sqrt{N^2 - 1}\\
\Rightarrow \hat\omega_c(M) = \frac{\hat\omega_c(1)}{\sqrt{N^2 - 1}} $$

This improves the approximation at small $\omega$ which includes the -3dB cutoff point, especially at small $N$:

Massimo's approximation always overestimates the cutoff frequency (see the error comparison), leaving room to improve it by altering the constant $1$. The error is the largest for $N = 2$. If its error is constrained equal to the (currently second largest) error at $N = 3$, we get an even better but just as cheap approximation:

$$\hat\omega_c(N) = \frac{2.78311475650302030063992}{\sqrt{N^2 - 0.8603235992780290790596}}$$

This and other tweaks of the constant, like Matt's constant $0.863031932778066$, work surprisingly well for large $N$ (see the error comparison). For large $N$ the error falls by a factor of 1000 for every increase of $N$ by a factor of 10. The explanation for these things is that the true cutoff frequency as a function of $N$ has a Laurent series:

$$\omega_c(N) = \displaystyle\sum_{k=0}^{\infty}\frac{a_k}{N^k},\ a_k = 0\text{ if }k\text{ even},$$

and the approximation and its Laurent series are:

$$\hat\omega_c(N) = \frac{a}{\sqrt{N^2+c}}\\= \frac{a}{N} - \frac{ac}{2N^3} + \frac{3ac^2}{8N^5} - \frac{5ac^3}{16N^7} + \frac{35ac^4}{128N^9} - \ldots,$$

such that:
$$a_1 = a = 2.78311475650302030063992\\
a_3 \approx -\frac{ac}{2}$$

If the approximate match in the $N^{-3}$-term was made exact, the approximation error should decrease by a factor of $10^5$ for an increase of large $N$ by a factor of $10$. The coefficients $a_k$ of the Laurent series $f(x) = \sum_{k=0}^{\infty}\frac{a_k}{x^k}$ of a function $f(x)$ as $x\rightarrow\infty$ can be found iteratively by:

$$f_0(x) = f(x)\\
a_k = \lim_{x\rightarrow\infty}f_k(x)\\
f_{k+1}(x) = \left(f_{k}(x)-a_k\right)x$$

When we do not have $f(x)$ in symbolic form, but can solve it numerically to any precision for very large $x$, we can do the equivalent of the above procedure numerically. The following Python script that uses SymPy and mpmath will calculate a given number (here 10) of the first coefficients $a_k$ in the desired precision for the Laurent series of the true cutoff frequency:

from sympy import *
from mpmath import *
num_terms = 10
num_decimals = 24
num_use_decimals = num_decimals + 5 #Ad hoc headroom
def y(omega):
    return sin(N*omega/2)/(N*sin(omega/2))-sqrt(0.5)

a = []
h = mpf('1e'+str(num_decimals))
for i in range(0, num_terms):
    mp.dps = 2*2**(num_terms - i)*num_use_decimals*(i + 1) #Ad hoc headroom
    N = mpf('1e'+str(2**(num_terms - i)*num_use_decimals))
    r = findroot(y, [2.7/N, pi/N]) #Safe search range
    for j in range(0, i):
        r = (r - a[j])*N
    a.append(r)
    mp.dps = num_decimals
    print 'a_'+str(i)+' = '+str(nint(r*h)/h)

On my computer the program runs for about 7 minutes. It prints the following, showing that the Laurent series consists only of odd negative powers:

a_0 = 0.0
a_1 = 2.78311475650302030063992
a_2 = 0.0
a_3 = 1.20103120331802052770586
a_4 = 0.0
a_5 = 0.767601805674195789947052
a_6 = 0.0
a_7 = 0.537174869947196574599614
a_8 = 0.0
a_9 = 0.387986305863698148870773

These numbers, shown to 24 decimal places, are not from an approximation in the sense that the Laurent series is unique; there is no other Laurent series that equals $\omega_c(N)$. Using only $a_1$ and $a_3$, a simple two-term truncated Laurent series approximation can be constructed:

$$\hat\omega_c(N) = \frac{2.78311475650302030063992}{N} + \frac{1.20103120331802052770586}{N^3},$$

and by $c=-\frac{2a_3}{a_1}$ the approximation:

$$\hat\omega_c(N) = \frac{2.78311475650302030063992}{\sqrt{N^2-0.863084212040982824494534}}$$

Both have $1/N^5$ error decay at large $N$, see error comparison, columns h) and i) respectively. A longer truncated Laurent series with more terms from the script's output decays even faster, $1/N^{11}$ for the 5-term approximation at column j) in the error comparison.


How to get a real signal from complex ifft values? I am trying to understand the difference between DCT and FFT in the context of compression. There is a nice practical example given here. (near the end of the page)

When I try it in matlab I am stuck at plotting the ifft values (which are complex). The plot command cannot plot the complex signal. One natural question comes to my mind, when is the ifft value complex, and when is it real? If it is complex, is it fine to discard the imaginary parts and keep only the real part?

What is the right way to reconstruct the real signal x from ifft complex values?

You can plot the real and imaginary components from the FFT result separately (or the magnitude and the phase results separately, as two plots).  Then you can recombine the data from the two plots into a complex spectrum required to do a reconstruction IFFT.  Both components are required, otherwise you've thrown away half the information about your signal, and can't reconstruct it (unless the original signal was purely symmetric or purely anti-symmetric to begin with).

The plot data isn't really doubled when you do two plots because the result of an FFT of strictly real inputs in conjugate symmetric.  Thus you can create the input for your reconstruction IFFT from only the first half of the data for the real plot and for the imaginary plot (by conjugate mirroring it to the other half).


Derive the Forward Euler substitution for transfer function In the book "The control handbook. Volume 1 " by Levine, the author shows that the transfer function:

can be aproximated and discretized in the transfer function:

using the forwar euler integral aproximation. In this way Levine shows that if i have a transfer function $G(s)=\frac{1}{s}$ i can simply substitute $s=\frac{z-1}{T}$ to apply the Forward Euler aproximation.
At this point the author doesn't show how to generalize this particular case to a generic transfer function $H(s)$. Are you able to show me ?

To the extent you can factor the transfer function into individual integrator sections of the general form $\frac{1}{s}$ you can make this substitution, which is an approximation of the Matched-$z$ Transform where you substitute every $s$ for $s=\frac{\ln(z)}{T}$. (map from $s$ to $z$ using $z =e^{sT}$).

This results in first order forms given by $H_\mathrm{a}(s) = k_\mathrm{a}\frac{s-b}{s-a}$ being replaced with $H_\mathrm{d}(z) = k_\mathrm{d}\frac{z-e^{bT}}{z-e^{aT}}$ where $k_\mathrm{a}$ and $k_\mathrm{d}$ are fixed gain constants. $k_\mathrm{d}$ can be determined by using $H_\mathrm{a}(s=0) = H_\mathrm{d}(z=1)$ for low-pass forms or $H_\mathrm{a}(s=j \pi/T ) = H_\mathrm{d}(z=\pi)$ for high-pass forms.

For all pole forms this is also equivalent to what you would get using the mapping from $s$ to $z$ using the method of Impulse Invariance. The Matched-$z$ Transform is not recommended for higher order polynomials with zeros and poles since neither the time or frequency domain response is preserved (The method of Impulse Invarance or Bilinear Transform would more often be used, each with their own advantages).

This post summarizes various mapping techniques further: Are there alternatives to the bilinear transform?


What are some face detection (not recognition) algorithms suitable for limited (embedded) processors? I am looking at a hardware based platform for surveillance class of applications. I want to identify faces of people as they appear in the scene.

I don't intend to really apply face recognition (then or later). Also orientation of faces can be assumed straight for practical purposes. Face detected once needs no further processing.

Is there is a fast (efficient) algorithm to detect faces from the live video source?

We know that the problem is solved because our cell phones and digital cameras can do it. You might be able to find out what algorithm they use by trawling the patent databases.

Computer vision libraries like OpenCV and SimpleCV offer face detection as a standard feature. See for example Near realtime face detection on the iPhone w/ OpenCV port. If you want to DIY see Simple Face-detection Algorithm Based on Minimum Facial Features.


Explanation of binning in frequency analysis post-FFT I am creating an Android app which records sound for t seconds at a sampling rate of 44.1kHz with a buffer size of 8192 (16 bit mono samples). I need to plot a graph of amplitude against frequency, and so first I need to run an FFT over the buffer each time it is full (which happens multiple times every second, more so with an increasing sampling rate). However, after I run the FFT I understand I need to 'bin' the frequencies and then apply an A-weighting to each bin. From the Nyquist criterion, I understand that I only need to look at the first half of the frequencies, and disregard the first corresponding to the DC component. However, I am unsure of what 'binning' actually means and would like some clarification as to what it is, and how to choose what goes in what bin.

Please advise :) .

In its simplest form, the 'binning' process consists in summing the energies (squared magnitude) within groups of adjacent FFT values. This will give you the total energy in a set of disjoint frequency bands.

A more elaborate form consists in using overlapping triangle filter banks - you compute a weighted sum of the energy in a range of FFT bins to get a number which can be interpreted as the energy measured at the output of a band-pass filter of a given center frequency/width. For example, one can use this technique to build a visualization which uses a 1/3 octave scale on the frequency axis.


Why would someone need a Spatial Light Modulator with a maximum phase shift $> 2\pi$ I am working with Spatial Light Modulators for a project and I was wondering why does the producer think someone would need a maximum phase shift $>2 \pi$ since the device works as a monitor and every phase-element monitored on the device could be reduced between $[0,2\pi]$ using a modulo?

Marko

That will be useful at least for pulse shaping, because then it is not just a single wavelength that you are dealing with, and as modulo $2\pi$ is actually implemented as modulo of displacement, it would work only for a single wavelength, creating jumps in the phase for nearby frequencies.


signals comparison metrics Am new to signal processing and was wondering when given two signals, what are the widely used statistical analysis methods to understand the relationship between them?

There are for instance coherence and (cross)correlation which are commonly used mostly because of their computational/practical attractiveness. If you need to go beyond second order relationship, that is higher order dependencies, you can look at (cross)cumulants. From information theory you may want to look at (cross)entropy and mutual information.


Find transfer function from root locus and step response diagram? I am given the response of a step of magnitude of 3 and the root locus and I have to find the transfer function of the system. The function I find gives me the step response(magnitude of 3 again) of the last diagram.

I'm a beginner at this so I've done something stupid probably but I have trouble finding answers regarding control engineering on the internet. This is what I tried doing:
I found the poles and the zeros from the root locus. $z=-5,+4, p=-6,-10,-3$. I think my transfer function is given from this formula but I'm not sure if we have an $H(s)$ in the feedback and it is not stated : $$ T(s)=  \frac{KG(s)}{1+KG(s)} $$
From the poles and the zeros my open-loop transfer function $G(s)$ is :
$$ G(s)= \frac{(s+5)(s-4)}{(s+10)(s+6)(s+3)}$$ Doing the calculations I find : $$ T(s)= \frac{Ks^2+Ks-20K}{s^3+(K+19)s^2+(108+K)+180-20K}$$
From the step response(final value is 4) and the final value theorem I find $\frac{-20K}{180-20K}=-4/3\implies K=5.14$.  I divided 4 by 3 because of the step magnitude. With this K the step response is the one in the third diagram.It's close to the first one but it's not the one I'm looking for.

What am I missing here?

The transfer function is

$$H(s)=\frac{K(s - 4)(s + 5)}{(s + 3)(s + 6)(s + 10)}=\frac{K(s^2 + s - 20)}{s^3 + 19s^2 + 108s + 180}$$
So we need to find $K$.

The LT of input step is $\frac{3}{s}$. Using the FVT for the step response $g(t)$:

$$\displaystyle\lim_{t\to \infty}g(t)=\displaystyle\lim_{s\to 0}s\cdot \frac{3}{s}\cdot\frac{K(s^2 + s - 20)}{s^3 + 19s^2 + 108s + 180}=-4$$
or
$$\frac{3\cdot K\cdot (-20)}{180}=-4$$
which gives $K=12$.
Then use the following code to plot the step response

K= 12;
sys = tf(K*[1 1 -20],[1 19 108 180]);
step(3*sys)

So there is no feedback.


Povey window formula In Kaldi, "povey" is a window  made to be similar to Hamming but to go to zero at the
edges, it's:

$$g(n) = \left(\frac{1-\cos\left(\frac{n}{N2\pi}\right)}{2}\right)^{0.85}$$

Please, can anyone mathematically explain why we are considering raising to the power of 0.85? Why not 0.8 or 0.9?

It is difficult to explain the exact reasoning behind this exponent quantitatively but more generally, raising a function to some power $\alpha$, sharpens (for $\alpha>1$) or "dulls", (for $\alpha<1$) the curve's slopes.

In this particular application, raising to a power, changes the slopes of the window function which, multiplied by some signal $x(n)$, will eventually shape its spectrum. This is the purpose of the window function.

(In the application of the Raised cosine filter for example, the "pulse shaping" is more intuitive).

Therefore, the short answer is that the coefficient is 0.85 because that coefficient gave the designers the sort of spectrum lobe shaping they were after.

Now, for a more practical demonstration, here is the exact same curve but at different exponents. On the left it is the window function itself, on the right it is the window's spectrum.

There are a few things worth noticing here:

I have put the windows and their spectra in the same graph because the differences for small variations of the exponent are subtle.
Because the window function is multiplied with the signal in the time domain, the resulting signal (from the application of the window) is the convolution of the signal's spectrum with the window's spectrum. So what you see on the right, the spectrum of the window function, is not multiplied with the signal. It undergoes the slightly more complicated operation of convolution. Given a signal $x(n)$ and a window function $h(n)$, the resulting signal is $y(n) = x(n) \cdot h(n)$. But this will result in $Y = \mathcal{F}(x) * \mathcal{F}(y)$ where $\mathcal{F}$ denotes the Discrete Fourier Transform.
This convolution operator is what causes the shaping of the spectrum. To understand broadly what effect does the window have on the spectrum, you need to understand the relationship between the time domain and the frequency domain. Broadly speaking, a "sharp" function in the time domain has a "wide" frequency domain representation and vice versa. For more information about why this is, consider the shape of the sinc filter. Now, if you convolve your spectrum with a "wide" function, you basically smooth the sharp transitions between harmonics and vice versa. But, to produce a wide convolving function in the frequency domain, you need a sharp function in the time domain (and vice versa). To achieve this, you raise the exponent.
This might be a little confusing at first but once you get the relationship between the two domains, along with the convolution theorem, you will be in a better position to understand what a window function achieves and how.

Hope this helps.


What filter is better with -1/3 2/3 -1/3 or with -1/4 1/2 -1/4 coefficients? Professor asked me what filter is better -1/3 2/3 -1/3 or -1/4 1/2 -1/4, the answer was second. But what those coefficients means? It was during discussion of linear filter but i'm not sure was it about it.

Both filters can be re-written in the form:
\begin{equation}
H(z) = \alpha ( 1 -2z^{-1} +z^{-2} )
\end{equation}
Clearly, in the first case, $\alpha=-\frac{1}{3}$ and in the second $\alpha=-\frac{1}{4}$.

In other words, the shape of the filter's frequency response is the same, the only difference is the scaling.

I can think of two reasons of why the second filter may be preferred.

In the case of the second filter, the maximum gain of the filter (occuring at the Nyquist frequency) is normalized to $1$, whereas the first filter has a maximum gain of $\frac{4}{3}$.
The other more compelling reason however, is probably that the coefficients of the second filter can be expressed as powers of two ($\frac{1}{4}=2^{-2}$, $\frac{1}{2}=2^{-1}$) and therefore multiplication by these coefficients can be quickly implemented as bitwise operations on a computer or digital signal processor.


How we can quantize a sampled signal in MATLAB? I have a continuous time signal $x(t)=\sin(2πft)$ where $0 \leq t \leq 3$.

I want to sample the continuous time signal and then quantize that sampled signal and then plot both sampled and quantized signals in MATLAB.

I have written code for sampling section/part but I don't know how to go with quantization part/section.

My MATLAB code is follow:

sampleRate = 500;
samplePeriod = 1/sampleRate;
signalFreq = 20;
nT = 0:samplePeriod:3;
signal = cos(2*pi*signalFreq*nT);

figure
stem(nT, signal)
grid on
xlabel('Time (sec)')
ylabel('Amplitude')

First you will need to determine the number of quantization levels. I am going to assume a power of two for digital convenience's sake.

 nbits = 8 % 256 qantization levels
 qLevels = 2^nbits

The next step will be to scale your signal to have the same magnitude as your number of bits.

 signalMin = -1
 signalMax = 1
 scalingFactor = (signalMax-signalMin)/qLevels
 signal = signal / scalingFactor

This gives you a signal ranging from -128 to 128.
You now have the choice of four functions to use:

floor()  %round down to the nearest integer
round()  %round to the nearest integer
ceil()   %round up to the nearest integer
fix()    %round towards zero

I will use round(), then scale the signal back to its original magnitude

signal = round(signal)
signal = signal * scalingFactor

You already seem to have the ability to graph down. I will leave graphing the function up to you.

Edit: To respond to abtj's comment about quant():
I wasn't familiar with quant, but it seems it would work just fine.
The second argument needs to be the value of the least significant bit. This is the same as the scalingFactor as calculated in the code above.

The scalingFactor is simply a way to scale the original signal to the size of the quantization. i.e. scale a signal from -1 to 1 volts to ±8 for a 4 bit quantization. This is to make a function like round() useful, which will only round to integer values.
This is done by taking the range of the original signal (signalMax-signalMin) and dividing it by the number of quantization levels desired.

Please point out any bugs in my code; I don't have access to MatLab to test right now.

Edited: Swapped the * and / when scaling the signal.
Edited: Added Parentheses per A_A's comment.


what is probability density function (PDF) of periodic signal Anybody knows what is the best pdf for periodic signal ?

I have signal with $480$ samples $\left(16000*0.03=480\right)$ but I could not find the exact pdf for that. sometimes it appear like 2 mixed Gaussian.

There is no best pdf for a periodic signal. There is also no way to find the 'exact' pdf of a measured signal. What you have to do is to measure the pdf from the data. Use a histogram to approximate the pdf of your data. Define a number of intervals within the amplitude range of your signal and simply count the numbers of data samples per interval. This gives you the histogram which is an approximation of the pdf of your signal.


Unit Problem in Designing a Filter for a Given Auto Correlation Function Given a WSS process with the following Auto Correlation Function:

$$ r\left ( \tau \right ) = {\sigma}^{2} {e}^{-\alpha \left | \tau \right |} $$

The Laplace Transform would be:

$$ R \left ( s \right ) = \mathfrak{L} \left \{ r \left ( \tau \right ) \right \} = \frac{-2 \alpha {\sigma} ^ {2}}{\left ( s - \alpha \right ) \left ( s + \alpha \right )} $$

Hence the filter would be of the form:

$$ H \left( s \right) = \frac{c}{s + \alpha} $$

My question is about units.
In the Auto Correlation Function the units of $ \alpha $ are [Hz].
While in the filter form, assuming $ s = j \omega $ the units are [Rad / Sec].
How this conflict can be resolved?
What a I missing?

Thanks.

Radians are considered to be dimensionless.  See Are angles dimensionless? and Dimensionless quantity.  They are considered to be pure numbers like pi.

So $\alpha$ is in Hz, which is a measure of 1/second, and $s$ is also considered to be measured per second.


Recommended literature to learn essential aspects of signals and its propagation Where can I learn about signal's features. For instance mean, power, amplitude, frequency, energy,  how they propagate in space or any other medium?.

What power means its amplitude increases and decreases to draw the sine wave form through a bit of time. at stable value (forming logical message) depending on medium resistance and transmitting signal power, human code model frequently repeated (hi) started with character ends with another completed word (world) another word stared with a header and ends with trailer completed word.

out of electricity network that it out of control compared with communication signals networks even both of them propagate as a flood of electrons in cables, i need to know its shape in a medium, MW or Electromagnetic waves or in cables.
voice fades away at the moment you shot a word. you already carried it on a carrier to reach ear of who you talk to but at low power so that its power and frequency decreased and lose its shape at a point seems to be disappeared or absorbed to be converted to another shape of energy.

We carries it on a carrier varies depending on the medium converting carrier to cover all values voice reach as possible, more values means more quality. filtering the result to get carrier and voice met points.

just comments needed, books and papers.

It seems that you are looking for books to learn about multiple topics. Skimming through your question I think these are the topics:

Aspects of signals (signal mean, its power,amplitude,frequency,etc )
Signals propagation in space or any other medium

I would recommend parts of the following books :

Bernard Sklar - Digital Communications Fundamentals and
Applications 2nd Edition
Bruce Carlson - Communications Systems 4th Edition
Simmon Haykin - Communication Systems  4th Edition
John G. Proakis - Digital Communications 4th Edition

A little bit more deep:

D.G.Manolakis - Statistical and Adaptive Signal Processing - Spectral Estimation, Signal Modeling, Adaptive Filtering and Array Processing
Fuqin Xiong - Digital Modulation Techniques 2nd edition

edition


MFCC delta: How to compute delta, in real time, with no future frames? I am trying to compute Delta coefficients for real-time audio analysis. As far as I know, a valid formula is:

$$d_t= \frac{\displaystyle\sum_{n=1}^{N}n\left(c_{t+n}-c_{t-n}\right)}{2\displaystyle\sum_{n=1}^{N}n^2}$$

It seems that I need $n$ future frames at each time step $t$ to compute both delta and delta-delta. However in real-time, I don't have future frames for several $\textrm{ms}$ (depending on my window segmentation). My methodology would be to compute, for timestep $t$, delta and delta-delta with an offset (i.e. compute the values from the previous frame or the one before). But then, these coefficients would correspond to one of the previous frames, bringing an offset to the feature vectors.

Is there any other way to get around this issue ? (maybe using a slightly different formula ?)

As long as you train and evaluate the same features then you will be perfectly fine with this delay.

In practice first two frames of the feature vector for both $\Delta$ and $\Delta\Delta$'s are zeros. So in fact the coefficients are delayed by two frames. These can be calculated according to the formulas:

$$\Delta[t] = c[t]-c[t-2]$$
$$\Delta\Delta[t] = c[t-2]-2c[t-1] + c[t]$$

Notice that the scaling factors are omitted intentionally. It's only a linear scaling of the feature space and should not change the results. For example in the case of GMM's the Mahalanobis distance takes care of differences in variance across dimensions.


How to use clustering file after training I trained the data using BOG with SVM , i trained it , now i just want to open it through xml , i don't want to train it again and again on every run ,
I trained the svm and now loading it using its loading function

CvSVM svm;
svm.load( "trainsvm.xml" );

But how to do the same with vocabulary file which i use for clustering using BOG algorithm , below is my code , which is training the data , but i want to open it now like svm

for (i = all_names.begin(); i != all_names.end(); ++i)
{
    Dir=( (count < files.size() ) ? YourImagesDirectory : YourImagesDirectory_2);

    Mat row_img = cv::imread( Dir +*i, 0 );

    detector.detect( row_img, keypoints);

    RetainBestKeypoints(keypoints, 20);

    extractor->compute( row_img, keypoints, descriptors_1);

    descriptors_1.reshape(1,1);

    bow.add(descriptors_1);

    ++count;
}

vector<string>::const_iterator k;

cout<<"CLUSTERING"<<endl;
Mat vocabulary = bow.cluster();

dextract.setVocabulary(vocabulary);

cv::Mat training_mat(num_img , dictionarySize,CV_32FC1);
cv::Mat labels(num_img,1,CV_32FC1);

This is clustering the data , now i don't want to run it again and again because it takes too much time , when it trained once , i just want to use it , so that i can speed up my program

You can use cv::FileStorage. You can follow this post here and this documentation here.


What is the Z-transform of $0.8^{n+2}u(n-1)$? I have 2 signals. One is $x(n)=(-0.5)^nu(n)$ and the other one is $y(n)=0.8^{n+2}u(n-1)$. I know that for the first one it is $X(z)= 1/(1+0.5z^{-1})$, but what about the other one? I know $y(n)$ is time shifted but i don't know how to find this z-transform. Any help would be appreciated!

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?


How to compare the syllable extracted from the audio I need to recognize just one bird. I am trying to compare the syllables extracted from the sound samples. As a result of the spectrogram function, the power spectral density is taken from the sound. They are 2 dimensional values. I like to compare this test signal with a known signal which is another 2 dimensional value. Could you please suggest what would be the best function to use for comparison. Please keep in mind that the syllable of the test bird could be time shifted or they closely match the original bird. I am attaching the sample spectrogram and its syllable extracted here after filtering..

Regards

A quick google for "bird voice recognition" yielded

https://peerj.com/articles/488/

which has an abstract describing that the authors refined the results they got using the old-school Mel-frequency cepstral coefficient comparison, as it has been common in (human) speech recognition for decades. I could start explaining what cepstrals and Mel-frequencies are, but there's so much literature out there that you literally just have to click read the first thing that hits you when you google for that. It's a relatively well-researched topic.


Proper generation of Complex white Gaussian noise using Matlab I am not comfortable in using the Matlab's awgn() function as it is hard to understand what is actually going on. So, I wanted to confirm if the following way is correct to generate noisy samples of a particular snr of 30 dB and a particular variance for data in complex domain.

Refer to the code below that generates some noise at a given SNR:

N = 100000;

% Generate some random signal
signal = randn(N, 1) + 1j*randn(N,1);

% Here, the signal power is 2 (1 for real and imaginary component)
signalPower_lin = 1/N*signal'*signal

% This corresponds to 3dB (assuming power=1 is 0dB)
signalPower_dB = 10*log10(signalPower_lin)

noisePower_dB = signalPower_dB - 30  % do 30dB SNR
noisePower_lin = 10^(noisePower_dB/10)

noise = sqrt(noisePower_lin/2) * (randn(N,1) + 1j*randn(N,1));

received = signal + noise;

The program outputs:

signalPower_lin =

   2.0042 + 0.0000i

signalPower_dB =

   3.0193 + 0.0000i

noisePower_dB =

 -26.9807 + 0.0000i

noisePower_lin =

   0.0020 + 0.0000i

Note that generating a complex noise of variance 1, you need to do

noise = sqrt(1/2) * (randn(N,1) + 1j*randn(N,1))

Since each component (real and imaginary) needs to have variance 1/2, such that their sum becomes 1.

To answer your points:
1) As a rule of thumb when to use 20 and when to use 10:
If you describe Powers or Energies, the factor is 10. If you describe amplitudes, the factor is 20. Then, in dB-scale both would have the same value, since:

$$
\begin{align}
P&=A^2 \text{ (Power = Amplitude squared)}\\
10\log(P) &= 10\log(A^2)\\
&=20\log(A)
\end{align}
$$

So, anytime you use a values which creates something meaningful when squared (like amplitude), use 20. Everytime, you have something that is meaningful when taken the square root from (like Power), use 10.

Note that in my code, I use noise = sqrt(noise_power/2) * randn(...) to generate the noise, i.e. I transform the noise power into the noise amplitude prior to multiplying the normal random variable (which corresponds to changing the factor 10 to 20 in the exponent). Also remember the rule
$$Var(aX) = a^2Var(X)$$
which is exactly the reason why you need to take the square root of the power, such that your resulting noise has the required power.

2) Here, I just transform noisePower_db into noisePower_linear, which is by definition linear=10^(db/10). Note that $10^{-SNR/10}$ is similar, since SNR=Signal/Noise, i.e. the noise is in the denominator (consider $10^{-1/x}=10^x$.


Units of FFT frequency I'm having some trouble understanding FFT plots and measures.

According to this, this will give me the magnitude plot of FFT, but I want the area under from 1 to 10 Hz for example.

Is this next formula I made right?

My doubt is, what's the relation between Hz and rad/s, f and w?

The formula you made will work if what you are after is simply "the area under the the magnitude plot from $\omega_i$ to $\omega_f$."

In many cases you can think of Hz as meaning "rotations per second."  In that case, it becomes more obvious that 1 rotation = 2$\pi$ radians.  Since $f$ has units of Hz and $\omega$ has units of radians, the conversion is simply:
$$
1 \, \text{Hz} = 2 \pi \, \frac{\text{rad}}{\text{s}}
$$
So evaluate your integration from $\frac{1}{2\pi}$ to $\frac{10}{2\pi}$ and you should have the result you are after.


FIR butterworth also possible or just IIR? can we use butterworth and other such filters like chebyshev,elliptic etc  with FIR

or they can be only used with IIR?

Filters according to those optimality criteria only exist as IIR filters. They are derived from the corresponding analog prototype filters via the bilinear transform, and this naturally results in IIR filters, i.e., filters with zeros and poles (away from the unit circle).

Of course, you can approximate these filters by FIR filters. The most straightforward way would be to truncate the infinite impulse response using some window function. But such an FIR filter is not optimal anymore since it can only approximate the optimal response.


Biorthogonal Wavelet Expansion Consider a multiresolution analysis (MRA) $\{V_j\}$  in $L^2(\mathbb R)$ generated by the scaling function $\phi(x)$ and such that for each $j\in\mathbb Z$, $V_j\subset V_{j-1}$. For a given function $f$ I will use the notation $f_{j,k}(x) = 2^{-j/2}f(2^{-j}x-k)$.

Let $\psi(x)$ be the wavelet corresponding to $\phi(x)$, then $\{\psi_{j,k}\}_{k\in\mathbb Z}$ spans the orthogonal complement $W_j$ to $V_j$ in $V_{j-1}$, namely, $W_j = \overline{span}\{\psi_{j,k}\}_{k\in\mathbb Z}$ is such that $V_{j-1} = V_j\oplus W_j$ and $V_j\perp W_j$. Thus for a fixed $J<0$ one can define the projection of a function $f\in L^2(\mathbb R)$ into $V_j$ through the expansion
$$
P_J f(x) = \sum_{k\in\mathbb Z} \langle f,\phi_{0,k}\rangle\phi_{0,k}(x) + \sum_{j=0}^{J+1}\sum_{k\in\mathbb Z}\langle f,\psi_{j,k}\rangle \psi_{j,k}(x)\qquad (1)
$$
where the equality is in the $L^2$-sense. That makes sense to me due to the orthogonality of $V_0$ and $W_0\oplus W_1\oplus ... \oplus W_{J+1}$.

Consider now a biorthogonal MRA. In this case one consider two dual MRAs $\{V_j\}$ and $\{\tilde V_j\}$ generated by the scaling functions $\phi(x)$ and $\tilde \phi(x)$. In a similar way as in the orthonormal case, one defines the dual wavelets $\psi(x)$ and $\tilde \psi(x)$ which span repsectively the (non orthogonal) complements $W_j$ of $V_j$ in $V_{j-1}$ and $\tilde W_j$ of $\tilde V_j$ in $\tilde V_{j-1}$. Namely one has $V_{j-1}=V_j+W_j$ and $\tilde V_{j-1} = \tilde V_j+\tilde W_j$ with $W_j = \overline{span}\{\psi_{j,k}(x)\}_{k\in\mathbb Z}$ and $\tilde W_j = \overline{span}\{\tilde \psi_{j,k}(x)\}_{k\in\mathbb Z}$.

In this case, for fixed $J$, one can write
$$
P_J f(x) = \sum_{k}\langle f, \tilde \phi_{J,k}\rangle\phi_{J,k}(x)
$$
and an analogous expansion exists for $\tilde P_j$. My question is: does an expression of the kind of (1) hold also in the biorthogonal case? i.e. can I write
$$
P_J f(x) = \sum_{k\in\mathbb Z} \langle f,\tilde\phi_{0,k}\rangle\phi_{0,k}(x) + \sum_{j=0}^{J+1}\sum_{k\in\mathbb Z}\langle f,\tilde\psi_{j,k}\rangle \psi_{j,k}(x)\;?
$$
It is not clear to me that this is possible since it is not true that $V_0$ is orthogonal to $W_{0}+...+W_{J+1}$. However the following argument appears to hold:

Let $Q_j := P_{j-1} - P_{j}$ (the detail operator) and define the dilation operator as $D_j f(x) := 2^{-j/2}f(2^{-j}x)$. We have $$P_{j+m}f(x) = D_{j} P_m (D_{-j} f(x))$$  and hence
\begin{align*}
Q_jf(x) &= P_{j-1}f(x) - P_{j}f(x) = D_j (P_{-1} - P_0) \big(D_{-j} f(x)\big)\\&= D_j Q_0\big(D_{-j}f(x)\big)
\end{align*}
Moreover it can be shown that $Q_0 f(x) \in W_0$ and from the above relation one obtain, for each $j\in\mathbb Z$
$$
Q_j f(x) = D_j\sum_k \langle D_{-j}f,\tilde \psi_{0,k}\rangle\psi_{0,k}(x) = \sum_{k} \langle f,\tilde \psi_{j,k}\rangle \psi_{j,k}(x)
$$
Thus by definition of $Q_j$
\begin{align*}
P_{j}f(x) &= P_{j+1}f(x) + Q_{j+1} f(x) = P_{j+1}f(x) + \big(P_{j+2}f(x)-P_{j+2} f(x)\big)+Q_{j+1} f(x)\\ &= P_{j+2}f(x) + Q_{j+2} f(x) + Q_{j+1}f(x) \\&\quad\vdots\\&=P_0 f(x) + \sum_{i=0}^{j+1} Q_i f(x)
\end{align*}
and thus (1) holds.

Is this correct?

Let me try to give you a short introduction to the concept of bi-orthogonal bases from a mathematical point of view. You are already using the concept, but your notation and nomenclature is obfuscating what is really going on. I will make this answer independent of the wavelet context you have, so that it is more general. You will certainly be able to apply it to your situation however.

You have a non-orthogonal basis, i.e. a set of linearly independent vectors spanning your $L^2(K)$ Hilbert space $\mathcal{H}$. That vector space induces another vector space, that of all linear forms $F:\mathcal{H}\to K$. We can make this space a dual space of $\mathcal{H}$ if we take the Hilbert space inner product $\langle \cdot, \cdot\rangle : \mathcal{H} \times \mathcal{H} \to K $ to generate a linear form $\langle \phi, \cdot \rangle : \mathcal{H} \to K$ for every $\phi \in \mathcal{H}$.

In an effort to simplify the notation a bit and make things more intuitive, we can introduce Dirac's Bra-Ket notation for vectors and their duals. Instead of $\phi \in \mathcal{H}$, we write the Ket $\left| \phi \right\rangle$ and instead of $\langle \phi , \cdot \rangle \; \phi \in \mathcal{H}$ we write the Bra $\left\langle \phi \right|$. Applying the linear form $\left\langle \psi\right|$ to $\left| \phi \right\rangle$ is then simply the notational fusion to a braket standing for the inner product $\left\langle \psi | \phi \right\rangle$.

In the bi-orthogonal basis business, this distinction is often ignored, making everything more confusing that it should be. Consider a basis that is not orthogonal with respect to the canonical inner product. We can of course use the usual knowledge of linear algebra to still construct everything of interest and expand vectors with respect to this basis, but it tends to get messy quickly.

That's where a neat trick comes in. We simply invent a new inner product, i.e. a non-degenerate, non-negative bilinear form $B: \mathcal{H} \times \mathcal{H} \to  K$, so that $B( \phi_n, \phi_m ) = \delta_{n,m}$ for your basis $\{ \phi_n, n\in\mathbb{N}\}$ and the Kronecker symbol $\delta$. Because the basis is complete, $B$ is entirely determined by this condition. We can also write this bilinear form using a linear operator $B:\mathcal{H}\to\mathcal{H}$ with the same name. Then we get $$\left\langle \phi_n | B | \phi_m \right\rangle := \left\langle \phi_n | B \phi_m \right\rangle = B( \phi_n, \phi_m ) = \delta_{n,m}$$.

At this point, it makes sense to define the adjoint $A^\dagger$ of a linear operator $A$ like follows: $$\left\langle A^\dagger \phi_n | \phi_m \right\rangle = \left\langle \phi_n | A \phi_m \right\rangle$$

Of course it follows from the definitions above: $$\left\langle B^\dagger \phi_n \right| = B(\phi_n,\cdot)$$
and therefore $$\left\langle B^\dagger \phi_n | \phi_m \right\rangle=\delta_{n,m}$$
which means $\left\langle B^\dagger \phi_n \right|$ is the dual of $\left|\phi_n\right\rangle$ under $B$, whereas the dual under the canonical inner product is still $\left\langle \phi_n \right|$.

Now what you can often see, and I cannot stress enough how much I think this is a bad idea, is that the dual under $B$ is mapped back to $\mathcal{H}$ as a dual under the canonical inner product. That means you construct a second basis $\left| \tilde{\phi}_n\right\rangle =\left| B \phi_n \right\rangle$ which is bi-orthogonal to the first basis with respect to the canonical inner product:
$$ \left\langle \tilde{\phi}_n | \phi_m \right\rangle = \delta_{n,m}$$
This obfuscates the real meaning of the second basis and puts both on the same footing and suggests that both bases live in the same space and transform in the same way, which they don't.

So, what is the alternative? Not taking the $B$-dual back to $\mathcal{H}$ makes things much more intuitive. You want to expand $\left|\psi\right\rangle$ in the basis $\phi_n$? Just use the $B$-dual: $\left|\psi\right\rangle=\sum_n \left| \phi_n \right\rangle \left\langle B^\dagger \phi_n | \psi \right\rangle $. Or you want to project onto a vector $\left| \kappa \right\rangle$? Just use the $B$-dual: $\left| \kappa \right\rangle \left\langle B^\dagger \kappa | \psi \right\rangle$. So whatever you do, just use $B$ as your inner product and you can work like with an orthogonal basis.

And to come back to your original question. You say that $V_0$ is not orthogonal to $\oplus_n W_n$. Yes, but only with respect to the canonical inner product. With respect to $B$, they are orthogonal. And so your argument for the orthogonal basis case works.

I hope all this does not confuse you. I've tried to be as clear as possible, but it's not a trivial topic, so don't be frustrated if you have to read the answer several times until it makes sense. If I can clarify anything, please let me know.


pitch extraction for speech I am building an LPC analysis tool, and am at the point of needing to do pitch period analysis.  I was initially interested in using the Gold & Rabiner algorithm for this, and got a hold of a FORTRAN implementation.  Before I start converting that implementation to a modern programming language, I was curious if anyone here had any suggestions for other pitch detection algorithms that are maybe a bit more modern and robust...  And even better, if there are any open source libraries that have apis for pitch period analysis that my tool can hook into?

EDIT:  I am building an LPC analysis tool for speech waveforms.

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$.


Adaptive Band pass filter for extracting respiratory waveform from accelerometer data I am recording data from an accelerometer attached to the chest (1000Hz). I need to extract the respiratory waveform. I tried an adaptive bandpass filter based on a dominant frequency in my signal based on https://cutt.ly/j24zdj.

Steps in brief: Take frequency spectra, find dominant freq f0, Make a bandpass filter as [max(0.1,f0-0,4), 0.4+rm]

The respiration rate could be between anything from 0.1 to 2 Hz. For deep and normal breathing, the filter works just fine. But for rapid and shallow breathings (low amplitude but high freq), it fails. It identifies the dominant freq incorrectly and thus the entire filter goes wrong.

For example in the image below: The dominant freq should have been around 0.6, but is identifies as 1.7 (Slight diff in amplitude). So filter works badly.

1) How do I design a filter for this? Any suggestions. I can't control the noise in this freq range. The person is as still as possible.

2) The red line is spectra after filtering. Why has the amplitude of certain frequencies increased? I am using the 'spec' function in R: uses the hanning window for computing.

It's not the issue with plotting. The amplitudes have actually increased.

Under the extremely controlled conditions required to measure respiratory rate using this method, you can try to use the autocorrelation (which is more robust in the presence of noise) to determine the "dominant" frequency in your signal.

  2) The red line is spectra after filtering. Why has the amplitude of certain frequencies increased? I am using the 'spec' function in R: uses the hanning window for computing.

If that is not a mere rescaling of the waveform as it gets plotted over the same figure (?), then it would be safe to assume that the filter included some overall gain.

Hope this helps.


frequency/phase offset estimation in satellite DOA beamforming applications In satellite QPSK demodulation, we first estimate and remove the phase offset introduced by the channel. Then, we demodulate the baseband information bits (I & Q channels).

In satellite DOA beamforming, Channel can also introduce phase / frequency offset. When we detect the desired signal from a certain direction in the frequency domain (spectrum), how could we detect this channel introduced phase / frequency offset?

Typically, you distribute a mixing oscillator to all ground stations, so that their respective phase is known and can hence be corrected upfront. Having known phase between these receiver parts typically requires a lot of RF design effort and external phase calibration.


Estimate the Filter Coefficients of 1D Filtration (Convolution) I have an output signal $y$ which is an input signal $x$ convolved $\star$ with an impulse response function $h$ with some added noise $n$ :

$$y(t) = h(t) \star x(t) + n(t)$$

I know the input signal $x$ and output signal $y$ and would like to calculate $h$ the impulse response function. I found that deconvolution is not as straightforward as convolution because the input signal contains zeros and then division in the frequency domain would not be defined. Looking around the internet for ways to "deconvolve" if found two methods: Wiener deconvolution and regularized deconvolution. The Wiener deconvolution seemed easier to understand so I wanted to try and implement it in Matlab (the Matlab function deconv gives me errors about the input signal having a zero at the first entry and if I read the help file it only seems to work correctly for polynomials?).

So per the wikipedia explanation you want to find $g$ so that:

$$\hat{x}(t) = g(t) \star y(t)$$

But then in the definition of how to calculate $G$ they use all variables in the original equation. Also they only show how to find $x$ but probably $x$ and $h$ can be exchanged because convolution is commutative, but I'm unsure about the correct length of both vectors. Currently they are the same length.

$$G(f) = \frac{H^*(f) S(f)}{|H(f)|^2 S(f) + N(f)}$$

where:

$H = {\tt fft}(h)$
$G = {\tt fft}(g)$
$S =$ power spectral density of $x$ ? is this ${\tt fft}(x)$?
$N =  $mean power spectral density of $n$, don't really understand what this is

My question is how do I get the impulse response function $h$ without already knowing it (as it is both in the definition of $G$ and in the original equation)? Since I know both input and output it should not be very different from finding the input with a known impulse response function. I want to do this in Matlab.

This is a nice question.
I will try solving it using 2 approaches (Which are basically the same).

The solution is the Least Squares Solution:

$$ \hat{h} = \arg \min_{h} \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$

We assume data is given in finite discrete form (As it is in practice).
The convolution is done in valid mode (Like in MATLAB valid property).

Direct Solution

One could write the above as:

$$ \hat{h} = \arg \min_{h} \frac{1}{2} \left\| X h - y \right\|_{2}^{2}  $$

Basically using the commutative property of the convolution one could the above as $ x $ was the filter and hence build a Convolution Matrix from $ x $.

Now, the solution is the usual Least Squares:

$$ \hat{h} = \arg \min_{h} \frac{1}{2} \left\| X h - y \right\|_{2}^{2} = \left( {X}^{T} X \right)^{-1} {X}^{T} y $$

Pay attention that for long signals the matrix is huge.
Hence this method is practical only for small signals.

Iterative Solution

Well, one could use the Gradient Descent Method to minimize the cost function:

$$ f \left( h \right) = \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$

The Gradient is given by (Correlation is the Adjoint Operator of Convolution):

$$ \frac{d}{d h} f \left( h \right) = \left( h \ast x - y \right) \star h $$

Where $ \star $ stands for Cross Correlation.

In MATLAB Code it is given by:

hObjFun = @(vH) 0.5 * mean( (conv2(vX, vH, 'valid') - vY) .^ 2 );
vG      = conv2(flip(vX, 1), (conv2(vX, vHEst, 'valid') - vY), 'valid');

Where all data is in Column Vector (The conv2() is used as it is faster since conv() is just a wrapper around it).

This method is fast and not limited by the size (Or less, to be accurate).

The full MATLAB Code is in my StackExchange Signal Processing Q18993 GitHub Repository.


How to interpret these FFT results? I have recently started to use MathNet.Numerics, and one of my needs is to perform FFT analysis of Electromyographical signals.

The first test I did was to get some 15 seconds of signals at 2kHz, this is how it plots:

I have a Python script I wrote some time ago, and the FFT spectrum generated is below:

lenfft = len(signal)/2.
nyq = 2000/2.
freq_axis = numpy.linspace(0, nyq, lenfft+1)
spectrum = numpy.abs(numpy.fft.rfft(signal))
plt.plot(freq_axis, spectrum)
plt.show()

The problem is, when I use MathNet.Numerics, I get a weird result (code and image below):

Complex[] signal = ArrayOfDoublesFromFile()
                        .Select(val => new Complex(val, 0))
                        .ToArray();

Fourier.Forward(signal, FourierOptions.Default); // inplace

SaveToFile(signal);

I know that if the inputs are pure real numbers, the output should be symmetrical or something like that, but since Fourier.Forward() requires a MatNet.Numerics.Complex[] array, I am creating an array with imaginary part set to zero, and also using only real part of result for plotting, is that correct?

So, the question is: whats the meaning of the output generated by MathNet.Numerics.Fourier.Forward(), and what should I make if I want a result in the visual form most commonly associated with FFT, that is, the middle figure?

Note the symmetry in your 3rd plot. That tells me the MathNet.Numerics FFT is computing a full-frequency FFT (both positive and negative frequencies), and plotting the FFT results where zero Hz should be the center of the horizontal axis. However, Problem# 1, for some reason the labeling of the horizontal axis does NOT have zero Hz in the center. I can’t explain that. Problem# 2, there’s something very wrong with the ‘flat-top’ spectral shape in your 3rd plot. Problem# 3, it seems impossible to me that your 3rd plot should have spectral components whose magnitudes are 10 raised to the 15th power. (Looking at your 1st plot I’d expect your spectral magnitude results to be on the order of 10 raised to the 7th power.) Stated in technical terms, your 3rd plot is “truly FUBAR.”

The answer, heltonbiker, is for you to determine precisely, in great detail, EXACTLY what processing is taking place in your MathNet.Numerics FFT and plotting routines. Good Luck.


How can I extract information about timbre from a .wav file. In which form can I study those features? To study the frequencies of a piece of audio, we can generate a spectrogram from that audio file. I'm wondering how the features of timbre are saved in such tiny wave files, and how we can extract those features and study the timbre? Thus we know what instruments are played in the audio?

I think the time history of the wave file is enough to obtain all sorts of features. There are mathematical formulae to extract the features from a given time series.

Some important features are:

Time Domain Audio features

(a) Short-term Energy, (b) Zero Crossing Rate, (c) Entropy of Energy, (d) Basic Statistics (short term mean and standard deviations).

Frequency Domain Audio Features

(a) Spectral Centroid and Spread, (b) Spectral Entropy, (c) Spectral Flux, (d) Spectral Rolloff, (e) Mel-Frequency Cepstrum Coefficients, (f) Chroma Vector

More details regarding audio analysis can be found in this book
https://www.elsevier.com/books/introduction-to-audio-analysis/giannakopoulos/978-0-08-099388-1

If you have enough data, using these features you can build a machine learning model to correctly classify an unknown sound.


Why is a precoder necessary for DQPSK and what does it accomplish? I've implemented a soft-decoder for DQPSK using the wonderful answers I received here:

How to soft decode DQPSK?

To get the soft-decoder working properly I needed to precode the data I was sending out.  I implemented the precoder mentioned in this paper:

$I_k=\overline{u_k \oplus v_k}*(u_k \oplus I_{k-1})+(u_k \oplus v_k)*(v_k \oplus Q_{k-1})$
$Q_k=\overline{u_k \oplus v_k}*(v_k \oplus Q_{k-1})+(u_k \oplus v_k)*(u_k \oplus I_{k-1})$

I'd like to know why this precoder is necessary -- what does that complicated expression of XORs actually accomplish?

Here's a table showing what the equation yields.  If "to_encode" is 00, the to_send symbol is the same as the previous ("prev") symbol.  If the "to_encode" is 11, the to_send symbol is the previous symbol xor 11.  What is the meaning in other cases?

to_encode prev  to send
[ 0 0 ] [ 0 0 ] [ 0 0 ]
[ 0 1 ] [ 0 0 ] [ 1 0 ]
[ 1 0 ] [ 0 0 ] [ 0 1 ]
[ 1 1 ] [ 0 0 ] [ 1 1 ]
[ 0 0 ] [ 0 1 ] [ 0 1 ]
[ 0 1 ] [ 0 1 ] [ 0 0 ]
[ 1 0 ] [ 0 1 ] [ 1 1 ]
[ 1 1 ] [ 0 1 ] [ 1 0 ]
[ 0 0 ] [ 1 0 ] [ 1 0 ]
[ 0 1 ] [ 1 0 ] [ 1 1 ]
[ 1 0 ] [ 1 0 ] [ 0 0 ]
[ 1 1 ] [ 1 0 ] [ 0 1 ]
[ 0 0 ] [ 1 1 ] [ 1 1 ]
[ 0 1 ] [ 1 1 ] [ 0 1 ]
[ 1 0 ] [ 1 1 ] [ 1 0 ]
[ 1 1 ] [ 1 1 ] [ 0 0 ]

Regard the dibits as Gray code representations of the
integers $0,1,2,3$, more specifically,

$$[0, 0] \leftrightarrow 0, ~~ [0, 1] \leftrightarrow 1, ~~ [1, 1] \leftrightarrow 2, ~~ [1, 0] \leftrightarrow 3.$$

Then, the precoding scheme is nothing but differential encoding for QPSK
with

to send = prev - to_encode modulo $4$.

For example the line [ 1 0 ] [ 0 1 ] [ 1 1 ] in the question says that to encode $3 = $ [1 0]  when the previous symbol was $1 = $ [0 1], we send $1-3 = -2 \equiv 2 \bmod 4$ where $2 =$ [1 1].
As stated in the answers and comments
on the previous question,

  Differential encoding for QPSK is not the same as differential
  (binary) encoding of the I and Q bit streams separately.

So the precoding scheme is just differential encoding for QPSK.
As to why it is included in the system, it must be that the
channel characteristics and system requirements are such that
coherent detection of the signal is not thought to be feasible,
and so the receiver is being implemented
is a differentially coherent detector for DQPSK (differentially encoded
QPSK).  Note that coherent detection requires that the receiver
phase lock loop acquire (and track) the incoming carrier in
both frequency and phase (including resolution of
the four-fold phase ambiguity) while for differentially coherent
detection, frequency acquisition and tracking suffices, so that
the receiver is simpler and more robust, though of course
for a given SNR, the BER of the differentially coherent
receiver is higher than the BER of the coherent receiver.


LTI systems inverse of each other So my question looks like this:

  Suppose that we have two LTI systems with impulse responses
  $$h_1(t) = \frac 12\delta(t-1)\quad\text{and}\quad h_2(t)=2\delta(t+1).$$
  Determine whether these systems are inverses of one another.

Here's what I've done but I'm stuck. I'd appreciate some help with this

$$h_1(t)\star h_2(t)\overset{?}{=}\delta(t)$$
\begin{align}
\int_{-\infty}^\infty h_1(\tau)h_2(t-\tau)d\tau &=\int_{-\infty}^\infty \frac12\delta(\tau-1)2\delta(t-\tau+1)d\tau\\
&=\int_{-\infty}^\infty \delta(\tau-1)\delta(t-\tau+1)d\tau\\
\end{align}

$$u = \tau - 1\implies \frac{du}{d\tau}=1\implies du = d\tau$$

$$\implies \int_{-\infty}^\infty h_1(\tau)h_2(t-\tau)d\tau = \int_{-\infty}^\infty \delta(u)\delta(t-u)d\tau$$

So you have
$$
\int_{-\infty}^{+\infty} \delta(u) \delta(t-u) du
$$

Try using the sifting property of $\delta(t)$:

$$
\int_{-\infty}^{+\infty} f(u) \delta(T-u) du = f(T)
$$


Multiple different length Windows for high resolution realtime FFT I am looking at trying to achieve a 'realtime' (as quick as possible to acquire, process and present data) FFT application to analyse audio.

I have setup my app to acquire a set of samples, apply a 50% overlap window function, zero pad, compute FFT, extract magnitude and phase and display these.  Currently I can do this with no issues, and minimal latency to the incoming signal.

However, as I have chosen a small sampling window (2048 samples) to achieve a fast time, I realise that this hinders resolution.  I am also aware that there is only so much zero padding can achieve, if the acquired data is not there to begin with.

So my question is to do with how to achieve a better resolution for lower frequencies.  It is my understanding that for accuracy, the time window should contain at least one cycle of the lowest frequency.

Is there a way that I can perform multiple FFTs on a signal with varying lengths and 'glue' the resultant FFT together. ie, a longer FFT will do the lower frequency range, and a shorter for the higher range.  And if I can do this, will this retain phase information, or will I lose this completely?

You can, and you don't need to zero pad, but a little overlapping does come into play.

Skipping the gnarly details and hitting this with a conceptual approach.  Applies to two or more "levels".

Start with a longer frame for your low tones.  The key is you need to smooth your signal fairly heavily.  Lots of ways to do this.  Then do a DFT (what a FFT really is) on the long interval.  The reason for the smoothing is to eliminate alias frequencies at this level.  The effects of your smoothing will be frequency dependent, so there are formulas for the different smoothing techniques.  Do the inverse of the FFT of the smoothed signal with attenuation compensation applied.  Subtract this from your original signal.  You have now documented and removed your lower frequencies.  Divide the interval into smaller chunks, and do DFTs as normal at the top level.

"Pasting" them together can also be done in different ways.  Easiest is just to display one spectogram above the other as "display bands".  You can also rescale the output from being linear in frequency to exponential if you want to put your results on an octave scale.


Music/Speech Classification: How are frame based features used in a feature vector? I'm a little confused on how frame based high dimensional features such as mfcc's and bfcc's are applied in the classification of musical instruments.

Are statistical measures of the features used e.g. using the standard deviation across the frames used as a feature? Or is each frame of coefficients used as independent features e.g. 2 frames, 13 coefficients --> 26 features? a combination of both or something completely different?

Thanks

Let us say you have a set of songs you want to classify by genre.

There are two approaches if the classification algorithm you use works on unidimensional vectors (eg: neural networks, SVM...):

Either you have one feature vector per song; and your feature vector contains statistics (mean, standard deviation, etc.) summarizing all the MFCCs computed at the frame level. In this case, your training set contains one vector per song; and during classification, the classification algorithm is run once per song.
Or you have one feature vector per frame. In this case, your training set is made of frames; and the classification algorithm is run once per frame. Voting is used to select the label for the whole songs. It is common to weight the votes to give more important to the middle of the signal (for example, for musical genre classification, the song might start with a slow/soft intro which might bias the classification).

A particularity of both solutions (the so-called "bag of frames" representation) is that they are invariant to permutations in the sequence of MFCC vectors - they do not take into account the temporal structure of the signal. Another approach is to use machine learning algorithms or models (such as HMMs) which can intrinsically deal with time-series - in which case your representation is the sequence of feature vectors.

The approach you describe in which MFCC vectors for each frame are concatenated to yield a gigantic feature vector is not practical: the large feature vectors would make the learning task computationally extensive, it would yield feature vectors of different dimensions for signals of different durations; and two signals identical modulo a delay would have very different feature vectors.


coded filter vs fdatool : inconsistency in the results for same parameters I tried creating a filter in 2 different ways in Matlab and I obtain different results: one filter is stable, whereas the other one is unstable.

Coded filter

delta_t = 1.53846e-04;
Fs = 1/delta_t;

cut_F = 8;
Wn = cut_F/(Fs/2);
ftype = 'high';
[b,a] = butter(8,Wn,ftype);
fvtool(b,a)

The result is an unstable filter

fdatool

This is what I obtain with the toolbox using the same parameters:

I am using already implemented functions in both cases...Why I get so different results?

Note that the fdatool gives you a structure with cascaded second-order sections, whereas the fvtool uses a direct-form II transposed structure. In combination with the filter specs this results in the difference you see, because the very low cut-off frequency combined with the relatively high filter order results in a cluster of poles and zeros very close to the DC point in the complex $z$-plane (i.e., close to $z=1$). Computing the exact locations of the poles and zeros is an ill-conditioned problem, and in the case of the direct-form II transposed structure, you need to find the roots of polynomials of degree $8$. Due to the round-off errors in the computation of the roots, some of the poles appear to lie just outside the unit circle, hence the conclusion that the filter is unstable. For the second order sections, the poles can be computed more accurately and they all turn out to lie inside the unit circle, hence the conclusion that the filter is stable.

So, is the filter now stable or unstable? Well, in theory it is stable. In practice, however, it is probably of no use because its poles are virtually on the unit circle. The tool probably offers you to compute the exact pole locations of the second order sections. You'll see that most of the pole radii are very close to $1$.


Trouble with understanding moving average filter as an oversampling filter I am struggling with understanding how an moving average filter can increase ENOB in a given sample.

I wrote a test for a simple MA (moving average) filter in matlab. The algorithm to create new MA value simply subtracts the new signal value from and old one (that is $n$ samples away), and adds the previous value of the filter to it. I run the fft and it seems to work as low pass filter  - exactly how expected.

I am aware that this same algorithm can be used to increase number of bits of resolution for a sample. Furthermore, if the old value that we used to build the difference with the new value is $n$ sample away, then the filter value is supposed to have the resolution increased by $n$ bits....
I just don't really understand how is that possible. To my understanding:

we take the old filter value ($m$ bit long),
we add a new value to it that is a difference of the current signal value and a past signal value, that is $n$ samples away - both values also $m$ bit long,
how is it possible that resulting filter value could be now $m+n$ bit long?

Lets say we have 12 bit($m$) signal values. How is it possible to 'gain' additonal 6 bits($n$) bits by adding the signal value that is 6 samples away from the current sample? Why should the outcome have the resolution of 18 bits?

Thanks in advance for your answers!

A moving average over 6 samples will not add 6 bits of precision to the ENOB (equivalent number of bits) as stated in the question. In general, a moving average over $N$ samples will increase resolution by $log2(N)/2$, and for 6 samples this would be approximately 1.3 bits.

This approximation holds to the extent that the quantization noise on each sample is represented by a uniformly distributed white noise process (which holds well when there are a sufficient number of quantization levels and the sampling clock is not correlated to the signal being sampled). From statistics, when you add N independent random variables, the standard deviation increases at $\sqrt{N}$. Thus in the average where we also divide that sum by $N$ the standard deviation is scaled by $\sqrt{N}/N = 1/\sqrt{N}$.

This results in a 3 dB reduction in quantization noise for every doubling of samples in the average (from $20Log_{10}(1/N)$ since dB is $20Log_{10}(\sigma)$ where $\sigma$ is the standard deviation. As already detailed in the post linked below, each bit of precision represents a 6 dB reduction in quantization noise:

What are advantages of having higher sampling rate of a signal?

A moving average over 4 samples would add 1 bit of precision, 16 samples would add 2 bits, etc.

To support the increased number of bits, the result from the addition of the samples must have that many number of bits which is also apparent to eliminate overflow conditions. Filtering a digital signal (which a moving average is) grows the signal! This is why we must use extended precision accumulators when filtering. We can scale the accumulated signal back after summing (for the average) while scaling the signal first prior to the summation, or equivalently scaling the coefficients, would result in further noise contribution.


Gaussian Blur Detection Let's say you have an original image and a version of the same image that may have been convoluted with a Gaussian blur. How could you demonstrate that the Gaussian blur has been applied and calculate the blur radius?

If you have the blurred and unblurred images, what I think you are asking is how you can recover the point spread function (PSF).

In the absence of noise, this is theoretically possible by considering the operation in the frequency domain. The blurred was introduced by a convolution of the unblurred image with a Gaussian kernel. In Fourier space, convolution is applied by multiplication of the image with the kernel, so if you have two of the operands (the source and the result), it should be possible to recover the blur kernel by division of the respective Fourier transforms.

This is effectively what some deconvolution operations attempt to do when the PSF is unknown or is a mixture of defocus and motion blur (e.g. blind deconvolution, Wiener filter, etc.).

So, in the frequency domain, the blur can be expressed as

G(u,v) = H(u,v)F(u,v) + N(u,v)

Where G is the Fourier transform of the blurred image, H is the blurring function, F is the source image and N is additive noise. Assuming the noise is negligible (unrealistic in most real world situations), the blur kernel can be recovered by calculating

H(u,v) = G(u,v) / F(u,v)

and then taking the inverse Fourier transform of H to obtain the PSF as an intensity image. From this you could measure the radius.


Fequency Domain Oversampling - decimating when going back to time domain I'm currently trying to simulate Frequency Domain Equalization, and the effect that over sampling in the frequency domain has on the affect of the equalization. If I understand things correctly this is the process that's supposed to happen using a single Carrier when transmitting

Modulate and transmit the symbols to create $x(n)$, each symbol being $N$ samples.
add noise, $z(n)$, and channel effects, $h(n)$, for the receiver, $ y(n) = h(n)\star x(n) + z(n)$
Zero pad, where $ZP$ is the number of added zeros, in the time domain to over sample in the frequency domain. $FFT(y(n),N+ZP)$. Note that increasing the order of the FFT should be the same as zero padding in the time domain
Do Equalization in the frequency domain
Go back to the time domain
Demodulate the symbols

The part I'm not sure about is step 5, where we go back to the time domain. In order for my number of samples per symbol to be correct, I need N samples, but I have $N+ZP$ samples if I take the same size $IFFT$. Doing the Zero padding and taking the $FFT$ should be interpolation... but i'm not sure how to do the decimation so that I have the correct number of samples. Does anyone have any idea about this? Thanks!

Your time domain waveform after 6 will be zero-padded as in step 3 with the additional "smearing" due to your equalization to remove the channel effects that will go into the zero padded area. The entire waveform in time will be delayed by the group delay of your frequency domain equalization, which may not necessarily match the group delay of the channel distortion.

If you applied no equalization (by multiplying in frequency by your frequency domain equalization); the IFFT would recover your time domain waveform in its original form with the noise and channel distortion with the additional zeros at the end that were introduced by the zero-padding when taking the FFT. (Along with boundary conditions at the transition such as Gibbs phenomena due to the effective rectangular window applied). Multiplying by your equalizer in frequency is identical to convolving in time the original time domain waveform with the Fourier transform of your equalizer fucntion, and therefore will delay this waveform based on the group delay of this filter (where group delay is $\frac{d\phi}{d\omega}$). Equalizers are typically not linear phase (unless your channel distortion was which is unlikely) so the group delay is not constant versus frequency; but the point is the resultant delay of the recovered original waveform will be completely dependent on the net delay of the channel and equalizer. With knowledge of the channel and specific equalizer function, the delay can be determined.


Wind Turbine Vibration Analysis: Which approach? I'm analysing vibration data from 4 wind turbines (WTs). 8 different sensors are sampled at 25.6 kHz for 10 seconds once a day. I have data from around 400 days (intervals). The plot underneath is for one of the four WT in the time domain.

The frequency domain from FFT of the four turbines (Gearbox HSS) in the the range 0-12kHz:

Current plan and questions:

I want to look for faults in the gearbox, and I am suspecting that WT 4 (bottom right plot above) is most degraded.
I want to apply some kind of high frequency technique for looking for early stage bearing damages in the higher frequency range, since these are not captured by FFT.

What approaches should can I consider here?  (I have looked a bit at Hilbert transform, but I'm not exactly sure how that would help me.)
I am also considering performing some kind of enveloping/wavelet transform.
Should I be looking at lower frequencies as well?

This is not my main study field, and I want to try building a machine learning classifier with variables obtained from the signal processing.

Plots of the lower frequency (0-2000Hz) development over the 400 intervals for all 4 turbines:

Thanks!

One technique that I've seen used in determining bearing faults is using the kurtosis of the vibration signal.

You can track as a function of time what Wikipedia calls the sample excess kurtosis. This is the kurtosis that is different from the kurtosis you would see if the signal was Gaussian distributed.

The sample excess kurtosis is defined as:

$$
\frac{m_4}{m_2^2} - 3 = \frac{\frac{1}{N} \displaystyle \sum_{n=0}^{N=1} (x_n - \bar{x})^4}{\left(\frac{1}{N} \displaystyle \sum_{n=0}^{N=1} (x_n - \bar{x})^2\right)^2} - 3
$$

where $\bar{x}$ is the sample mean, $m_2$ is the sample second order moment, and $m_4$ is the sample fourth order moment.

Brüel & Kjær have a nice writeup about using kurtosis, though they use the actual kurtosis rather than the excess kurtosis.

Why?

Bearing faults tend to show up when grit or dirt or broken bearing pieces get into the rotating elements.  When that happens, the vibration signals tend to be more "spiky" as the fault-inducing material gets ground against the rolling elements and bearing races.

Other approaches

You can also look at removing what you know about the harmonic content of the signals and then looking at the remainder.  There's some stuff about that here.  However, try something simple like kurtosis first.


QAM constellation coder - transmitting a file The file with size 6 MB is transmitted in 100 seconds using a constellation coder QAM with a constellation of 6 bit/symbol

$ { M = 10^6, B = 8} $

Tasks:

1. What is the signal bandwidth if the band exceedance factor is 0.5?

$W_m = (1+α)R_s$

$W_m = (1+0.5)*80 000 = 120k$

2. What is the signal sampling frequency at the modulator output if the oversampling ratio equals 3.

$ T_S = 6*0.0000020833 = 0.0000125 $

$ n = T_s / T_p $

$ 3 = 0.0000125 / T_p  = 0.00000416666 = 4.17 µs$

$ F_p = 1/T_s => F_p = 239808.15 Hz = 240 kHz  $

3. How numerous is the set of symbols of the constellation used?

$ M = 2^k $

$ k = 6 $

$ M = 2^6 = 64 $

4. What is the modulation performance?

$ R_b / W_m = 480k / 120k = 4 $

5. How long will the signal transmission last if the band exceedance factor is 1.

And this is the only task I have no idea how to do. Any suggestions?

I would also be very grateful I you check the other tasks, if they were done correctly, or if I made any mistakes. Thank you in advance!

Since the bandwidth stays constant at 120 kHz, the symbol rate needs to be reduced: $$R_s = W_m/(1+\alpha) = 120000/2 = 60000 \text{ symbols per second.}$$

Then the bit rate is $R_b = 6R_s=360,000 \text{ bits per second.}$

The transmission duration is then $$T = \frac{48 \times 10^6 \text{ b}}{360 \times 10^3 \text{ b/s}} = 133.3 \text {s}.$$


Constrained LASSO Problem - $ {L}_{1} $ Regularized Least Squares with Linear Equality Constraints I have an optimization question.

I want to solve the following problem:

$$
\arg\min_S\frac{1}{2}\|s-c\|_2^2 +\lambda\|\Phi s\|_1 \mbox{ s.t. } As = 0
$$

in which $\Phi$ is the wavelet transform operator.

My strategy is to find the close form solution of the lasso problem by using the orthogonal invariance property of $\Phi$, and then we project the solution onto the affine space.

I think this method will not converge to the optimal solution, but the numerical results shows good.

Is there any efficient method to solve this kind of problem?

Indeed you can not solve the problem ignoring the equality constraints and then project the solution onto the set of solution for the constraint. It is easy to build real world example which shows that.

Yet, it might be that in most cases it will work reasonably well.

You didn't mention how you solve the LASSO Problem but one of the easiest ways to solve it is using Sub Gradient Descent (Though it is not the fastest, but the idea is simple).

The Sub Gradient is a generalization of the Gradient Descent for cases which are not differentiable like the $ {L}_{1} $ Norm.

The Sub Graident of your function above is given by:

$$ s - c + \lambda {\Phi}^{T} \operatorname{sign} \left( \Phi s \right) $$

The nice thing about the Sub Gradient method that it has the Projected Sub Gradient variant which fits exactly what you need.

What you need to do is each iteration of the Sub Gradient Descent to project the result of the step onto the set of the constraint.

It will be something like:

$ {s}^{ \left( k + \frac{1}{2} \right) } = {s}^{ \left( k \right) } - \alpha \left( {s}^{ \left( k \right) } - c + \lambda {\Phi}^{T} \operatorname{sign} \left( \Phi {s}^{ \left( k \right) } \right) \right) $.
$ {s}^{ \left( k + 1 \right) } = {s}^{ \left( k + \frac{1}{2} \right) } - {A}^{T} {\left( A {A}^{T} \right)}^{-1} A {s}^{ \left( k + \frac{1}{2} \right) } $.

When you need to iterate between those steps until the solution is converged.
It is easy to see that step (2) is projection onto the convex set defined by $ A s = 0 $.

Of course you can replace the Sub Gradient by something more modern and efficient as FISTA and this will still work like a charm (Iterating by Graidnet Step and Projection Step).


Determine minimum sample rate for continuous sinusoid Consider a signal $$ x(t) = \cos(175\pi t) $$ which is sampled to produce discrete time signal $$ x[n] = x(nT_s) $$ The fundamental period of $x[n]$ is $$ N_0 = 7 $$

Given this, what is the smallest possible sample rate $T_s$? (Ans: 1.6327 ms).

I would assume that this is related to finding the Nyquist frequency. I was thinking:

Since, $$N_0 = 7 \implies f_0 = \frac17 \implies f_{\mathrm{Nyquist}} = 2 \frac17\implies T_s = \frac72 $$ However, this is obviously incorrect. I am not even using any information of the original signal. Any suggestions on what I could be doing wrong here?

Remember that the cosine function repeats itself every $2\pi$. The discrete signal can be written as $$x(n)=\cos(175\pi n T_s)$$

So when the argument reaches $2\pi$, we'll have a period. As we know that the fundamental period is $7$, then

$$175\pi \cdot7T_s=2\pi \implies T_s =\frac{2}{7\cdot 175}=1.6327\cdot 10^{-3}$$


Why is $e^{i2\pi Nf_snT}=1$? Why is $e^{i2\pi Nf_snT}=1$ in snippet from book?

Because $Tf{_s}$=1 and $N$ and $n$ are integers.  So the exponent becomes 2$\pi$ times an integer, say $K$. But $e{^{i\theta}} = \cos(\theta)+i\sin(\theta)$. Therefore $e{^{i2K\pi}} = 1$.


How do I perform a time domain phase shift in the frequency domain? I've seen some information on this topic around, but I don't quite understand it.

I have a time domain signal. I understand that if I want to time shift this signal, I can do so by multiplying its Fourier transform by $\exp(-j\omega\delta t)$, where $\delta t$ is the time delay.

However, what if I want to introduce a phase shift in the time domain signal? What would I multiply its Fourier transform by? Or do I have to do this transformation in the time domain?

For context, I am trying to apply both a time delay and phase difference to a time series to correlate with another one, given that the detectors that produced the data are a distance apart and have different orientation towards the source of the signal.

Let’s be clear on what we will refer to as time delay and phase shift. Due to the common association of individual frequencies as sinusoids many confuse delay and phase shift as being equivalent. However a phase shift in time is simply multiplying a time sample by $e^{j\phi}$ while a time delay is displacing the sample in time such as done with $x(t-\tau)$.

The same thing occurs in frequency, where a phase shift of a particular carrier is done by multiplying that carrier by $e^{j\phi}$, while a frequency displacement is done by moving the carrier to a new frequency such as $X(\omega-\omega_c)$.

As you showed, a delay in time is a linear phase versus frequency in the frequency domain.

Similarly due to the bidirectional properties of the Fourier Transform, a displacement in the frequency domain would be a linear phase versus time in the time domain, specifically $e^{j\omega t}$. Multiplying a time domain signal by $e^{j\omega t}$ with $\omega t$ as a linear phase ramp, results in a frequency translation, while multiplying a time domain signal by $e^{j\phi}$ with $\phi$ as a constant, results in a static phase shift (or if it makes it clearer: a phase displacement, or a phase rotation). This phase displacement is not directly a time displacement.

Further, any multiplication in the frequency domain would be equivalent to a convolution in the time domain (and vice versa). If you want to apply a phase change in the time domain; if the intention is equivalent to a multiplication of the time domain samples- this would need to be done by a convolution effectively in the frequency domain.

To provide further details we need to understand more precisely what is meant or intended by a “phase shift in the time domain”, specifically if this means all time samples are rotated by the same phase (which is done by convolving the frequency domain function with a complex impulse at f=0 with magnitude = 1 and associated phase such as $\delta(\omega)e^{j\phi}$ or if the time samples have a linearly increasing or decreasing phase versus time (which is done by convolving the frequency domain function with an impulse at a frequency offset $\delta(\omega-\omega_c$).


Fixed point saturation and round Another question related to fixed point saturation and rounding. Say using notation S<WL, FW> where S means it is signed number, WL is the full bitwidth, and FW is the fractional bitwidth. So S<16, 15> means a 16-bit signed number that has 15-bit fractional bits.
Normally in a N-tap filter, denote input as S<x_WL,x_FW> and coefficients are S<c_WL,c_FW>. Then each multiplication would be S<x_WL+c_WL, x_FW+c_FW> and the N-tap output y would grow S<x_WL+c_WL, x_FW+c_FW> by log2(N) bits to S<x_WL+c_WL+log2(N), x_FW+c_FW>. Finally y's bitwidth S<x_WL+c_WL+log2(N), x_FW+c_FW> would be saturated and round off to S<y_WL, y_FW> where y_FW = x_FW+c_FW-round_off, and the part to the left of decimal point of y would be y_WL-y_FW.
Beside this routine approach, what other ways to limit internal bit growth, for example, in individual multiplication output S<x_WL+c_WL, x_FW+c_FW> and final output? There were some papers on designing the coefficients so that each multiplication can guarantee only grow 1-bit. But what else can be done?
Thanks

I think this overly conservative. If you apply it to an IIR filter, you would conclude infinite bit growth which isn't the case.
For a filter with impulse response, $h[n]$, the "worst case" gain the filter can provide to the input is given as the absolute some the coefficients, i.e.
$$G_{max} = \sum |h[n]| $$
This means that if $x[n] < x_{max}$ then $y[n] < x_{max} \cdot G_{max} $ for all $n$.
You are guaranteed not to overflow as long as you accommodate this maximum gain.
Even that's too conservative for most applications. The input signal that will create this maximum gain is basically $sign(h[-n])$ which is extremely unlikely to occur. So in practice a statistical approach works better: Run a sizable number of representative input signals through the filter and look at the probability distribution function of the output.
At the end of the day this is a trade off between the noise from occasionally clipping an output sample vs the quantization/rounding noise you create on every output sample.


How does noise reduction for speech recognition differ from noise reduction that is supposed to make speech more "intelligible" for humans? this is a question that has interested me for some time now, mainly because I'm working on noise reduction for an existing speech recognition system myself.

Most papers on noise reduction techniques seem to focus on how to make speech more intelligible for humans, or how to improve vague terms like "speech quality".

I'm sure that, using criteria like these, you can identify filters that make noisy speech signals easier to listen to for humans.
However, I am not sure that these criteria can simply be adapted when trying to evaluate speech signals that have been denoised to improve the accuracy of speech recognition system.

I don't really find papers that discuss this difference. Do speech intelligibility and speech quality correlate with the accuracy of speech recognition systems? Are there objective measures that can evaluate how "good" a denoised speech signal will be for a speech recognition system, for example if also given the original clean speech? Or is the only way to find out how good your noise reduction technique is, to train the speech recognition system on the denoised data and look at the accuracy?

I'd be happy if someone could point me into the right direction, or maybe give some papers that discuss this. Thanks in advance!

I don't really find papers that discuss this difference.

There are whole books on the subject:

Robust Automatic Speech Recognition 1st Edition

  Do speech intelligibility and speech quality correlate with the accuracy of speech recognition systems?

Usually no, usually noise reduction corrupts features in unpredictable way and reduces speech recognition accuracy.

  Are there objective measures that can evaluate how "good" a denoised speech signal will be for a speech recognition system, for example if also given the original clean speech? Or is the only way to find out how good your noise reduction technique is, to train the speech recognition system on the denoised data and look at the accuracy?

Second. Moreover feature-based noise reduction actually removes important information from the spectrum altogether so you can not repair an accuracy of the clean system. For that reason modern approach is to perform multi-style training on noisy data instead of using noise reduction algorithm beforehand. It ends in more accurate recognition.