Using MFCCs for acoustic machine failure prediction MFCCs are ubiquitously extracted for speech processing tasks, but I would like to know how suitable they are for non-speech processing tasks. Intuitively, it is my understanding that MFCCs are powerful features because they are more representative of how the human ear "hears" things. However, for tasks that are seemingly dissimilar from those that would have warranted such a biological adaptation, should MFCCs still be used? For example, in acoustic machine failure prediction, what is the logic behind using MFCCs, if any? Such a domain appears to be one that the human ear shouldn't have any significant advantages in.

Perhaps the human ear has evolved to better identify subtle changes in sound for hunting, or detecting predators, making MFCCs relevant to anomaly detection tasks.

Any thoughts? Should MFCCs be used universally?

The biological foundation for MFCCs are pretty weak. The mel-scale, which the mel-spectrogram and thus MFCC uses, is a rough approximation of how frequency differences are perceived along the frequency axis. The log transform of amplitudes is better than linear, but not really a good model of the ears loudness perception - which is non-uniform both with amplitude and frequency, as well as state-dependent (temporal masking).
The choice of hop-lengh is also maybe roughly close to what the human ear has, but vastly simplified - the ear temporal resolution is not uniform either. And no binaural representation either, a critical part of the human auditory system.
The primary reason for using MFCCs is its numerical convenience. The log mel-spectrogram (when computed with appropriate parameters) preserves a large amount of the information in an audio signal - but reduces the feature space considerably and introduces useful structure (time-frequency). Using a 32 ms hop, 512 samples at 16kHz = 32 frames per second, and 64 mel bands one can reduce one second of audio from 16000 dimensions to 64x32=2048. Then when computing the DCT-II to get the MFCC, most of the energy of this ends up in the lower coefficients. This lets us drop the higher coefficients without loosing much information. Typical would be 20 (or 13) coefficients, leading to 20x32=640 dimensions. The information in the different coefficients also tends to be much less correlated than the input spectrogram, an advantage for some models - such as the commonly used GMM-HMM (Gaussian Mixture Model, Hidden Markov Model).
Bottom line is that the feature space with 640 MFCCs per second, instead of 16000 samples per second, is much easier to learn. This dimensional reduction and feature transformations of MFCC is very useful, even if the input is not speech. So one can make OK performing systems for many non-speech tasks using MFCC.
That said with the advent of deep learning and larger amounts of data and compute - MFCC based models are being outperformed by log mel-spectrograms processed by 2d Convolutional Neural Networks, as well as raw waveforms processed by 1d CNNs.
Note also that in the majority of cases, we are still struggling to achieve human level performance. And with supervised learning, our labels are typically provided by humans - a performance roof which the model cannot easily surpass. In this scenario it is not that problematic that our audio reception is limited by loose inspirations from the human auditory system.
There are notable exceptions though, such as ultrasound (for bat detection etc) and infrasound (for condition monitoring) - there it is critical to go beyond the frequency range that human hearing has.


sampling a sinusoidal signal at the Nyquist rate To sample a sinusoidal signal at the Nyquist rate, we have to take two samples per cycle. Is there any condition on the position of the two samples? is it necessary to take them at the peaks?

Noha, when the classical sampling theorem (Nyquist-Shannon) for baseband signals
is applied to an ideal sinusoidal of infinite extend of frequency $f_m$ Hertz
$$x(t) = \cos(2 \pi f_m t)$$ it's found that the signal must be sampled at a sampling frequency greater than $2 f_m$; i.e.
$$f_s > 2 f_m$$

Note that even though the selection of $f_s = 2 f_m$ would sufficiently work for most of the typical engineering signals, it will not work for an ideal sinusoidal signal. The reason is that the ideal sinusoidal has an exact impulse at the $f_s / 2$ frequency whereas most typical signals will not.

The practical solution is to use a sampling frequency slightly greater than the theoretical value given by the sampling theorem unless you are sure that there wont be any components at $f_s/2$.


Random Process at a particular time instance I was studying Random Process and I thought I understood what it was all about until I came across this example. Consider a random experiment of tossing a coin with sample space S = {H, T} The sample functions are given by:

$X(t, H) = x_1(t) = \sin(\omega_1t) \tag{1}$

$X(t, T) = x_2(t) = \sin(\omega_2t) \tag{2}$

Now, I understand that random process at any particular time instance say $t_k$ gives me a random variable. Also, the sum of probabilities of a random variable should be 1 i.e.
$$X(t_k,H) + X(t_k,T) = 1 \tag{3}$$
subsequently in general no matter where I sample the random process I should have,
$$\sin(\omega_1t) + \sin(\omega_2t) = 1 \tag{4}$$

But the Equation 4 won't be true for all values of t.

Question 1. Where am I going wrong above?

Question 2. If I take the value of $X(t,\lambda)$ at a particular time $t_k$and at a particular $\lambda_i$ the value $X(t_k,\lambda_i)$ will be a number. This number is equal to the probability of $X(t,\lambda)$ taking value $\lambda_i$ at time instance $t_k$. Am I correct?

Your equation (3) and (4) are wrong.

When you say

  Also, the sum of probabilities of a random variable should be 1

you actually mean the total sum of the values of probability mass function $P_{X_k}[x]$ equals one; i,e.

$$ \sum_{k} P_X[x_k] = 1 $$

(or the integral of the probability density funtion $f_{X_k}(x)$ equals one as
$ \int_{-\infty}^{\infty} f_X(x) dx = 1 $ )

But in your equation (3) and (4), you are not summing the values of probability mass function but the values of the sample function itself, which leads to te wrong conclusion and confussion.

To reach the correct conclusion, in your example RP, the pdf associated with each RV $X_{t_k}$ is
$$
f_{X_k}(x) = P(H)  \delta( x- \sin( \omega_1 t_k) ) + P(T) \delta(x- \sin( \omega_2 t_k))
$$

which integrastes (sums) to one as

$$
\begin{align}
\int f_X(x) dx &= \int \left( P(H)  \delta(x- \sin( \omega_1 t_k)) + P(T) \delta(z- \sin( \omega_2 t_k)) \right) dx \\
&= P(H) \int \delta(x- \sin( \omega_1 t_k)) dx + P(T) \int  \delta(x- \sin( \omega_2 t_k)) dx \\
&= P(H)+P(T) = 1
\end{align}
$$


Calculation of inverse impulse response in the frequency domain I want to calculate the inverse impulse response of a LTI system in the frequency domain. I generate a simple impulse response $g$. For this I generate a vector of 100 zeros. I set the value of the 20th sample to 0.5.

I calculate the frequency response function by calculating the Fourier-Transformation $h=\mathcal F(g)$.
I invert every single frequency bin in $h$ where $h_{inv}[k]=\frac{h[k]^*}{|h[k]|^2}$ (This is the formula for the inverse of a complex number)
I do a inverse Fourier-Transformation $g_{inv}=\mathcal F^{-1}(h_{inv })$
I plot $g$ (blue) and the calulated inverse $g_{inv}$ (red)

Question: Why is the impulse of $g_{inv}$ at 82 samples? Is $g_{inv}$ my inverse impulse response?
Matlab code:
g=[zeros(1,19) 0.5 zeros(1,80)];
h=fft(g);

h_inv=conj(h)./(abs(h).^2);
g_inv=ifft(h_inv);

figure
plot(g)
hold on
plot(g_inv)

Your result is correct, even if it seems a bit counterintuitive. Note that you don't compute the Fourier transform but the discrete Fourier transform (DFT). The system you want to invert has the impulse response
$$g[n]=\frac12\delta[n-19]\tag{1}$$
Obviously, the inverse system must advance the incoming signal by $19$ samples:
$$g_{inv}[n]=2\delta[n+19]\tag{2}$$
Clearly, that system is non-causal. Taking into account the implicit periodicity of the DFT (and IDFT) with length $N=100$, the signal $g_{inv}[n]$ cannot be distinguished from a signal $g_{k,inv}[n]=g_{inv}[n+kN]$, $k\in\mathbb{Z}$. Your result is just $g_{-1,inv}[n]$:
$$g_{-1,inv}=2\delta[n+19-N]=2\delta[n-81]\tag{3}$$


Region of root locus not inside the feasible region for a given performance specification I wanted to plot the root locus for a system, whose plant transfer function is
$$
P(s)= \frac{k}{Ts+1}e^{-ts}.
$$
A PI controller was suggested for this system, and the task is to find $K_p$ and  $K_i$ values, such that it satisfies the following performance specs:

Maximum peak overshoot < 5%
Settling time (2% band) < 0.8 seconds

This is how I modified the open loop transfer function to suit the format of root locus finding.
$$
G(s)H(s) = C(s)P(s)= K_p\frac{k(s+A)(2-ts)}{s(Ts+1)(2+ts)},
$$
where $A = K_i/K_p$ and Pade approximation was used. For starters, $A$ was assumed to be 1.
But, when I plot the root locus of the above, this is what I get.

Here, the white region is the feasible region which satisfies the performance specs.
As we can see, the root locus of the pole introduced by the integrator (at Im axis) ends up being completely inside the unfeasible region. This is the case for any $A$ value I take.
What should I do? Is my approach correct? And more importantly, how will this pole affect my performance?
I'm new to this field, and I'm a bit confused by this. Please help me out here. Thanks in advance.

The root locus will start at the open loop poles and move toward the open loop zeros as the gain is varied. Your pole cannot get past the zero that is located in the yellow region. This pole that is close to the imaginary axis will be dominant, meaning it will dominate your response times, slowing down the entire system. One approach is to cancel that zero with an additional open loop pole at it's location, unless that zero was introduced by your controller in which you have the ability to move it to be out in the white area, or eliminate it entirely.
It appears that the zero that terminates the closed loop pole as shown is due to $(s+A)$ in the numerator, in which case choosing a value for A that is within the white region should allow the pole to extend further to that point as the gain is increased. Also, the pole and zero due to the Pade approximation does not appear to be visible in the root locus.
To see how poles close to the imaginary axis are dominant, consider the simple case of $x(t) = e^{-at}$ which has a Laplace Transform $X(s)=\frac{1}{s+a}$  (which is a pole at $s= -a$. The time domain equation is a decaying exponential with time constant $a$, if a is very small this will take a very long time to decay. If $a$ is very large, the decay rate will be very fast.


Inverse Fourier transform - porting from MATLAB to Python I am porting a script from MATLAB to Python, but I am failing when it comes to the inverse Fourier transform.

The function in MATLAB (ifft) includes a 'symflag', which treats the data as conjugate symmetric and ensures that the output is real. However, the equivalent function in numpy does not have a similar option, and therefore my outputs are different.

The 'symflag' clearly covers up for some of my mathematical inadequacies, and I would be grateful for assistance in resolving this issue. My aim would be to have a solution in Python that closely matches the solution in MATLAB.

In Matlab

Y = [complex(1), complex(2), complex(3), complex(3), complex(2), complex(1)]

y = ifft(Y, 'symmetric')

y = 2.33    -0.50   -0.16   0   -0.16   -0.50

In Python

Y = np.array([1, 2, 3, 3, 2, 1], dtype='complex')

y = ifft(Y)

y = ?

Many thanks in advance.

For a vector Y of length $N$, where $N$ is even, conjugate symmetry requires the following property to hold:
$$
Y[k] = Y[N-k]^*
$$
for $k=1,...,N-1$ and $Y_0$ is real. (See this.) Here $^*$ denotes the complex conjugate.

The array [1,2,3,3,2,1] this conjugate symmetry property doesn't hold. What did Matlab do when you used 'symmetric' with an input vector that wasn't conjugate symmetric? It pretends [1,2,3,3,2,1] is nearly conjugate symmetric$^\dagger$ and just replaces it with the closest conjugate symmetric vector, which is, [1,2,3,3,3,2].

If you try this in Python:

Y = np.array([1,2,3,3,3,2])
y = np.fft.ifft(Y)

your answer will match Matlab's answer with the 'symmetric' option for [1,2,3,3,2,1]. Note there there will be some rounding issues and you can verify that the imaginary part of y in Python is practically zero. Here's my output:

In [8]: np.fft.ifft(np.array([1,2,3,3,3,2]))
Out[8]:
array([  2.33333333e+00 +0.00000000e+00j,
        -5.00000000e-01 -1.48029737e-16j,
        -1.66666667e-01 +0.00000000e+00j,
         1.11022302e-16 +2.90348171e-16j,
        -1.66666667e-01 +0.00000000e+00j,  -5.00000000e-01 -1.42318434e-16j])

In [9]: np.real(np.fft.ifft(np.array([1,2,3,3,3,2])))
Out[9]:
array([  2.33333333e+00,  -5.00000000e-01,  -1.66666667e-01,
         1.11022302e-16,  -1.66666667e-01,  -5.00000000e-01])

In [10]: np.imag(np.fft.ifft(np.array([1,2,3,3,3,2])))
Out[10]:
array([  0.00000000e+00,  -1.48029737e-16,   0.00000000e+00,
         2.90348171e-16,   0.00000000e+00,  -1.42318434e-16])

And now to your main question: How do you port this to Python?

Step 1: Take the input vector Y=[1,2,3,3,2,1]

Step 2: Force it to be conjugate symmetric by looping through and changing the right half of the vector to be equal to the complex conjugate of the left half. Call this modified vector Y_new = [1,2,3,3,3,2].

Step 3: Compute np.fft.ifft(Y_new) and throw away its imaginary part. (You can do an assert and ensure the imaginary part has a very tiny norm.)

Alternatively, when $N$ is even, Steps 2 and 3 can be replaced by a quicker np.fft.irfft(Y[0:N/2+1]) as endolith suggests in the comment below.

See also: https://blogs.mathworks.com/steve/2010/07/16/complex-surprises-from-fft/ for an example when $N$ is odd.

$^\dagger$ AFAIK, mathematically, there is no such thing is "nearly conjugate symmetric." A vector is either conjugate symmetric or not! I wish Matlab's documentation was more explicit about what it is doing under to hood.


Is there an equivalent of Parseval's theorem for wavelets? Parseval's theorem can be interpreted as:

  ... the total energy of a signal can be calculated by summing power-per-sample across time or spectral power across frequency.

For the case of a signal $x(t)$ and its Fourier transform $X(\omega)$, the theorem says:

$$
\int{|x(t)|^2 \; dt} = \int{|X(\omega)|^2 \; d\omega}
$$

For the case of discrete wavelet transform (DWT), or wavelet packet decomposition (WPD), we get a 2D array of coefficients along the time and frequency (or scale) axis:

        |
        | c{1,f}
        | ...
freq    | c{1,2}
        | c{1,1} c{2, 1} ... c{t, 1}
        |______________________________
                 time

Can a sum of this series somehow be understood as a signal's energy? Is there an equivalent rule to Parseval's theorem?

Yes indeed! In theory as long as the wavelet is orthogonal, the sum of the squares of all the coefficients should be equal to the energy of the signal. In practice, one should be careful that:

the decomposition is not "expansive", i.e. the number of samples and of coefficients is the same.
wavelet filter coefficients are not re-scaled, as happens in some applications (like lifting wavelets, to keep integer computations).
wavelets are orthogonal (this is not the case in JPEG2000 compression).

You can verify this indirectly, looking at approximation coefficients. At each level, their number of samples is halved, and their amplitudes have around a $1.4$ scale factor, which is just $\sqrt{2}$. This feature is used for instance to estimate the Gaussian noise power from wavelet coefficients:

$$ \hat{\sigma} = \textrm{median} (w_i)/0.6745$$

A little further, there is a notion that generalizes (orthonormal) bases: frames. A set of functions $(\phi_i)_{i\in \mathcal{I}}$ ($\mathcal{I}$ is a finite or infinite index set) is a frame if for  all vectors $x$:

$$ C_\flat\|x\|^2 \le \sum_{i\in \mathcal{I}} |<x,\phi_i>|^2\le  C_\sharp\|x\|^2$$

with $0<C_\flat,C_\sharp < \infty$. This is a more general Parseval-Plancherel-like result used for general wavelets.

In other words, it "approximately preserves energy" by projection (inner product). If the constants $ C_\flat$ and  $C_\sharp $ are equal, the frame is said to be tight. Orthonormal bases are non-redundant sets of vectors with $ C_\flat=C_\sharp = 1 $.

For those using Matlab, you should care about the native border extension, which is obtained by dwtmode('status'). Some add tails to the data to help inversion with little border artifacts. With a periodic mode dwtmode('per') and a number of samples that can be divided by $2^L$ where $L$ is the wavelet level, you can get a good match in energy, with tiny differences:


What is the difference between the V-OFDM and OSDM I am trying to understand the difference(s) between V-OFDM (Vector Orthogonal Frequency Division Multiplexing) and OSDM (Orthogonal Signal Division Multiplexing) which are used as a multicarrier modulation techniques in wireless communication. Based on what I understand both are dividing the sub-carriers into  groups of sub-carriers vertically, and then process the $DFT$ horizontally. But when checking online and read different papers and books, some are talking about V-OFDM and others OSDM !!

Could anyone please explain further about differences between them ? design and differences  of OSDM from V-OFDM? Or they are similar, only two different names with same technology !!!

Thanks in Advance!

They are little bit different in performing the $DFT$ and the adding of $CP$. her you can read more details about those differences : https://sci-hub.tw/10.1109/ICC.2009.5198846

Good luck


Characterising a system from its input and output I have 2 data sets sampled at regular intervals which represent the input and the output of a system. A quick analysis of the plots shows the following:

There is a proportional factor between the input and the output.
There is a delay before a change in the input is observable in the output. This delay depends on the scale of the change: larger changes in the input are observable in the output after a longer delay than smaller changes in the input.

I would like to approximate the system with a function, but I am novice in signal processing and I have no idea how to get started. I don't even know if this is a solvable problem. Any idea? I am comfortable using Matlab/Octave, and computer programming in general.

I can get more data from the system if required, but I cannot control the input. This is because the system is only observable when a PID controller keeps it stable.

Use the matlab cpsd (cross power spectral density) function on your two signals.
That will give you the frequency response.
Do an xcorr (cross correlation) to find the time delay between the two signals.
Do this for several pairs of input and output files and see if the results match between different pairs.
If they are (reasonably) similar then you could create an FIR filter that reproduces the frequency response and delay.

You could then apply the filter to a set of input data to approximate the expected output.

If there is a PID controller in the middle, however,  I don't know if this procedure will be enough.


Fourier Transform of $u(t)$ I am just unable to find the correct Fourier transform of these signals (unit step, sine and cosine functions) which are containing delta functions in their Fourier transform.
For unit step function,
I proceeded as -->
$$ U(ω)=\int_{-\infty}^{\infty}u(t)e^{-jωt}dt $$
$$=\int_{0}^{\infty}e^{-jωt}dt$$
$$=\frac {e^{-jωt}}{-jω}\Big|_0^{\infty} $$
$$=\frac {1} {jω}$$
But as I checked this is just a part of the answer as the Fourier transform also includes $\pi δ(ω) $.Now how do I get this? Same is the case for sine and cosine signals. Can anybody derive the answer in simple terms?

Your derivation doesn't give the correct result because the limit $\lim_{t\to\infty}e^{-j\omega t}$ doesn't exist. You seem to have confused that limit with $\lim_{t\to\infty}e^{-\omega t}$, which does exist and which equals zero for $\omega>0$.
There are many ways to derive the Fourier transform of the unit step function $u(t)$ (see e.g., this answer), but I'll show you one which doesn't involve more than just engineering math. At least if you accept some intuitively pleasing facts such as that the derivative of the unit step equals the Dirac delta impulse:
$$u'(t)=\delta(t)\tag{1}$$
It's important to note that $(1)$ holds for any function with a discontinuity at $t=0$ that is otherwise constant, satisfying $u(0^+)-u(0^-)=1$. Hence, $(1)$ only defines $u(t)$ up to a constant. Taking the Fourier transform of $(1)$ gives
$$j\omega U(\omega)=1\tag{2}$$
It would be tempting to conclude from $(2)$ that $U(\omega)=1/j\omega$ holds, however, this is not correct. The reason is the additive constant that got lost due the derivative in $(1)$.
Hence, the correct conclusion from $(1)$ is
$$U(\omega)=c\delta(\omega)+\frac{1}{j\omega}\tag{3}$$
where $c$ is the constant that can be added to any function satisfying $(1)$ without violating Eq. $(1)$. We just have to find the specific value of $c$ for which $(3)$ holds, with $U(\omega)$ being the Fourier transform of the unit step.
It's easy to see for which function $c=0$ holds. Since $1/j\omega$ is purely imaginary, we know that the corresponding time domain function is odd, it has a jump of height $1$ at $t=0$, and it is otherwise constant. The only function satisfying all of these requirements is a scaled signum function, i.e.,
$$\frac12\textrm{sgn}(t)\Longleftrightarrow \frac{1}{j\omega}\tag{4}$$
Since the unit step function satisfies
$$u(t)=\frac12+\frac12\textrm{sgn}(t)\tag{5}$$
we finally arrive at the Fourier transform of $u(t)$ by combining $(4)$ and $(5)$:
$$\bbox[#f8f1ea, 0.6em, border: 0.15em solid #fd8105]{
U(\omega)=\pi\delta(\omega)+\frac{1}{j\omega}}\tag{6}$$


Invariance of $T[x[n]]=\sum_{k=n-1}^{n+2}x[k]$ I have to test whether the following system is invariant or not: $T[x[n]]=\sum_{k=n-1}^{n+2}x[k]$, so I want to verify that, if $y[n]=T[x[n]]$, then $y[n-N]=T[x[n-N]]$.
$$T[x[n-N]]=\sum_{k=n-1}^{n+2}x[k-N]$$
$$y[n-N]=\sum_{k=n-N-1}^{n-N+2}x[k]$$
In the first equation, if $k'=k-N$, we have that $T[x[n-N]]=\sum_{k'=n-N-1}^{n-N+2}x[k']$, so both sums are equal and thus the system is invariant. Is this correct or am I making wrong assumptions? Thanks in advance!

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

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

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

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


Created Matlab echo signal: Is it an actual "echo?" Or am I just manipulating the signal to make it seem like an echo? This is just an outline of what I did (assume that all my dimensions fit):

    (code that imports audio to read audio as function x)
    X=fft(x);
    (code to define time t and frequency f)

    H=L+M*exp(-1i*2*pi*f*N)+O*exp(-1i*2*pi*f*P);
    % This is where it gets tricky
    % In the frequency domain I have the function above where f is the frequency
    % variable and L M N O P are arbitrary values

    Y=X.*H;
    % again assume dimensions fit

    y=real(ifftshift(ifft(Y)*length(t)));
    %Taking the inverse fft and fftshift was the only way I could get back a signal
    % Once again assume the dimensions fit

Here L M N O P = [.9 .75 1.0 .5 .5]

Is this legit? or am I just pulling my own leg?
Is it an actual "echo?" Or am I just manipulating the signal to make it seem like an echo?

To generate a repetitive echo all you need to do is pass your signal through an appropriate AR filter with a delay corresponding to the echo period. For instance the following code generates (for a damped 1 Hz sinusoid sampled at 100 Hz) an echo which repeats every 1000 samples. The strength of the echo can be adjusted by the AR filter coefficient (0.5 here)

Example Matlab code

% generate 10000 samples of a damped sinusoid at 100 Hz
sig = sin(2*pi*(0:9999)/100).*exp(-(0:9999)/200);

% AR filter to generate a repetitive echo, period 1000 samples (10 seconds)
b = 1;
a= [1 zeros(1,999), 0.5];

echo = filter(b,a,sig);

You may want to experiment with your AR filter to match your echo decay. Also look into 'gapped deconvolution'


How best to digitize a noisy phone recording for filtering? How do I digitize an analog telephone recording and filter it with a Chebyshev filter?

Unless you have tight space constraints, I would stick to lossless formats like good old WAV. Phones do not pass much high-frequency content so you should be fine with a 8KHz Mono WAVE file. You will however want to keep the bit depth high, say to 24 bits, since you want to process the audio. Once you are finished with your processing you can compress with a lossy codec, such as OGG.

All codecs (coders/decoders) may be classified as lossy or lossless. Lossless codecs, as their name implies, allow for perfect reconstruction of the original signal. Lossy codecs, by contrast, merely approximate the source. By forgoing perfect fidelity, lossy codecs enable higher compression ratios, commensurate with the amount of degradation one is prepared to accept.

Perhaps related

A non-lossy aka lossless codec without patent-restrictions for audio?


Radix-4 FFT versus Radix-2 Is a radix-4 implementation faster than a equivalently well coded radix-2 FFT?  And if so, why would it be faster?

It depends. Theoretically you can save a few multiplies with a radix-4 as radix-4 has a 1/4th the number of butterflies and 3 mpy + 8 adds per butterfly (if properly structured) and the radix 2 has 1 mpy + 2 adds per butterfly.

So in terms of multiplies it's a bit better, however there is higher complexity in terms of code structure, exception handling, coefficient management, register management, digit-reverse addressing, etc.

So it's only an advantage if the number of mpy is the limiting factor which for most hardware these days is not the case.


Is Parks-Mcclellan algorithm guaranteed to stop? Is there any guarantee that Parks-Mcclellan eventually stop? If so, is there any formula that shows number of maximum interations that ensure that the algorithm stops?

You have to distinguish between the algorithm and a specific implementation. The Parks McClellan algorithm is based on the Remez exchange algorithm, which is guaranteed to converge. So in principle any correct implementation of the algorithm should converge. However, convergence of the algorithm is only guaranteed if all intermediate results are computed with infinite precision, which is not the case in a practical implementation. Some design specifications can result in numerically ill-conditioned problems that will prevent correct convergence of the algorithm.

In short, yes, the algorithm is guaranteed to converge. However, in some cases, numerical problems can prevent convergence to the optimal solution. But even in the usual case  of convergence, the number of necessary iterations cannot be determined beforehand.


Why should the cut off frequency and the transition bandwidth (BW) be expressed as a fraction of sampling rate? I came across this confusion in Chapter 16,"Book - The Scientist and Engineers guide to DSP" and Topic is "Windowed Sinc-Filters". I am doing some calculations for designing the sinc filter and came across these doubts:

                               M = 4/BW

where "M" is the filter kernel length and "BW" is the transition bandwidth. The author says that "The BW must be expressed as a fraction of sampling rate with a value between 0 to 0.5 and same goes with the cut off frequency". For the values of M = 20, 40, 200, the BW is 0.2, 0.1, 0.02 and like we can see, it satisfies the condition of 0 to 0.5.

After I searched online, I understood that the value of cut off frequency and BW to exist between 0 to 0.5 has to do something with the power (dB) and presumed that 0 is 0db and 0.5 is -3 dB, since 10log(0.5) = -3dB.

My question is, what is the relationship between the cutoff of frequency and BW being a fraction of the sampling rate and between the values 0 to 0.5? What does it mean?. Is it something to do with the frequency components present beyond the cut off freq or BW being attenuated?

For example: If my sampling rate is 100sps, so the cut off frequency and BW must be a fraction of 100 sps with a value between 0 to 0.5?

It is just a normalization so that the apparent sampling frequency is 1 Hz. Say for example you have a 101 tap lowpass filter with with a cutoff of 10 Hz at a sampling frequency of 50 Hz.

If you don't change the filter coefficients but apply the exact same filter at a sampling rate of 500 Hz then the cutoff will be 100 Hz (not 10 Hz) because the frequencies scale according to the sampling rate.

An important point to note - consider the original filter has a transition width of 10 Hz i.e. from 10 Hz to 20H at the 50 Hz sampling rate. If the filter coefficients remain unchanged when applied in a system with a sampling rate is 500 Hz, the resulting transition width is now from 100Hz - 200 Hz, so the transition width has also increased. This is the reason that if you are trying to design a filter with a narrow transition width at high sampling rates it  requires more filter taps i.e. the impulse response of the filter gets longer in terms of the number of samples.

In most of the filter design software they accept normalized frequency values for the specification - but be careful some use slightly different normalizations e.g. I believe Matlab expects the frequencies to be between 0 and 1 rather than 0 and 0.5, their filter design interface has changed over the years.

Also, if you look at the formulas in the literature for estimation of the FIR filter length for routines like the Parks- McClellan algorithm, you'll see that the frequencies are normalized with respects to the sampling frequency.


How do I construct input to neural network from audio signals? Input: Microphone recordings of digits from 0 to 9 from different speakers.

Output: The digit from 0 to 9.

I am doing this for fun. So first I will train my neural network using some samples and then use it to classify digits.

Problem is every person takes different amount of time to say the digit also every person takes different amount of time to say different digits. This gives me sound signal of different time duration. But input to the neural network  must be of some fixed size. How can I make all my input signals of same length so that I can give input it to neural network? Is there any standard technique? Fourier transform?

So the problem is that you want to classify some audio samples but each sample has different sizes. AFAIK there is no special way for doing this with (simple) neural networks.

For general classifiers I am sure there are lots of ways but these are some very simple ones so you can begin playing with the problem:

Divide your signal by frames/segments of equal size and use each frame as if it were a training example. At the end you will use a voting strategy. Since you are using a neural network, you can use the probabilistic outputs of the last layer instead of the hard classes to weight this voting.

Using the average across all frames/segments in your signal. This is a very light-weight approach but can have good results sometimes.

Either way in the "feature extraction" part you will want to use some spectral-based transformation. I recommend either the logarithm of the spectrum or computing MFCC vectors. You can also train your network from raw audio but this is a bit unusual.

You should note that these two methods completely spoil the dynamics of your system because of the granularity. Try playing with the frame sizes or aggregating neighbor frames to increase or decrease the granularity.


Noise modeling and reduction in Video when we have just noisy video We would like to send a video (sequence of images), we have limit bandwidth. If there is a lot of noise, the required bandwidth (bitrate) is needed to send the video+noise becomes high that is not suitable. We have just a noisy file (i.e. we do not have noise separately). The solution is to model the noise. we do not have any other information rather than a noisy video. How we can reduce the noise? How we can model the noise?

Reduction of Noise in a video depends on what type of noise is present in the video.
The noise can be reduced by filters which are all present [There are software for this] and many type of filters such as median ,mean,wiener filter.
And The noise model can be build to add noise to the video there are many matlab code as well as software also there.
Gaussian ,salt and pepper,poisson and speckle noise model can be build and by adding the image sequence with the video you will get a noisy video.

imnoise(video,noisename);

this will give noisy video as output.

If you have any doubt ask me.


How to decipher this Microphone sensitivity Chart for Choosing a SPL meter I have a mic having sensitivity chart as below which has relative dB on y-axis showing

  1V/0.1Pa

what does this convey how to read my mic's sensitivity

a document on microphone sensitivity from analog says the microphone sensitivity is measure with 1Pa pressure as reference

confused in half way after seeing my datasheet of mic

usually i was of the opinion that the microphone sensitivity is measure with reference to 20micropascals,but i read with this SPL measurements
so i am wrong

so if i have to buy an SPL i have to correlate my SPL with the microphone i am using

because if my microphone i am using in my application picks up some sound according to its sensitivity and the spl meter i am purchasing may not give me the required SPL level that i have detected

when i have seen the frequency Weighting curves of the SPL meter they are shown as below

what do these convey ? are they trying to attenuate some frequencies or they are trying to follow human auditory model

so if i have to buy an SPL meter having correlative sensitivity as of my microphone would be better,what is your opinion am i going the right way in selection of an SPL

Microphones are typically specified in dB re 1 V/Pa. So am microphone with a sensitivity -44 dB (like this one http://www.panasonic.com/industrial/electronic-components/parametric-search.aspx?src=/www-cgi/jvcr13pz.cgi?E%2BSS%2B2%2BABA5010%2B0%2B4%2BAM) would produce 6.3mV RMS at 1 Pa of sound pressure. 1 Pa also corresponds to 94 dBSPL. So you can easily convert dB SPL in dB Pa by subtracting 94.

The datasheet for our microphone looks strange. I have never seen unit of 1.0V/0.1Pa before. This is potentially a print error.

The weighting curves do indeed mimic the frequency response of the human auditory system at different overall levels. However they are all 0 dB at 1 kHz, so if you calibrate at 1 kHz the weighting doesn't matter.


Fourier Transform negative amplitude meaning I am reading this example http://www.thefouriertransform.com/pairs/truncatedCosine.php

What does it mean to have some of the frequency components be negative in its amplitude ? I am not talking about the negative frequencies.

The fft returns complex values, to get the amplitude you need to take the abs( ). The real and imaginary portion tell you about the signals phase. Remember the fft is changing the basis by projecting your signal onto a complex sinusoid:
$$e^{i \omega t} = \cos(\omega t) + i \sin(\omega t)$$

and thus your signal is now a set of complex sinusoids which have some phase and amplitude. Think about the phase of a vector $v = [a, \, i\cdot b]$ on the complex plane and what this would mean.


Intuition behind impulse response terms in convolution We're learning about convolution in my signals and systems class right now. I have been able to do all of the problems by simply working out the respective sum/integral, but I'm still having trouble gaining the intuition behind it.

Consider the following example. Let $x[n]$ be a discrete-time signal and input it into some LTI system with impulse response $h[n]$. Then,

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

Let's plug in some values to make this more concrete. Suppose that we want to compute $y[6]$. Well,

$$
y[6] = \sum_{k=-\infty}^{\infty} x[k] h[6-k] = \cdots + x[4]h[2] + x[5]h[1] + x[6]h[0] + x[7]h[-1] + \cdots.
$$

I understand the shifting, but I feel as if the multiplications should be in a different order. Namely, why are we multiplying $x[7]$ by $h[-1]$. I feel as if we should be multiplying it by $h[1]$, since we've essentially shifted everything to the right by $6$ units, to $6$ is the new $0$, which would mean that $7$ is the new $1$ (under the shifting). I have it backwards, and understand why when I work out the math, but why?

More concretely, my question is the following:

With respect to the above example, what exactly is the meaning of $h[2]$? $h[-1]$? $h[k]$ in general?

The coefficient $h[n]$ is the value of the system's response at time $n$ when the input signal was an impulse at $n=0$. Obviously, that's why we call $h[n]$ the impulse response. From this you can see that for a system to be causal, $h[n]$ must be zero for $n<0$, otherwise the system would "know" in advance that an impulse will come at $n=0$.

Note that any discrete-time signal $x[n]$ can be written as a sum of unit impulses:

$$x[n]=\sum_{k=-\infty}^{\infty}x[k]\delta[n-k]\tag{1}$$

Since the response to a shifted impulse $\delta[n-k]$ is $h[n-k]$, and since the system is linear and time-invariant, from $(1)$ the output signal must be given by

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

which is of course the discrete-time convolution. Consequently, positive indices of $h[n]$ correspond to the memory of the system, that's why $h[1]$ is multiplied with the past input sample $x[n-1]$, etc.


What is an "annason" (VLF whistler)? I am watching the following video on software defined radio and at the beginning the narrator points out some "whistlers" on the spectrum display and says maybe they are caused by an "annason". He goes on to say that he knows "commercial annasons scan at a slower rate than the whistlers". What is an "annason". I don't know the spelling because it is verbal.

http://www.youtube.com/watch?v=ZuRcaxpbYCw

I figured out the answer. It is an "ionosondes".


mne raw plot looks like vertical lines instead of the usual signal I have a data from a single eeg channel. When I plot it with plt, it looks like a regular signal.

However, when I'm trying to plot with mne raw object's method plot, it looks very strange, like vertical lines:

Does anyone know why it looks like that?
Thanks!

You need change the scalings parameter. In this case, you need increase it.


Simulating the impulse response of a multirate digital filter I have designed a CIC filter in Verilog and I would like to test its impulse response, but I am not sure where to start from.

My input sampling frequency is 51200Hz, the filter is a third order decimating CIC with downsampling ratio = 256, so the output sampling frequency is 200Hz. The differential delay was set to one.

The input word is one bit wide, while the output will be 16 bits. I have already calculated all the internal registers dimensions and I am satisfied with that, my current Verilog design does not use true registers but integer numbers so that is not a problem anyway.

How am I supposed to compute its impulse response?

First of all there's the problem that the sequence '1000...' is not an impulse, because 1 corresponds to +1 while 0 corresponds to -1. I have solved that doing two simulations, one with '10101010...' and the other with '11010101010...' as input sequences. The impulse response is then half of the sum of the two responses. So far so good.

Since I want the frequency magnitude response referred with the output sampling frequency I figured that my input impulse should be 256 high frequency samples wide, is that correct? I have done such a simulation and computed the specter in Matlab, the result is qualitatively satisfying but I have a lot of noise, especially in the low frequency part of the specter. I have also tried some windowing without much luck. I simulated my sistem for about half a second, that's 25000 samples.

I must note that my simulation output gives out samples at an higher non constant rate, I've tried resampling my signal with no luck.

Given my limited experience with digital signal processing I am out of ideas. Is there something fundamentally wrong in my approach? Is it ok but there's something I should poke to get better results?

By definition, the Impulse Response of a Filter is the response the filter has to an Impulse: https://en.wikipedia.org/wiki/Impulse_response

Since your Cascaded Integrator Comb Filter is an optimized FIR filter, it's response will be Finite in length (as opposed to IIR - Infinite Impulse Response) and you can calculate this response by sending a single maximum value followed by zeros through it.

You mentioned [1, 0, 0, 0] is not an impulse response since 0 represents -1.  I believe the value that represents 0 would be half of the unsigned maximum value in this situation.  If your values are Unsigned integers from 0 to 255 with 0 representing the largest negative value (-1) and 255 representing the highest positive value (+1), you should put a maximum value (255) followed by zero values (128).  After the filter, you will probably shift/scale/truncate back to the previous bit-width to account for bit growth/gain of the filter.  This will give you the values of your Filter's Impulse Response.

You have a right idea for calculating the Impulse Response when you tried input '10101010...'  and '11010101010...' (these are really like putting in [1 -1 1 -1 ...] and [1 1 -1 1 -1 1 ...]).  Since FIR filters (CIC in this case) are linear systems, you can use the Superposition rule and add the inputs and then filter.  You will get the same answer as filtering the signals separately and then adding. That is why your experiment gives you twice the filter impulse response.  If you add after you filter, you will still get twice the filter impulse response.

I proved this to myself with this example (MATLAB):

If a is your filter:

a = [1 -1 1 1 1 -1]

and b and c are your inputs:

b = [1 -1 1 -1 1 -1 1 -1 1 -1 1]
c = [1 1 -1 1 -1 1 -1 1 -1 1 -1]

filter(a,1,b) + filter(a,1,c) = [2 -2 2 2 2 -2 0 0 0 0 0]

which is two times your filter impulse response.

However, you still get the filter impulse response by inputting a single max input:

filter(a,1,[1 0 0 0 0 0 0 0 0 0 0 ]) = [1 -1 1 1 1 -1 0 0 0 0 0]

For a very good treatment of CIC filter analysis, see:
http://dspguru.com/sites/dspguru/files/cic.pdf

To calculate the Frequency Magnitude Response of your filter, you can take the FFT of the coefficients.  To plot it, you should do something like this in MATLAB:

impResp = [1 -1 1 1 1 -1];

freqDomain = fft(impResp);
fSamp = 100;
frequencies = -fSamp/2:fSamp/(numel(impResp)-1):fSamp/2;

figure
plot(frequencies,db(abs(fftshift(freqDomain))))
xlabel('frequencies')
ylabel('db')

This filter does not have a very interesting Impulse response, though.


Using a Wiener Filter to Estimate a Transfer Function As a follow-on to this question about estimating a transfer function of an unknown system using a Wiener filter,

How would you put a minimum MSE criteria on how well the estimated filter weights matched the actual transfer function of the system? [Suppose you needed the MSE to be no more than -50dB]?
How would you change his formulation if you wanted poles as well as zeroes (an IIR rather than an FIR filter)?

The desired MSE is application dependent, so there can be no general
rule. If the approximation doesn't satisfy your needs you can
increase the filter length to obtain a better match.
There is no straightforward way to change the FIR Wiener filter solution to an IIR solution because the IIR formulation results in a set of nonlinear equations which have no closed-form solution. The IIR solution might also be unstable, so FIR filters are a much more practical choice when computing a Wiener filter.


Understanding sampling rate conversion I am a little bit confused about the concept of sampling rate conversion. Consider the following question:

Wideband telephony signals are sampled at a sampling frequency fS1 = 16kHz. To process these signals with an automatic speech recognizer, the sampling rate is converted here to fS2 = 11.025kHz, i.e., the wideband signal is ﬁrst upsampled by a factor U, lowpassﬁltered, and ﬁnally downsampled by a factor R. Determine the lowest possible values for R and U.

I am not sure whether U / R = 16kHz / 11.025 kHz, or vice-versa?

You would upsample by 441 and downsample by 640 since $11025 \cdot 640 = 16000 \cdot 441$

That would be a rather expensive operation, so in practice you would either fugde it a bit: 31/45 would get you to 11022Hz and 11/16 to 11000 Hz or you would use a suitable poly-phase irrational conversion algorithm


Relationship between FFT amplitude and sample size Would anyone know if there is any relationship between the FFT amplitude vs. data length?

For instance, in Matlab,
if signal_30k had 30,000 data points,
sampled at 1000 Hz,
and I do this:

fft_amplitudes = abs( fft( signal_30k ) );

then dominant frequency's resulting amplitude is some number.

but the same signal with the same frequency components,
signal_500k had 500,000 data points,
also sampled at 1000 Hz,
and had the same operation done on it:

fft_amplitudes = abs( fft( signal_500k ) );

the same dominant frequency's resulting amplitude would be another number but much larger.

Is there any relationship between the number of data points and the resulting amplitude that comes out of the FFT?
The goal is to ensure that regardless of the number of samples, the amplitudes have some means of being compared with one another.

Any help would be greatly appreciated.

If you want to have the same amplitudes, you need to scale the FFT result by the length of the FFT:
Fs = 100.0
f = 5;

t1 = np.arange(0, 5, 1/Fs)
t2 = np.arange(0, 50, 1/Fs)

s1 = np.sin(2*np.pi*f*t1)
s2 = np.sin(2*np.pi*f*t2)

plt.figure(figsize=(10,10))
plt.subplot(2,2,1)
plt.plot(t1, s1)

plt.subplot(2,2,2)
plt.plot(t2, s2)

S1 = np.fft.fft(s1) / len(t1)
f1 = np.linspace(0, Fs, len(t1))
S2 = np.fft.fft(s2) / len(t2)
f2 = np.linspace(0, Fs, len(t2))

plt.subplot(2,2,3)
plt.plot(f1, abs(S1))

plt.subplot(2,2,4)
plt.plot(f2, abs(S2))

If you want time and frequency domain signal to have the same power, you need to scale the FFT output by 1/sqrt(length).


1D Deconvolution with Gaussian Kernel (MATLAB) Suppose that I know the output and the transfer functions of a system and I would like to calculate the input function using deconvolution.

To get a grasp of the idea I have created a simple demonstration using Gaussians.

As a sanity check I first calculated the output function y using conv and then used deconv to replicate the original function. This works as expected. The input function calculated by deconvolution (green) is the same as the original input function (blue).

x = -2:0.1:2;
% the input function
a = 1; b=0; c=.5;
u = a*exp(-((x-b).^2)/(2*c^2));
figure(1)
plot(x, u), hold on,

% the transfer function
a = 0.1; b=0; c=.5;
h = a*exp(-(x-b).^2/(2*c^2));

% calculate the output function using conv
y = conv(h, u);
y2plot = conv(h, u,'same');
figure(1)
plot(x, y2plot, 'r')
legend('input func', 'output func')

u1 = deconv(y ,h);
figure
plot(x, u1,'g')

Fig 1.Input and Output functions

Fig 2.Transfer function

Fig 3.Input function calculated by deconv

Going back to the original problem, suppose that I know the output and transfer functions (the red and black curves in figures 1 and 2) and  I want to estimate the input function u.

x = -2:0.1:2;
% the transfer function
a = 0.1; b=0; c=.5;
h = a*exp(-(x-b).^2/(2*c^2));

% the output function
a = 0.9; b=0; c=.7;
y = a*exp(-(x-b).^2/(2*c^2));
figure;
plot(x, y, 'r')

u = deconv([zeros(1,20) y zeros(1,20)], h);
figure
plot(x, u,'g')

Fig 4.Input function calculated by deconv

The result I am getting is not what I would expect. Although I would expect a Gaussian similar to the green curve in Fig. 3 I am getting something different. Is this a matter of padding? Any ideas would be really appreciated.

Approaches

There are many methods for Deconvolution (Namely the degradation operator is linear and Time / Space Invariant) out there.
All of them try to deal with the fact the problem is Ill Poised in many cases.

Better methods are those which add some regularization to the model of the data to be restored.
It can be statistical models (Priors) or any knowledge.
For images, a good model is piece wise smooth or sparsity of the gradients.

But for the sake of the answer a simple parametric approach will be takes - -Minimizing the Least Squares Error between the restored data in the model to the measurements.

Model

The least squares model is simple.
The objective function as a function of the data is given by:

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

The optimization problem is given by:

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

Where $ x $ is the data to be restored, $ h $ is the Blurring Kernel (Gaussian in this case) and $ y $ is the set of given measurements.
The model assumes the measurements are given only for the valid part of the convolution. Namely if $ x \in \mathbb{R}^{n} $ and $ h \in \mathbb{R}^{k} $ then $ y \in \mathbb{R}^{m} $ where $ m = n - k + 1 $.

This is a linear operation in finite space hence can be written using a Matrix Form:

$$ \arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| H x - y \right\|_{2}^{2} $$

Where $ H \in \mathbb{R}^{m \times n} $ is the convolution matrix.

Solution

The Least Squares solution is given by:

$$ \hat{x} = \left( {H}^{T} H \right)^{-1} {H}^{T} y $$

As can be seen it requires a matrix inversion.
The ability to solve this adequately depends on the condition number of the operator $ {H}^{T} H $ which obeys $ \operatorname{cond} \left( H \right) = \sqrt{\operatorname{cond} \left( {H}^{T} H \right)} $.

Condition Number Analysis

What's behind this condition number?
One could answer it using Linear Algebra.
But a more intuitive, in my opinion, approach would be thinking of it in the Frequency Domain.

Basically the degradation operator attenuates energy of, generally, high frequency.
Now, since in frequency this is basically an element wise multiplication, one would say the easy way to invert it is element wise division by the inverse filter.
Well, it is what's done above.
The problem arises with cases the filter attenuates the energy practically into zero.  Then we have real problems...
This is basically what's the Condition Number tells us, how hard some frequencies were attenuated relative to others.

Above one could see the Condition Number (Using [dB] units) as a function of the Gaussian Filter STD parameter.
As expected, the higher the STD the worse the condition number as higher STD means stronger LPF (Values going down at the end are numerical issues).

Numerical Solution

Ensemble of Gaussian Blur Kernel was created.

The parameters are $ n = 300 $, $ k = 31 $ and $ m = 270 $.
The data is random and no noise were added.

In MATLAB the Linear System was solved using pinv() which uses SVD based Pseudo Inverse and the \ operator.

As one can see, using the SVD the solution is much less sensitive as expected.

Why is there an error?
Looking at a solution (For the highest STD):

As one could see the signal is restored very well except for the start and the end.
This is due to the use of Valid Convolution which tells us little on those samples.

Noise

If we added noise, things would look differently!
The reason results were good before is due to the fact MATLAB could handle the DR of the data and solve the equations even though they had large condition number.

But large condition number means the inverse filter amplify strongly (To reverse the strong attenuation) some frequencies.
When those contain noise it means the noise will be amplified and the restoration will be bad.

As one could see above, now the reconstruction won't work.

Summary

If one knows the Degradation Operator exactly and the SNR is very good, simple deconvolution methods will work.
The main issue of deconvolution is how hard the Degradation Operator attenuates frequencies.
The more it attenuates the more SNR is needed in order to restore (This is basically the idea behind Wiener Filter).
Frequencies which were set to zero can not be restored!

In practice, in order to have stable results on should add some priors.

The code is available at my StackExchange Signal Processing Q2969 GitHub Repository.

This answer was given to a different question yet both of them deals with deconvolution of a signal with a Gaussian Kernel.

If there are issues not answered, please point me and I will broaden.


Algorithm correcting chromatic aberration on hyperspectral images In this case, the picture below shows a color fringe detected on its edges.

I am experiencing the same difficulty on a hyperspectral image, so I am looking for a good algorithm correcting chromatic aberration for hyperspectral images. I understand there are several algorithms for RGB images, but it seems that it's difficult to extend any of them for hyperspectral images.

Having worked with hyperspectral images from experimental setups, I know it's hard to fix these little details.

The best way my group has to fix alignment and chromatic shifts is to use calibration images: Take a picture of something simple for which you know how the result should look like, and use the same correction on every pictures you take hereafter.

A simple example I'm thinking would be to take a picture of a small white dot against a black background. Every misaligned color will create a colored dot at a different place. You then find the mathematical operation necessary to put all these colored dots back together, creating a single white one (it will probably be a translation with a color dependant magnitude).

IF this shift is always the same from picture to picture, you can implement it in your data reduction code. IF not, you will have to take a calibration picture before every picture you take.

If you're lucky, it will be a translation. If you're not, the transformation might also depend on where you are in the picture. In that case, use a grid of white dots.


Spectral density interpolation: Applications? Let $\Phi_1$ and $\Phi_2$ be two matrix-valued power spectral densities. I wonder whether the problem of interpolating $\Phi_1$ and $\Phi_2$ has been studied in the literature and/or applied in practical scenarios.

In case of 1-D spectra, I know that a possible application could be speech/voice morphing, however I cannot imagine a practical application in which the considered spectra are matrix-valued functions.

Any suggestion/pointer to the literature is really appreciated!

  Addendum. The motivation behind my question comes from the fact that I'm currently working on the geometry of the space of power spectral densities (mostly from a mathematical viewpoint) and I would like to apply some of the mathematical tools I analyzed to some "practical" situations. One of these tools is a method for interpolating spectral densities. I'm especially interested in the case in which spectral densities are functions that take matrix values, since in this case the geometry of the space of spectral densities becomes particularly "rich" and "interesting".

TL;DR:
The basis of this response is that matrix valued PSDs imply vector valued quantities that are expressed as Tensors. The Tensor is explored, briefly, to provide a basis to think (and possibly search further, beyond this post, for) practical applications. Practical applications explored here include DWTI-MRI and Quaternion RGBs.

  Let $\Phi_1$ and $\Phi_2$ be two matrix-valued power spectral densities.

How does a Power Spectral Density ends up being matrix-valued?

By applying it over a Tensor field or other sequence of Tensors.

What is a Tensor? And why do I need one?

The answer to "What is a Tensor" varies wildly depending on who one dares to ask. That's because different specialities are interested in different aspects and applications of Tensors. From my point of view, I thought I would provide a "necessity" based 2 min definition to cover how do we get to that "matrix valued" part:

A "scalar field" type of signal, associates a single quantity ($x \in \mathbb{R}^d$ ) with some point in space. That is, independently of the dimensionality ($d$) of the scalar:

One dimensional signals such as a discrete sampled sound recording ($x[i]$) are represented as a sequence of real ($x \in \mathbb{R}$) samples of sound pressure at distinct points in time. In this case, all we are interested in is how does the Sound Pressure Level vary over time. Notice here that $x$ is assumed to be one continuous recording without "gaps". In this case, successive samples $i, i+1$ are assumed to have been obtained at some $T_s$ time interval called "Sampling Period" (and of course, $i$ is an integer ($i \in \mathbb{Z}$)) . That is, we don't need to convey that $x[i]$ contains successive samples. It is implied.
Two dimensional signals, such as grayscale images, are still scalars. In visible photography, a sensor (or film surface) records how does light intensity vary within its Field of View. Again, we associate a single measurement ($x[i,j]$) with a point in space ($i,j$).
Three dimensional signals, such as volumetric data coming out of structural Magnetic Resonance Imaging scanner, are still scalars. In this case, all we are interested in is the concentration of a substance ($x[i,j,k]$) in three dimensional space ($i,j,k$). This concentration implies that $x$ might be more dense at $i_1,j_1,k_1$ but less dense around $i_2,j_2,k_2$ as it happens with fatty tissue around muscles for example.

Three dimensional signals such as video are still scalars. That is, we associate light intensity $x[i,j,n]$ with a point in 2D space $i,j$ and a point in time $n$.

Four dimensional signals such as those coming out of functional Magnetic Resonance Imaging are still scalars. In this case, we associate the concentration of a substance $x[i,j,k,n]$ with a point in 3D space $i,j,k$ and a point in time $n$.

But, what do scalar fields look like and why are they not enough?

This is what a three-dimensional scalar field looks like:

And here is a more realistic one. I am showing this one because the effect is more pronounced.

What do we observe? Well, it's not smooth is it? Why is it not smooth? Because it is only composed by little cubes that represent sampled space without any regard for direction. I only know that a particular point in space had some value but I don't know where was this value coming from.

Contrast this previous image with something like this:

This skull looks like it is composed of tiny little surfaces that tell us not only that some point in space $i,j,k$ is occupied by some bone substance but also what the orientation of this substance, in space, is.

To be able to represent this at the point of measurement, we cannot use a scalar field. We need something different. We need something that can encode direction. We need a vector.

Now, someone might say here: Hold on, you can calculate the local slope of a scalar field by obtaining its Gradient. Which, by the way, is a vector field.

And, for some applications, that could be enough but that would also mean that we make a particular assumption about the underlying phenomenon. That is, that "things" flow from high pressures to low pressures. In other words, we would still sample a scalar field of pressure (or density) values and from that infer flow. What if these regions of "high pressure" don't move? (in reality). Or what if only a small amount of their "density" flows?

What if what we need, genuinely, is to be able to not only record the concentration of a substance $x$ at some point in three dimensional space ($i,j,k$) but also its flow towards some direction $\rho, \theta, \alpha$?

Such is the case with Diffusion Tensor Magnetic Resonane Imaging. MRI is a beautiful "sensor" that can return a huge amount of information about a biological sample in its field of view.

It samples the human body directly in the frequency domain and to produce an "image" it transforms the acquired signals, from the frequency domain to the time domain.

Because of the way magnetic fields interact with magnetic objects (even as tiny as atoms), it is possible to selectively "excite" and "listen" for signals at a particular point in three dimensional space $i,j,k$.

This is how structural and functional MRI works which return scalar fields. You "excite" a particular region and then listen to the signal of atoms (from that particular region) "reacting" to the excitation (typically in the region of MHz, e.g. for H the frequency is ~43 MHz per Tesla). The more atoms you "listen" the higher the $x[i,j,k]$ at that point. Or more generally, $x$ is proportional (even if negatively) to the concentration of specific atoms at that $i,j,k$

But, by modifying the way atoms are excited you can also infer if and how did they move. How do you do that? You "excite" a particular region but "listen" from another region in quick succession. In other words, you "tag" some atoms, wait and then look where did the tag go (but very quickly because the "tag" is "wearing out" with time).

The geometrical orientation of these "excite" / "listen" regions now forms a vector. If you found "a lot of tagged atoms" towards a particular direction then this direction has higher magnitude of "flow", otherwise it does not.

So, in this case, we need to represent both the magnitude of the quantity ("How much substance moved?") and its direction ("In what direction did it move?").

We do this using a Tensor Field. It looks like this:

$$A[i,j,k] = \begin{bmatrix}
 A_{ii} & A_{ij} & A_{ik} \\
 A_{ji} & A_{jj} & A_{jk} \\
 A_{ki} & A_{kj} & A_{kk} \\
  \end{bmatrix}
$$

So, imagine a three dimensional grid such as this one:

This grid has a "node" (a point where three axes meet) at distinct locations $i,j,k$. But this time around, instead of having a single number associated with that $i,j,k$ we have a Tensor ($A$).

The MRI populates this matrix by obtaining six sets of measurements corresponding to "geometrical displacements" (i.e. directions) $A_{ij}, A_{ik}, A_{ji}, A_{jk}, A_{ki}, A_{kj}$. Pairs of $indexVar_1 indexVar_2$ denote directions. So $ij$ can be left to right, $ji$ right to left, $ik$ top to bottom, $ki$ bottom to top, $jk$ front to back $kj$ back to front WITH RESPECT TO $[i,j,k]$.

It is important to understand that $A$ is a matrix and it is one thing. Just as some vector is a "packet" (tuple) of two numbers ($\vec{a} = (b,c))$ but we treat it effectively as one thing, in exactly the same way, $A$ is one thing. One, indivisible, inseparable object. ONE thing.

This is important for what is about to follow.

Now, if you check back our $A$, you will notice that it contains a huge amount of information. So much in fact, that we can infer flow towards any direction.ANY direction. No, seriously, look back at $A$ again. $A$ (and of course, not $A$ but the MRI itself) captures "symmetric" flows. It captures flow from left to right and at the same time it captures a flow from right to left. Therefore, we can "see" exchange phenomena. Something travelling from left to right (perhaps mostly) but at the same time another flow (perhaps a low amount) travelling backwards in the same direction at the same site.

Most of the times however what we are interested in is just one vector out of all the possible directions that $A$ can encode (The vector of "most flow"). We extract this via the eigenvectors and eigenvalues of $A$.

The last step in this process is to visualise the flow. In this case, we simply "follow the arrows" to reconstruct an orbit. This orbit coincides with neuronal axons (in the brain) and this is what the DTI images depict. About this image, it is important to understand that this is a static image. This is a snapshot of the tracts, although it has been created by an extremely large amount of measurements that may have taken some time. If we did DTI over time, then we would have one $A$ associated with each $i,j,k,n$.

OK, that's fine, but how do you end up with PSDs that look like matrices?

Remember that $A$ that is one thing? Just like some $x[i]$ (e.g. some integer number) is one inseparable thing?

Find the periodogram of $A[i,j,k]$.

This is where it pays off to think about $A$ as one thing. Because you don't treat it differently than any other vector, although, yes, it certainly is not "just a vector".

Now, I will use Welch's method of deriving the periodogram here which specifies splitting the "signal" into overlapping segments, taking the DFT of those segments, deriving the spectrum and then averaging those segments to derive the final power spectrum.

BUT!!!

The most challenging thing here is NOT the computational part. The most challenging thing here is the connection with reality. What does it mean? What does it look like? In fact, the most challenging thing in the aforementioned process of deriving the periodogram is not the "indexing" (chopping the Tensor field down to overlapping segments) because that's just like slicing a scalar array. Only that now, each one of the entries of this Tensor field is a Tensor. So, this would look like an Array of Tensors, which we cut down into overlapping segments. It is not the averaging, because, averaging is averaging. Tensor addition works as you would expect, you just add the individual components.

The most challenging part, in terms of comprehension is The Discrete Fourier Transform. Again, it's not about calculating it, because, the Discrete Fourier Transform of a Tensor field is the Discrete Fourier Transform of the individual components of the Tensors that compose it. In terms of notation, you still have sinusoids but they are packed in a Tensor "package" this time.

The challenging bit here is the meaning of it and what does it look like given that it describes so much information.

One (useful) view would be to treat it like the DFT of complex data. Because in that case, it's not $x \in \mathbb{R}$ anymore but rather $x \in \mathbb{C}$ and therefore $x$ is now actually a vector (i.e., it has magnitude and direction).

When $x \in \mathbb{C}$ then what the Fourier transform actually does is decompose $x$ to the sum of a set of "squiggly" Lissajous symbols. That is, the complete repertoire of "doodles" that two interacting (orthogonal) sinusoids can create, including a perfect circle (at DC) but more importantly ellipsoids at various orientations.

This concept can be applied as is to the Fourier Transform of the Eigenvector of $A$ described above. In that case what we would be getting after $\mathcal{F}(\lambda \cdot v), A = \lambda \cdot v$ would be a decomposition of $v$ down to a set of "squiggly" three dimensional surfaces including a perfect sphere and also ellipsoids with various orientations (and everything in between).

I suppose that something similar can be obtained to interpret $\mathcal{F}(A)$ although it probably does not look exactly the same because each frequency of the PSD is now a complex matrix and returns a huge amount of directional information about that particular point (which we kind of "discard" when we focus only on the eigenvector).

  ...a practical application in which the considered spectra are matrix-valued functions

This has been partially covered by the DTI and discussion in the previous section. The only thing that I would like to add is that in terms of interpolating, you can interpolate incredibly complex curves in surprising ways with results similar to this one.

And also that in the signal / image processing space, there is a "Tensor" expression of the DFT that would return spectra you can work with but I am not really sure about interpolation.

Perhaps you can define some sort of "smoothing" filter. That is, given a Tensor field with missing spectral values, fill them in by interpolating the known values around it.

Other than this, I suppose that you can define a metric of "distance" between two spectra (?) to "match" tensor fields.

Beyond these applications in DSP, the next stop would be the Physics SE. They use Tensors a lot in cosmology. I don't quite understand this stuff,...I mean, what does it look like...and sometimes the "why". I guess they like tensors because they express that oneness required to describe space and time...but how do you go from that to, "maybe "time" does not exist", is impenetrable for me.

Anyway, this is why I asked "what is it that you do"? What is your "necessity" for ending up with Matrix valued PSDs? What is the "problem"?

Hope this helps.


Where do RBJ's analog shelving filters come from? In RBJ's Audio EQ Cookbook, the analog prototype for a low shelf filter is given as such:
$$H(s)=A\cdot\frac{s^2+\left(\frac{\sqrt{A}}{Q}\right)\cdot s + A}{A\cdot s^2 + \left(\frac{\sqrt{A}}{Q}\right)\cdot s + 1}$$
Where does this prototype come from? From what I understand, the lowpass and highpass filters are standard Butterworth filters. But how was this particular transfer function derived? Is it some transformation of a Butterworth filter?

If you remove (for the time being) that leading factor $A$ as a constant gain factor:
$$H(s)=\frac{s^2+\left(\frac{\sqrt{A}}{Q}\right)s + A}{As^2 + \left(\frac{\sqrt{A}}{Q}\right)s + 1}$$
what you get then is a symmetric, but otherwise general shelf that could be equally described as "LowShelf" or "HighShelf".  In dB, the gain at the low end is the negative of the dB gain at the high end.  And, in log gain (dB) vs. log frequency (octaves or decades), the gain function always has odd symmetry about the "resonant frequency" (more accurately described as the "significant frequency" which is the shelf mid-point frequency, a little different from the biquad resonant frequency), which in the normalized case, is $\Omega_0=1$.  At that shelf mid-point frequency this general low/high shelf has a gain of $1$ (or 0 dB).
$$\begin{align}
\Big|H(s)\Big|^2 &= \left| \frac{s^2+\left(\frac{\sqrt{A}}{Q}\right)s + A}{As^2 + \left(\frac{\sqrt{A}}{Q}\right)s + 1} \right|^2 \\
\\
\Big|H(j\Omega)\Big|^2 &= \left| \frac{(j\Omega)^2+\left(\frac{\sqrt{A}}{Q}\right)(j\Omega) + A}{A(j\Omega)^2 + \left(\frac{\sqrt{A}}{Q}\right)(j\Omega) + 1} \right|^2 \\
\\
 &= \left| \frac{A-\Omega^2+j\frac{\sqrt{A}}{Q}\Omega}{1-A\Omega^2+j\frac{\sqrt{A}}{Q}\Omega} \right|^2 \\
\\
 &= \frac{\Big|A-\Omega^2+j\frac{\sqrt{A}}{Q}\Omega \Big|^2}{\Big|1-A\Omega^2+j\frac{\sqrt{A}}{Q}\Omega\Big|^2}  \\
\\
 &= \frac{(A-\Omega^2)^2+\left(\frac{\sqrt{A}}{Q}\Omega\right)^2}{(1-A\Omega^2)^2+\left(\frac{\sqrt{A}}{Q}\Omega\right)^2}  \\
\\
 &= \frac{ A^2-2A\Omega^2+\Omega^4 + \frac{A}{Q^2} \Omega^2}{1-2A\Omega^2+A^2\Omega^4 + \frac{A}{Q^2} \Omega^2}  \\
\\
 &= \frac{ A+\frac{1}{A}\Omega^4 - (2 - \frac{1}{Q^2}) \Omega^2}{\frac{1}{A}+A\Omega^4 - (2 -\frac{1}{Q^2}) \Omega^2}  \\
\\
 &= \frac{ \frac{A}{\Omega^2}+\frac{\Omega^2}{A} - (2 - \frac{1}{Q^2})}{\frac{1}{A\Omega^2}+A\Omega^2 - (2 -\frac{1}{Q^2})}  \\
\end{align}$$
Remember that the "significant frequency" is set to one.  If you express the normalized frequency in terms of log frequency, then
$$ \Omega \triangleq 2^p $$
where $p$ is the log frequency, in octaves, centered at the "significant frequency".
$$\begin{align}
 \Big|H(j2^p)\Big|^2 &= \frac{ \frac{A}{2^{2p}}+\frac{2^{2p}}{A} - (2 - \frac{1}{Q^2})}{\frac{1}{A2^{2p}}+A2^{2p} - (2 -\frac{1}{Q^2})}  \\
\\
 &= \frac{2^{2p-\log_2(A)} + 2^{-2p+\log_2(A)} - (2 - \frac{1}{Q^2})}{2^{2p+\log_2(A)} + 2^{-2p-\log_2(A)} - (2 - \frac{1}{Q^2})} \\
\\
 &= \frac{2 \cosh\big(\log(2)(2p-\log_2(A))\big) - (2 - \frac{1}{Q^2})}{2 \cosh\big(\log(2)(2p+\log_2(A))\big) - (2 - \frac{1}{Q^2})} \\
\\
 &= \frac{\cosh\big(2\log(2)p-\log(A)\big) - (1 - \frac{1}{2Q^2})}{\cosh\big(2\log(2)p+\log(A)\big) - (1 - \frac{1}{2Q^2})} \\
\end{align}$$
Now if we express this as log amplitude (let's say in dB) vs. log frequency (in octaves w.r.t. the "significant frequency"), this becomes
$$\begin{align}
\tfrac{10}{\log(10)} \log \left(\Big|H(j2^p)\Big|^2 \right) &= \tfrac{10}{\log(10)} \log \left(\frac{\cosh\big(2\log(2)p-\log(A)\big) - (1-\frac{1}{2Q^2})}{\cosh\big(2\log(2)p+\log(A)\big) - (1-\frac{1}{2Q^2})}\right) \\
\\
&= \tfrac{10}{\log(10)} \log \left(\cosh\big(2\log(2)p-\log(A)\big) - (1-\tfrac{1}{2Q^2})\right) \\
& \qquad - \tfrac{10}{\log(10)} \log \left(\cosh\big(2\log(2)p+\log(A)\big) - (1-\tfrac{1}{2Q^2})\right) \\
\end{align}$$
Since $\cosh(\cdot)$ is an even-symmetry function, I think you can show that the function above is an odd-symmetry function in $p$.  Substitute $-p$ in for $p$ and you will see that the + and - terms swap.
$$\begin{align}
\frac{\mathrm{d}}{\mathrm{d}p} \tfrac{10}{\log(10)} \log \left(\Big|H(j2^p)\Big|^2 \right) &= \tfrac{10}{\log(10)} \frac{\mathrm{d}}{\mathrm{d}p} \log \left(\cosh\big(2\log(2)p-\log(A)\big) - (1-\tfrac{1}{2Q^2})\right) \\
& \qquad - \tfrac{10}{\log(10)} \frac{\mathrm{d}}{\mathrm{d}p} \log \left(\cosh\big(2\log(2)p+\log(A)\big) - (1-\tfrac{1}{2Q^2})\right) \\
\\
&= \tfrac{10}{\log(10)}  \frac{2\log(2)\sinh\big(2\log(2)p-\log(A)\big)}{ \cosh\big(2\log(2)p-\log(A)\big) - (1-\tfrac{1}{2Q^2})} \\
& \qquad - \tfrac{10}{\log(10)} \frac{2\log(2)\sinh\big(2\log(2)p+\log(A)\big)}{ \cosh\big(2\log(2)p+\log(A)\big) - (1-\tfrac{1}{2Q^2})} \\
\end{align}$$
When $p=0$ (at the midshelf) this shelf slope is:
$$\begin{align}
\frac{\mathrm{d}}{\mathrm{d}p} \tfrac{10}{\log(10)} \log \left(\Big|H(j2^p)\Big|^2 \right) \Bigg|_{p=0} &= \tfrac{10}{\log(10)}  \frac{2\log(2)\sinh(-\log(A))}{ \cosh(-\log(A)) - (1-\tfrac{1}{2Q^2})} \\
& \qquad - \tfrac{10}{\log(10)} \frac{2\log(2)\sinh(\log(A))}{\cosh(\log(A)) - (1-\tfrac{1}{2Q^2})} \\
\\
&= - \tfrac{10}{\log(10)} \frac{4\log(2)\sinh(\log(A))}{ \cosh(\log(A)) - (1-\tfrac{1}{2Q^2})} \\
\\
&= - \tfrac{10}{\log(10)} \frac{2\log(2) \left(A - \tfrac{1}{A} \right) }{ \tfrac{1}{2} \left(A + \tfrac{1}{A} \right) - (1-\tfrac{1}{2Q^2})} \\
\\
&= - \tfrac{40\log(2)}{\log(10)} \ \frac{A^2 - 1}{A^2 + 1 - 2A(1-\tfrac{1}{2Q^2})} \\
\\
&= - \tfrac{40\log(2)}{\log(10)} \ \frac{A^2 - 1}{A^2 - 2A + 1 + \tfrac{A}{Q^2}} \\
\\
&= - \tfrac{40\log(2)}{\log(10)} \ \frac{(A-1)(A+1)}{(A-1)^2 + \tfrac{A}{Q^2}} \\
\\
&= - \tfrac{40\log(2)}{\log(10)} \ \frac{(A-1)(A+1)}{(A-1)^2 + \tfrac{A}{Q^2}} \\
\end{align} $$
I'm still not done with this.  I'm doing this from memory.  I can't find my notes from a quarter century ago anywhere.


Using the linear operator to check for time invariance of a differential equation? I have a differential equation $$\frac{d^2y(t)}{dt^2}+y(t) = \frac{dx(t)}{dt} + x^2(t)$$ and I need to see if this system with input $x(t)$ and output $y(t)$ is time invariant and linear. I tried to break up the equation by substituting a linear operator D as the derivative.

So I get: $$y(t) = \frac{Dx(t)+x^2(t)}{D^2 + 1}$$

But the problem is, this equation is meaningless since you are dividing by the derivative operator. However, mathematically, it works out, so can I get $y(t)$ in this form and test for linearity or time invariance using this method?

Let $y(t)=\mathcal{H}(x(t))$ describe the system by relating its input $x(t)$ to the output $y(t)$. Then, to check linearity simply check with $x(t)=1$ (i.e. a constant). You will see, that

$$\mathcal{H}(x(t)=1)+\mathcal{H}(x(t)=1)\neq \mathcal{H}(x(t)=1+1)$$

which is one requirement for the linearity. So, the system is non-linear. Actually, you can immediately see this from the fact that $x(t)$ counts quadratically, which is a non-linear relation.

Regarding time-invariance, you see that the time-argument is the same for all functions and the time-value is not used in the equation itself (such as $t\cdot x(t)$ would ne not time-invariant). So, the system is time-invariant.


Calculating the total energy of a signal I have a vector of signals and need to calculate the total energy. I have read somewhere that you can do this by summing up the entire vector and then finding the absolute value of the result. Is this correct? Could anyone suggest any tutorials, or, equations if this is the wrong method, thanks.

The energy of a discrete-time signal is the sum of the squares:

$$
E_x = \sum_n |x[n]|^2
$$

If you need the average power of your signal, divide the energy by the number of samples in your vector. Power is the average energy per sample of your signal.

If your signal vector is a sample of a stationary process (i.e. a recording clip of someone talking), then average power probably has more meaning than energy (which will grow with the sample length).


Suggest me some good references for sonification I need good reference books/tutorials/doccumentations on sonification.

Please suggest me some books or other references from which I can learn about sonification and algorithms of sonification.

Thanks.

1)

Proceedings of annual International Conference on Auditory Display. The annual ICAD conference is the prime forum for academia and industry to exchange ideas and discuss developments in the field of auditory display.

https://smartech.gatech.edu/handle/1853/49750

2)

http://sonification.de/handbook/

This book is a comprehensive introductory presentation of the key research areas in the interdisciplinary fields of sonification and auditory display. Chapters are written by leading experts, providing a wide-range coverage of the central issues, and can be read from start to finish, or dipped into as required (like a smorgasbord menu).


disadvantages of FFT, it can not extract enough frequencies without enough samples Let's say sampling rate is $Fs = 44\mathtt{kHz}$, now I have $N = 2048$ samples, then I can get $N/2 + 1 = 1025$ frequencies.

I'm confused by Matlab's FFT documentation that says the frequencies are present at:

(Fs/2) * linspace(0, 1, N/2+1)

It is very probable that the real frequencies are not among  the frequencies generated. What should we do to get the real frequencies?

Any measurement has a finite precision. You can't find real frequency exactly in general. You can evaluate it with necessary precision instead. In the case of FFT if the real frequency doesn't match FFT frequency grid exactly, power spectrum is spreading mostly among neighbouring frequency bins (among all bins in general, but most of it is focused in a few adjacent ones). Nevertheless you can estimate the real frequency by finding maximum of such a spectrum. Dealing with FFT you can improve the precision of frequency search by increasing the number of samples. If you can't, you should use another power spectrum estimation techniques but it's the separate question.

Hope this helps.


What is importance of the conjugate reciprocal and none reciprocal zeros in FIR Type 2 Filters? In DSP class we have this slide:

So i like to know why this classification is so imoortant?

I have seen this post :

  Conjugate reciprocal pairs of zeros and poles in ...

  Assuming the impulse response h[n]h[n] of an FIR filter is real for
  all nn,

  Why are zeros and poles in FIR design found in reciprocal and
  conjugate pairs?Is the assumption necessary for this phenomenon to
  take place?

But i dont get my answer?( mybe because of being newbi in this field)

I would like to provide a short answer.

Consider

$$ h[n] = \sum_{k=0}^{M} b_k \delta[n-k] ~~~~\longleftrightarrow ~~~~H(z) =  \sum_{k=0}^{M} b_k z^{-k} $$

Left: LTI FIR filter impulse response coefficients $b_k$, and right: its Z-Transform.

If the FIR filter is real; i.e., coefficients $b_k$ are real, then the Z-transform is a polynomial of real coefficients, and the roots of such polynomials are either real complex-conjugate pairs. This explains why a real FIR filter's Z-domain poles & zeros are in conjugate pairs when they are complex.

To explain why, also, they are in reciprocals, we need to consider linear phase property (without losing generality) : $h[n] = h[-n]$; then since the Z-transform of $h[-n]$ is $H(1/z)$; we have the following :

$$ h[n] = h[-n]  ~~~\implies ~~~ H(z) = H(1/z) $$

And this means that if $z_0$ is a zero (or pole) of $H(z)$; i.e., $H(z_0) = 0$ then its reciprocal $1/z_0$ is also a zero of $H(z)$.

To sum up; for a real and linear phase FIR filter, if $z_0$ is a zero of $H(z)$ then $z_0^*$, $1/z_0$ and $1/z_0^*$ are also zeros; hence zeros and poles are in pairs of four when they are complex.


Manipulating with frequencies through fftw I'm using portaudio for I/O audio streaming, and I wanted to apply some equalization on audio input from microphone. I found that most people refer to fft for task like that, so I use FFTW library in my project. My first idea was to select frequency bins in some range, from fft resulting complex array, and  increase their amplitudes by scaling them with some factors. But this resulting only in noise,  there is no sign of any equalization. So my question is, is it possible to increase volume at specified frequencies, through fft, and if so, how to do that. My experience so far makes me think that any change of those values in fft complex array just destroys the sound. Is there a proper way to manipulate with that complex data from fft resulting array?
Any tips for equalization other than fft are also welcome.

You applied the FFT to the full signal? Then modification of the amplitudes will likely destroy any time locality.

Use the so-called STFT, short term Fourier transform, or overlap-add method to get a time-localized frequency analysis.


Noise PSD and sampling rate relation Let's consider generating samples of a random process like white Gaussian noise (AWGN). Let's assume I am generating $N$ samples of AWGN with variance $\sigma ^2$ in MATLAB by using randn() funtion i.e.:

\begin{equation} \tag{1}\label{Eq:n}
    n_1 = \sigma * \text{randn}(N,1) \\
    n_2 = \sigma * \text{randn}(N,1)
\end{equation}

Note that, here I am not sampling a noise signal to get these samples. I am simply generating some random numbers to consume them later.

Now I want to add these noise samples to my clean signal samples. Let's assume I have two cases of clean signal samples, namely $y$ which has a sampling frequency of $f_s$ and $z$ which has a sampling frequency of $2f_s$. So we have

\begin{equation} \tag{2}\label{Eq:y_n}
    y_n = y + n_1
\end{equation}

\begin{equation} \tag{3}\label{Eq:z_n}
    z_n = z + n_2
\end{equation}

This means that in $y_n$ I am using noise samples with the rate $f_s$ and in $z_n$ and I am using noise samples with the rate $2f_s$.

Since here we are dealing with random process, instead of frequency spectrum we have the power spectral density (PSD). Let's name the PSD of $y_n$ and $z_n$, $S_{yy}(\omega)$ and $S_{zz}(\omega)$ respectively. Also let's consider the PSD of noise in Eq. \eqref{Eq:y_n} which is read at the lower rate is $S_{nn1}(\omega)$ and the PSD of noise in Eq. \eqref{Eq:z_n} which is read at the higher rate is $S_{nn2}(\omega)$

Question 1) How are $S_{nn1}(\omega)$ and $S_{nn2}(\omega)$ different? More precisely, if you want to draw them on the frequency axis how do they look? Does $S_{nn1}(\omega)$ span from 0 to $f_s$ and  $S_{nn2}(\omega)$ span from 0 to $2f_s$?
And if you calculate the area underneath them (which shows the power or variance of the noise) which one is greater? Do they have the same power?

Let's consider the auto-correlation and PSD for random processes. We know that PSD is the DTFT of the auto-correlation, more specifically:

\begin{equation} \tag{4}\label{Eq:DTFT}
    R_{nn}(m) = \frac{\Delta t}{2 \pi}\int_{0}^{\frac{2\pi}{\Delta t}} S_{nn}(\omega) e^{j\omega m \Delta t} d \omega,
\end{equation}
where $R_{nn}$ is the auto-correlation and $S_{nn}(\omega)$ is the PSD.

From Eq. \eqref{Eq:DTFT} we have:
\begin{equation} \tag{5}\label{Eq:Rnn}
    R_{nn}(0) = \frac{\Delta t}{2 \pi}\int_{0}^{\frac{2\pi}{\Delta t}} S_{nn}(\omega) d \omega,
\end{equation}
which means the area under PSD (which is the power of the signal) is given by the auto-correlation at zero.

Question 2) Which one gives the power of the signal, just the output of the integral (without the multiplier) in Eq. \eqref{Eq:Rnn}, or the whole Eq. \eqref{Eq:Rnn}?

Question 3) This is the continue of question 1, or maybe similar to it. I basically want to know how Eq. \eqref{Eq:Rnn} changes for noises $n_1$ and $n_2$? Cause obviously in Eq. \eqref{Eq:Rnn}, $\Delta t$ is changing for the noise when the sampling rate is changed. But how $R_{nn1}$ and $R_{nn2}$ or $S_{nn1}$ and $S_{nn2}$ change when sampling rate changes?

I appreciate if someone can mathematically explain these to me.

Important Information : Sampling at $f_{s}$ will map $[-\frac {f_{s}}{2}, \frac {f_{s}}{2}]$ to digital frequency $\omega=[-\pi, \pi]$, and similarly sampling at $2f_{s}$ will map $[-f_{s}, f_{s}]$ to digital frequency $\omega=[-\pi, \pi]$. Also, we need to look only at digital frequency $\omega = [-\pi, \pi]$ as the digital spectrum is a $2\pi$-periodic function.

Now, lets look at discrete White Gaussian Noise Samples which you have generated.
$$n_{1} = N(0,\sigma^{2}), n_{2} = N(0, \sigma^{2})$$

This means auto-correlation function for both $n_1$ and $n_2$ will be:
$$r[k] = \sigma^{2}\delta[k]$$
and hence, the PSD of $n_1$ and $n_2$ in digital domain will be :
$$P(e^{j\omega}) = \sigma^2, \omega \in [-\pi, \pi]$$
and hence total Noise Power will be :
$$T_{NoisePower} = 2\pi.\sigma^2$$

Basically, once you have discrete Noise Samples for a White Gaussian Process, then autocorrelation function is fixed and hence the PSD in digital domain is fixed too. Sampling frequency does not matter for determining Noise Power.

Though, oversampling can help in increasing SNR of the received Signal. Think about it!

Hint: Oversampling followed by filtering!!


Why $H(A)=H(C)$ where $C$ is $A$ with an additional parity bit? Let $A=\{00,01,10,11\}$ with equal probabilities for each symbol, and $B=\{0, 1\}$ be a parity generator such that
$$
b=\begin{cases}
0, & \text{if} \,\, a=00 \quad \text{or} \quad a=11 \\
1, & \text{if} \,\, a=01 \quad \text{or} \quad a=10
\end{cases}
$$
Now assume we transmit $(a_i,bj)$ where $a_i$ is a symbol in $A$ and $b_j$ is the parity bit associated with it. I calculated the entropy for $A$ as follows:
$$
H(A)=4 [4 \log_2(4)]=2
$$
and to calculate the entropy of the new symbols, call it alphabet $C=\{000,011,101,110\}$, we need the probabilities:
$$
p_c(0)=\mathbb{P}(a=00 \, \text{and} \, b=0)=\mathbb{P}(a=00)\mathbb{P}(b=0 \mid a=00)=\mathbb{P}(a=00)=1/4
$$
Similarly, $p_c(1)=p_c(2)=p_c(3)=p_c(0)=1/4$, the entropy of $C$ is
$$
H(C)=4 [4 \log_2(4)]=2=H(A)
$$
How is this the case? We have more bits per symbols in $C$.

Because the entropy represents information quantity, or if being measured in bit, the smallest number of bits per symbol we need to represent a source.

The source $A$ contains $4$ equiprobable symbols hence it is obvious that we need $\log_2(4) = 2$ bits per symbol to represent the source.

The source $C$ is simply $A$ with its parity bits hence no new information is added. We have more bits per symbols, but the number of information bits is unchanged.

In communications, these additional bits are added to scope with error of transmissions.


Usage of Digital Cosine Transform in lossy compression of images I was reading about Digital Cosine Transform (DCT Type-I) and it's use in JPEG.
I have the following doubt. My teacher offered an explanation which wasn't very convincing.
What property of Digital Cosine Transform makes it attractive for compression of images? Why Digital Sine Tranform is not used?
More specifically, what is the intuition behind the usage of DCTs?

A rectangular window of data can have a discontinuity between the front and back edge.  An FFT of that window represents that circular edge discontinuity with energy in a bunch of high frequency bins.  Take that segment of signal and mirror it around one edge.  Note that a circular or periodic extension of this no longer has a discontinuity across the edge.  The FFT of that doubled-by-mirroring signal might thus have less information energy in the high bins.  Thus there may be more of these high bins that can be "thrown away" by a compression algorithm without major damage (information loss).

A DCT is an FFT of the doubled-by-mirroring product.  Since the mirrored intermediate vector is symmetric, some efficiency can be gotten by throwing away any computations only required for the sine components since they will be zero anyway.

A DST would be the FFT of a doubled-by-mirroring-the-inverse signal, which would make the discontinuities worse, which would spread even more information/energy into the bins you would really like to throw away, thus making things far worse for simply but lossy compression.


Could you describe the effects for varying different parameters of a canny edge detector? Last couple of questions touched upon Canny edge detector

What are the limitations of a Canny edge detector?  and
Best way of segmenting veins in leaves?

The basic outline of the algorithm is as follows:

a. Apply Gaussian Convolution. (Choice of $\sigma$ to be made here)
b. Apply 2D derivative
c. Tracking through the ridges of this edges and thresholding (set pixels to zero which are not on the edge) with Hysteresis Lower and Higher T0 and T1 (Choice of $T0$ and $T1$ to be made here).

Read this for more background.

While, it is claimed that Canny is optimal; when getting practical results matters, tweaking factors  as listed above $\sigma, T0,$ and $T1$ does make a lot of difference.

So how does one select these (tweaking) parameters practically? Even if, there is no definite approach or value, what is the general technique to know this?

Following http://www.kerrywong.com/2009/05/07/canny-edge-detection-auto-thresholding/ is one of the few resources that shows how to choose thresholds Tlow and Thigh

According to this, for a picture which is sufficiently spread in historgram, one can choose T_low = 0.66 * mean value of image and T_high = 1.33 * mean value.

However, when the image is not sufficiently spread, one should use median as opposed to mean value of the image.

If the gap between T_low and T_high is very small, the resultant edges will be smaller in continuity and hence there will be more fractions. As the gap increases, you will have more single line edges.

As regards to sigma, as the sigma increases, the smoothing increases and noisy edges will go away, but at the same time, the location of edges may also move little. See this document, http://www.cse.unr.edu/~bebis/CS791E/Notes/EdgeDetection.pdf and results at page 29 shows this effect.


Filtering vs. decimation I have a signal collected via an RTL-SDR at 2.4Msps.  I want to narrow in on a small signal within this (an FM broadcast station, for example).  I am multiplying my signal by a complex exponential to shift to the target frequency.  Now, I'd like to narrow in.

I see two ways of doing this.  The first is to decimate the signal down to the target signal's bandwidth.  The second is to filter.  In my experiments filtering takes much longer to perform.  Is there any reason to prefer one technique over the other?

Let us define the terms first:

Downsampling means reducing the sampling rate of a signal. If there is energy outside the new Nyquist frequency, there will be aliasing.
Filtering means eliminating all energy in some frequency bands.
Decimation means anti-alias filtering a signal and then downsampling it.

In a typical SDR application, your captured bandwidth is larger than the actual band of interest. There are two things you want to do:

Eliminate all signals outside the band of interest (by filtering).
Downsample to the minimum required rate, in order to reduce the computational complexity of the baseband processing algorithms.

One common way to accomplish this is the way you describe: downconvert the band of interest to DC, or to a low intermediate frequency if DC offset is a concern, and then decimate.

There could be a few reasons why you're seeing that filtering is slower than decimation. One is that maybe your filter order is larger than necessary. Another is that there are algorithms that perform low-pass filtering and downsampling simultaneously and with very low complexity. See for example half-band filters.


How to filter sudden spikes and baseline changes in microphone signal? The orange is a signal I recorded, the blue - my first restoration attempt, after which I decided to ask.

The signal is audible, however there are 2 types of defects 1) baseline changes 2) sudden drops/spikes (which are not single point):

This doesn't look like smoothing/IR filtering problem, but rather some detection + arithmetic. Maybe I can code solution myself, but are there some ready-made solutions?

Assuming some processing delay is acceptable (100 ms should be more than sufficient), and given there is such a substantial difference between the spike peaks and the signal level, a robust solution would be a non-linear puncture based filter by removing the spikes based on a threshold detector (that then replaces the the required number of samples before and after the threshold with a prior running mean of the filtered signal) and then pass the resulting signal through a low pass filter that has a duration in time at least ten times longer than the interval of time removed due to a spike. Use a running mean and rms of the filtered output for dynamic processing. The mean is described above to replace removed samples as a reasonable minimum error estimate of the sample otherwise hidden and the standard deviation is used to set the threshold. A threshold of 3-sigma above the mean should work well, but ultimately you could also track the statistics of the peak and place the threshold based on optimized probability of false alarm versus probability of detection trades.


Discrete Time Fourier Transform to Frequency Domain I failed to do this question on the exam and finding it very difficult, I would be glad if you can help me solve it. How shall I start?

Most introductory texts on DSP have solved examples of this kind, so you should probably get some good text and work with it. But I'll give you a few hints anyway.

You know that

$$X(e^{j\Omega})=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\Omega}\tag{1}$$

So for (a) you get

$$X(e^{j\Omega})=\sum_{n=-\infty}^{\infty}\left(\frac34\right)^n u[n-4] e^{-jn\Omega}\tag{2}$$

Since $u[n-4]$ equals zero for $n<4$, (2) simplifies to

$$X(e^{j\Omega})=\sum_{n=4}^{\infty}\left(\frac34\right)^n e^{-jn\Omega}=
\sum_{n=0}^{\infty}\left(\frac34\right)^{n+4}e^{-j(n+4)\Omega}=\\=
\left(\frac34\right)^4e^{-j4\Omega}\sum_{n=0}^{\infty}\left(\frac34\right)^{n}e^{-jn\Omega}\tag{3}$$

where you can evaluate the final sum using the formula for a geometric series with $r=\frac34e^{-j\Omega}$.

For (b) you have

$$X(e^{j\Omega})=\sum_{n=-\infty}^{\infty}a^{|n|}e^{-jn\Omega}=\sum_{n=0}^{\infty}a^ne^{-jn\Omega}+\sum_{n=-\infty}^{-1}a^{-n}e^{-jn\Omega}=\\=
\sum_{n=0}^{\infty}a^ne^{-jn\Omega}+\sum_{n=1}^{\infty}a^{n}e^{jn\Omega}=
\sum_{n=0}^{\infty}a^ne^{-jn\Omega}+ae^{j\Omega}\sum_{n=0}^{\infty}a^{n}e^{jn\Omega}$$

Now you have both sums in the desired form to apply the formula:

$$X(e^{j\Omega})=\frac{1}{1-ae^{-j\Omega}}+\frac{ae^{j\Omega}}{1-ae^{j\Omega}}$$

You can combine the two terms, which should finally result in

$$X(e^{j\Omega})=\frac{1-a^2}{1-2a\cos\Omega+a^2}$$

Please remember that this is basic stuff that you must learn to do yourself. The only way to learn it is by solving the problems which you find at the end of each chapter of any introductory DSP text.


Why does a simple delay system with $h[n]=\delta[n-n_0]$ have a linear phase? Let's consider a simple delay system with impulse response $$h[n]=\delta[n-n_0]$$ We also know that it has linear phase with group delay of $n_0$
.

In general, a FIR filter with symmetric (or anti-symmetric) coefficients with filter order $M$ always has linear phase and its (group) delay is $M/2$. Now my question is:
1) What is the filter order of the simple delay system?
    My guess is that the order is $n_0$ (since it is number of unit-delays)
2) Does it have symmetric filter coefficients?
    Apparently it is not, but not quite sure.
3) Why does it have a linear phase?

In response to "1) What is the filter order of the simple delay system?"

The order of a filter is the power of the highest nonzero coefficient of the Z-transform. For the simple delay system given, the Z-transform is given as $H(z) = z^{-n_0}$ so it has order $n_0$.

In response to "2) Does it have symmetric filter coefficients?"

The coefficients are not symmetric if you stop at $n_0$, but adding additional zero coefficients to the signal will not alter it and by doing so, one may describe an equivalent filter that has coefficient symmetry. For example, if the delay is 2 and our original filter is [0, 0, 1], then we can add two more zeros to obtain [0, 0, 1, 0, 0]. These two filters have identical behavior, but the second one can be said to have symmetric coefficients. I would argue that a simple delay does have symmetric coefficients since it can always be put into an equivalent form that has symmetric coefficients.

In response to "3) Why does it have a linear phase?"

The simplest way to understand why a delay results in a linear phase is to examine the output for a sinusoidal input that is delayed by a constant amount. If we have some sinusoidal signal such as $$x[n] = \cos\left(2 \pi f_0 n\right)$$ and we consider what happens after it is passed through our system $$y[n] = x[n] * h[n]$$ which in our case is simply
$$y[n] = x[n-n_0] $$
or if we substitute in $x[n]$ we obtain
$$y[n] = \cos\left(2 \pi f_0 n - 2 \pi f_0 n_0\right).$$
Notice that the resulting output of the system is a sinusoidal signal with the same frequency as the input, but with an additional phase shift. In this case, the phase shift is $2 \pi f_0 n_0$. However, if we are to consider the phase as a function of the input frequency of the sine wave (that is instead of fixing the frequency as $f_0$, we allow it to vary and rename it $f$), then we obtain a phase function
$$ \phi(f) = 2 \pi n_0 f $$
and from this, you can see that the phase response varies linearly in $f$ with a slope of $2 \pi n_0$ radians. We can use this phase function to predict the resulting phase offset at various frequencies. And this prediction is related to the frequency through a simple linear relationship. That is why we say that it has linear phase.