Morphological Gradient Versus Linear Gradient in Image Processing To obtain the gradient from an image, one can use a linear filter that for example is derrived from the derrivative of a gaussian. Then convolving the image with this linear filter gives you the gradient image.
The gradient magnitude can also be approximated by use of (grayscale) morphological gradients, such as: (dilation-erosion), (opening-closing), or ((dilation-erosion)-(opening-closing)).
Why would you resort to using a (grayscale) morphological approximation of the gradient magnitude? Is it any faster or is it different in some cases?
Thanks,
Tom

I've spent few years working with specialists on mathematical morphology, and I didn't get a mathematical answer why morphological gradient was better than a classical convolution. It was more a choice of ideology.

However, the results are a little bit different and according to the problem, you can choose one over the other.

But morphological gradient is super fast, because it is usually performed with a small structuring element (radius 1 or 2) of type square or hexagon, so you can use SIMD/SSE optimization.


Is there an algorithm to compute the phase for a single frequecy? If you have a function $f(t)=A \cdot \sin(\omega t+\phi)$, and reference sin wave $\sin(\omega x)$ what would be a fast algorithm to compute $\phi$?

I was looking at Goertzel algorithm, but it doesn't seem to deal with phase?

The problem can be formulated as (nonlinear) least-squares problem:

$$F(\phi) = \frac{1}{2}\sum^{n}_{i=1}\left[ A \cdot \sin(\omega i+\phi) - f_{i}(\omega) \right]^{2}$$

where $F(\phi)$ is the objective function to minimize with respect to $\phi$.

The derivative is very simple:

$$F'(\phi)=\sum^{n}_{i=1} A\cdot \cos(\omega i+\phi)  \left[A \cdot \sin(\omega i+\phi) - f_{i}(\omega)\right]$$

The above objective function can be minimized iteratively using Gradient descent method (first order approximation), Newton method, Gauss-Newton method or Levenberg-Marquardt method (second order approximation - $F''(\phi)$ need to be provided in these).

Obviously, the above objective function has multiple minima because of periodicity, hence some penalty term can be added to discriminate other minima (for example, adding $\phi^{2}$ to the model equation). But I think the optimization will just converge to the nearest minima and you can update the result subtracting $2\pi k, k\in N  $.


Is this phase response linear phase? Would the following frequency response be considered linear phase? I know that linear phase is when the phase response is a linear function of frequency. However, the phase response is only linear from ~0.3 so I'm not sure if it is considered "linear".

The phase response is clearly non-linear, as mentioned in a comment. Note that phase jumps at the zeros of the frequency response occur depending on how the phase is defined, but these phase jumps don't make the phase non-linear. However, in the given example the phase is clearly non-linear at low frequencies.

You could also look at the corresponding impulse response. If the impulse response has a finite length and if the coefficients are symmetrical or anti-symmetrical with respect to the center, then the phase is linear.

Another method is to look at the pole-zero diagram of the transfer function. If the system is causal then all poles must lie at the origin and all zeros must either lie on the unit circle, or they must occur in quadruples $z_0$, $z_0^*$, $1/z_0$, $1/z_0^*$.


Compute time constant of discrete exponential decay pulses I want to retrieve (within a FPGA) the time constant $\tau$ of incoming pulses of the form $$x[n]=Ae^{-n/\tau}u[n]+C$$
In addition to the offset, inputs are subject to noise.

To measure the integral of the pulses, I used a trapezoidal filter (2 successive differentiator+accumulator) whose properties allow for DC offset removal and noise attenuation. Starting from that, I was wondering if there is were a smart way to compute an approximation of the time constant.

A few articles mention Fourier Transform or Corrected Successive Integration technics but I could not find any free material.

Any tip to put me in the right direction would be appreciated.

Take four sample values.

$$ y_1 = x[p] $$
$$ y_2 = x[q] $$
$$ y_3 = x[p+d] $$
$$ y_4 = x[q+d] $$

Subtract the third from the first.  Notice that the $C$s cancel out.

$$ y_1 - y_3 = x[p] - x[p+d] = Ae^{-p/\tau}-Ae^{-(p+d)/\tau} =  Ae^{-p/\tau} \left( 1 - e^{-d/\tau} \right) $$

Similarly, subtract the fourth from the second.

$$ y_2 - y_4 = Ae^{-q/\tau} \left( 1 - e^{-d/\tau} \right) $$

Take the quotient of these two differences.  Notice that the $A$s cancel out.

$$ \frac{y_1 - y_3}{y_2 - y_4} = e^{(q-p)/\tau} $$

Take the log of each side.

$$ \ln \left( \frac{y_1 - y_3}{y_2 - y_4} \right) = \frac{q-p}{\tau} $$

Solve for $\tau$.

$$ \tau = \frac{q-p}{ \ln \left( \frac{y_1 - y_3}{y_2 - y_4} \right) } $$

I did this more thoroughly in another answer, but this should give you start.  You want the differences to be large compared to the noise level, meaning d should be large.  $(q-p)$ being larger will help too.  I am guessing that four evenly spaced samples are as good as anything.  Something to investigate.

If noise is a problem you can spread the calculation across more points.  You can either do an average of neighbor points for each $y$ value as long as you use the same weighting or employ more points directly in the quotient as in my referenced answer.

The other answer will also show you how to find the best fit values of $A$ and $C$.

Using more points to mitigate noise.

Here is a way you can fully utilize all your samples to reduce the effects of noise I believe as much as possible.

Select the number of samples ($N$) to be a multiple of four.  Partition the samples into four quarters.  Then set the $y$ values to be the sum of the samples of the corresponding quarters.  You don't need to take an average (an extra division) because that washes out in the taking of the quotient.

$$ q-p = N/4 $$

In addition to taking one measure of the entire curve, you might want to try a sliding window taking many smaller measures.  The values you get should be fairly constant.  If not, your assumption of exponential decay may not be correct.

Sample Python code:

import numpy as np

#================================================
def main():

#---- Set Parameters

        A   = 1.234
        C   = 5.678
        tau = 24.08016032

        N = 100

#---- Construct the pulse

        x = np.zeros( N )

        for n in range( N ):
            x[n] = A * np.exp( -n / tau ) + C

        print x

#---- [Add noise here]

#---- Example Using a Single Point

        p = 10
        q = 20
        d = 20

        y1 = x[p]
        y2 = x[q]
        y3 = x[p+d]
        y4 = x[q+d]

        S = ( y1 - y3 ) / ( y2 - y4 )

        #++[Check S here]

        tau_calc1 = ( q - p ) / np.log( S )

        print tau_calc1

#---- Example by Quarters

        y1 = np.sum( x[p:p+10] )
        y2 = np.sum( x[q:q+10] )
        y3 = np.sum( x[p+d:p+d+10] )
        y4 = np.sum( x[q+d:q+d+10] )

        S = ( y1 - y3 ) / ( y2 - y4 )

        #++[Check S here]

        tau_calc2 = ( q - p ) / np.log( S )

        print tau_calc2

#================================================
main()


What is the difference between normalized peak of correlation, versus peak of correlation divided by its average? Given a template, and a signal, the question arises as to how similar the signal is to the template.

Traditionally a simple correlation approach is used, whereby the template and a signal are cross-correlated, and then the entire result normalized by the product of both their norms. This gives a cross-correlation function that can range from -1 to 1, and the degree of similarity is given as the score of the peak therein.

How does this compare to taking the value of that peak, and dividing by the mean or average of the cross-correlation function?
What am I measuring here instead?

Attached is a diagram as my example.

In order to get the best measure of their similarity, I am wondering if I should look at:

Just the peak of the normalized cross-correlation as shown here?
Take peak but divide by the average of the cross-correlation plot?
My templates are going to be periodic square waves with some duty cycle as you can see - so should I not also somehow exploit the other two peaks we see here?

What would give the best measure of similarity in this case?

Thanks!

EDIT For Dilip:

I plotted the cross-correlation squared VS a cross-correlation that isnt squared, and it certainly does 'sharpen' the main peak over the others, but I am confused as to what calculation I should be using to determine similarity...

What I am trying to figure out is:

Can/Should I use the other secondary peaks in my calculations of similarity?
We have a squared cross-correlation plot now, and it certainly sharpens the main peak, but how does this help in determining final similarity?

Thanks again.

EDIT For Dilip:

  The smaller peaks don't really help in similarity calculations; it is the main peak that matters. But the smaller peaks do provide support for the conjecture that the signal is a noisy version of the template. "

Thanks Dilip, I am a little confused by that statement - if the smaller peaks do in fact provide support that the signal is a noisy version of the template, then doesnt that also aid in a measure of similarity?

What am I confused about is whether I should simply use the peak of the normalized cross-correlation function as my one and final measure of similarity and 'not care' about what the  rest of the cross-corr function does/looks like, OR, should I take the peak value and some_other_metric of the cross-cor into account as well.

If only the peak matters, then how/why would squaring the function help, since it just magnifies the main peak relative to the smaller ones? (More noise immunity?)
Long and short: Should I care about the peak of the cross-correlation function only as my final measure of similarity, or should I also take the entire cross-correlation plot into account as well? (Hence my thought about looking at its mean).

Thanks again,

P.S. Time delay in this case is not a problem, in that, it is 'not cared' about for this application.
P.P.S. I do not have control over the template.

To add some perspective, you could go back to the probabilistic interpretation of correlation.
Remember that the correlation is used because it measures the degree of similarity under some linear generative model of the signal knowing the template (additive gaussian noise). The autocorrelation gives a measure of the log-Probability.

To go back to your question, there are a number of free parameters, notably the variance of noise, and this relates to your choice of normalisation.


Why does this Moiré pattern look like this? I was making some gifs of Mobius transformations in Matlab, and some strange patterns began to appear. I'm not sure if a deeper knowledge of the filetype/algorithm is needed to understand this phenomenon, but I thought that there could perhaps be a purely mathematical explanation. The image is obtained by coloring the complex plane like a checkerboard, and then inverting it by taking the reciprocal of the complex conjugate. Here is the math psuedocode for the image with a given zoom $k$:

\begin{align}
&\mbox{checkerboard}:\mathbb C \to\{\mbox{black},\mbox{white}\}\\
&\mbox{checkerboard}(z):=\begin{cases}
\mbox{black} & \mbox{if }\lfloor\Im(z)\rfloor+\lfloor\Re(z)\rfloor\equiv 0\mod 2\\
\mbox{white} & \mbox{if }\lfloor\Im(z)\rfloor+\lfloor\Re(z)\rfloor\equiv 1\mod 2
\end{cases}\\
&\mbox{image} = \{z\in\mathbb C:|\Re(z)|,|\Im(z)|\leq 1\}\\
&\mbox{color}:\mbox{image}\to\{\mbox{black},\mbox{white}\}\\
&\mbox{color}(z):=\mbox{checkerboard}(k/\overline{z})
\end{align}

And here are the pictures for $k=1$, $k=50$, and $k=200$. The resolution of each picture is $1000\times 1000$. I have no background in signal processing, but I would be eager to learn things!

EDIT:

More specifically, why does the Moiré pattern 'sync up' with the resolution of the picture at certain points?
Can the Moiré pattern be predicted?

You'll need to understand the sampling theorem. In short, each signal has what we call a spectrum¹, which is the Fourier transform of the signal as it comes in time domain (if it is a time signal), or spatial domain (if it is a picture. Since the Fourier transform is bijective, a signal and its transform are equivalent; in fact, one can often interpret the Fourier Transform as change of basis. We call that "conversion to frequency domain", since the Fourier transform's values for low ordinates describe the things that change slowly in the original (time or spatial) domain signal, whereas high-frequency content is represented by Fourier transform values with high position.

Generally, such spectra can have a certain support; the support is the minimal interval outside of which the spectrum is 0.

If you now use an observing system whose ability to reproduce frequencies is limited to an interval that is smaller than said support (which often is infinite, by the way, and always is infinite for signals that have finite extension in time or space), you can not represent the original signal with that system.

In this case, your picture has a certain resolution – which is, in the end, the fact that you evaluate the value of your function at discrete points in a fixed, non-infinitesimal spacing. The inverse of that spacing is the (spatial) sampling rate.

Thus, your picture cannot represent the original signal – it's simply mathematically impossible that the mapping of underlying function to pixels is truly equivalent to the original function, since we know that in this case, the total range of frequencies representable by your evaluation at discrete points ("sampling") is half the sampling rate, and thus, something must go wrong with the part of your signal's spectrum that is above half the sampling rate.

What happens is, in fact, that the spectrum gets aliases – every spectral component at a frequency $f_o \ge \frac{f_\text{sample}}{2}$ gets "shifted" down by $n\cdot f_\text{sample},\, n\in \mathbb Z$, so that $|f_o-nf_\text{sample}| <  \frac{f_\text{sample}}{2}$. In effect, that leads to "structure" where there (feels like) shouldn't be some.

Take the "large" structures from your picture that I've painted green:

It certainly looks like there is low-frequency content here - but in reality, it's just the high-frequecy content at frequencies $>\frac{f_\text{sample}}2$ that got aliased to low frequencies, since it was close to an integer multiple of the sampling rate.

So, yes, you can predict the artifacts that happen to a 2D signal when being sampled by comparing its Fourier transform to the bandwidth offered by the sampling rate.

¹ this might be different from the spectrum as used in linear algebra to describe the Eigen-properties of operators.


Machine Vision - Infinity Focus effect I am working on a research project that involves machine vision. The camera we are using implements infinite focus - I have been informed that this can affect the accuracy of computer vision algorithms, why is this?

At the moment we are implementing ORB for feature matching, out of a total of 175 possible detection events 10 of those appear to be missed due to the object not being clearly focused - but it was my understanding that the point of infinity focus was that everything between the lens and "infinity" would be in acceptable focus?

Is there any other effect of infinity focus that could account for this?

"Infinity focus" refers to the Hyperfocal Distance of the camera's lens. This is the distance beyond which all objects will appear acceptably focused.

Therefore, if your objects come closer (to the camera) than its hyperfocal distance they will appear blured and detection of sharp features such as corners will become more challenging.


Geometric Transformation in Gradient Domain Possible? I work with images in gradient domain because I sometimes need to manipulate image gradients (i.e. Laplacian filter responses) rather than intensity values.

However, I need to apply a geometric transform (distortion) to an image already in gradient domain.

Is is then meaningful to reintegrate the distorted image back to spatial domain?

Here is my earlier question related to gradient domain reconstruction with some images.

The geometric transformation needs to be independent from the pixel coordinates $x$ and $y$. assume that $G$ is your transformation matrix:
$$
\int\frac{d^{2} (G I)}{dx dy}  = G \int\frac{d^{2} (I)}{dx dy}
$$
is true only when $G$ is independent from the coordinates.

The restoration will succeed with constant (i.e. global) image rotations, scalings, and shifts.


Reason of effect of continuity and discreteness on periodicity of signal? We know that a sine signal in continuous time is a periodic signal with a period $2\pi$ whereas the same sine signal in discrete time is aperiodic.
My question is how changing just a type of signal (by sampling) can affect the signal's periodicity and why this happens whereas we have just omitted some real values between samples and how can we make use of this fact practically?

To claim that "sampling changes the continuity and periodicity of a signal", implies that given a signal $x(t)$ and its sampled version $q(n) = x(n \cdot T), n \in \mathbb{N}, T = \frac{1}{Fs}$ there is a way to "transform" $x(t)$ in an aperiodic $q(n)$ just by periodic sampling.

This is impossible, especially for linear systems that cannot go into oscillations out of their own accord.

The only thing that you do with sampling is "tuning in" to a part of a signal's spectrum. Here is why, given a periodic signal at some frequency:

$$x(t) = sin(2 \pi f t)$$

We are going to sample it with a periodic train of pulses:

$$q(t) = \text{III}_T(t) \cdot x(t) \implies \frac{1}{T} \sum_{n=-\infty}^{n=\infty}e^{i 2 \pi n \frac{t}{T}} \cdot x(t) \implies $$

$$\frac{1}{T} \sum_{n=-\infty}^{n=\infty}\left( \cos\left(2 \pi n \frac{t}{T}\right) +i \cdot \sin\left(2 \pi n \frac{t}{T}\right) \right) \cdot x(t)$$

This expression by itself, does not guarantee that $q(t)$ is not going to be aliased and this is why $x(t)$ is first going through a low pass anti-aliasing filter. After this operation, sampling, is putting a "canvas" over that antialiased spectrum (discretises it).

How it is then possible that the operation described above can generate an aperiodic $q(t)$ for a periodic $x(t)$?

The only thing that this operation does is sample from the already existing components of $x(t)$ which are periodic functions. If $x(t)$ contained discontinuities, so would $q(t)$. If $x(t)$ was aperiodic, so would $q(t)$.
At the end of the day, what you do here is a straightforward mapping. You don't change any "type of signal".

The case that robert-bristow-johnson mentions is a result from number theory, transposed to periods of sinusoids. You cannot compose an irrational number out of integers and this tiny little difference that remains after each period presents itself as a "running phase" of the sinusoid.

There are however cases where sampling can have an effect on the signal, BUT that is entirely resonable. In the integration of non-linear differential equations by naive schemes such as Euler's method, the "sampling" that is effected by choosing a fixed step size can lead to large errors in the reconstructed waveform and more often than not in cases like this, the sampling becomes irregular. More dense around areas of large errors and vice versa. But even in that case, periodic regions (in the input) would remain periodic in the output.

Hope this helps.


Compare between JPEG and JPEG2000 JPEG image compression is Fourier based  DCT while modern image compression technique like JPEG2000 is based on more multi-scale technique like Wavelets.I want to know how Fourier  and Wavelets are useful in Image compression.

So can anybody  explain advantages and disadvantages of JPEG and JPEG2000 with the help of (characteristics of)  transforms they use?

JPEG is far simpler.  It divides the image into 8x8 pixel blocks, and processes each using a Discrete Cosine Transform.  The results are quantised and then encoded.  The quality is fixed by the encoder.

JPEG2000 uses a 2D wavelet function, the output of which is four "images", each a quarter the size of the original.  One of those is actually an image, while the others are high-frequency components that can be added to it to re-construct the full-resolution image.  The wavelet process may be repeated multiple times.  The result is a tiny image, and a series of high frequency components that may be combined with it.  Each resulting component is quantised and encoded.

JPEG is fine for high quality and modest compression, which is why it is still very widely used.

JPEG2000 offers several advantages:

To achieve very high compression, it is possible to throw away or heavily quantise the high frequency components.  This gives a poor quality - but usable - image where JPEG would fail completely.
Images can be re-constructed progressively at ever improving quality.  This can be either in terms of increasing resolution or bit depth as required.
It supports the JPIP protocol for progressively transmitting images to a client.  The client may retrieve low-resolution thumbnails and then just the parts of the image they want at better resolutions.


Matched Filters for attenuation detection? I have been reading about matched filters and I'm interested in possible uses for them,

I have seen they are used very often in radar, in such a way that you send out a signal, and then you read the incoming signals to look for the reflection of the pulse you sent out using matched filters.

I think I'm right in saying that by the time difference from when you sent out the signal to when you found it again, you can detect how far the object is from you.

What I am interested to know is, can you detect the attenuation of the signal from matched filters and not just the time delay?

Hence giving you some information on the material properties of the reflecting object

Your interpretation of the matched filter for radar seems almost right, just to be sure, let me reiterate a bit: Think of it this way: we're sending a pulse function $p(t)$. If a target in distance $d$ reflects the signal, we receive $\gamma \cdot p(t-\tau)$, where $\tau = 2d/c$ is the total travelled delay.

The problem is: how to sample such a signal? Let's say your pulse was rectangular, you would need a lot more samples than the pulse width to know its extent and therefore, the proper delay. This is where the matched filter helps.

Applying matched filtering to your received signal, you filter it with $p(-t)$. This gives a signal comprised of the autocorrelation function $\varphi_p(t) = p(t) * p(-t)$, i.e., we obtain $\gamma \cdot \varphi_p(t-\tau)$, more or less. Now, the ACF always has its maximum at $t=0$, which greatly helps for the sampling. Also sampling the ACF at zero gives the optimal SNR: you gather the whole energy of the pulse into the sampling point. It's also called pulse compression for this reason.

Now, back to your question: You can extract both from the received signal, the delay $\tau$ and the reflection coefficient $\gamma$. In fact, this is often done. People call it RCS (radar cross section). Note that it is not so easy to infer material properties from $\gamma$ as $\gamma$ is typically highly direction dependent: the geometry of the object has a lot of effect on how strong the reflected wave is.


Squared magnitude of System Function H(z) If:

$$
|\alpha|^2 = \alpha \alpha^*
$$

Then, why does:

$$
|H(z)|^2 = H(z) H(z^{-1})
$$

instead of:

$$
|H(z)|^2 = H(z) H^*(z)
$$

It's important to understand that the equation
$$|H(z)|^2=H(z)H\left(\frac{1}{z}\right)\tag{1}$$
is only valid on the unit circle, i.e., for $z=e^{j\omega}$. For complex-valued systems the general form of $(1)$ is
$$|H(z)|^2=H(z)H^*\left(\frac{1}{z^*}\right)\tag{2}$$
For real-valued systems we have $H(z)=H^*(z^*)$ and, consequently, $(1)$ is valid for real-valued systems. Note that both $(1)$ and $(2)$ are only valid on the unit circle.
The reason why we use $(1)$ or $(2)$ as a generalization of $|H(e^{j\omega})|^2$ to the whole complex plane is that we want an analytic function that is actually a $\mathcal{Z}$-transform of some sequence. And the function $H(z)H^*(z)$ does not satisfy that requirement.
Note that the IDFT of $H^*(e^{j\omega})$ is $h^*[-n]$. So if we find the $\mathcal{Z}$-transform of $h^*[-n]$ we have what we need:
$$\sum_nh^*[-n]z^{-n}=\left(\sum_nh[n](z^*)^{n}\right)^*=H^*\left(\frac{1}{z^*}\right)\tag{3}$$


Why is Scipy implementation of Hilbert() function different from Matlab implementation of the function? I am trying to fit Hilbert envelop to a high frequency ultrasonic signal of frequency 250 KHZ and sampling rate 12000000. Raw signal looks like below: .

I used hilbert() function from scipy.signal package in python this is what it looks like .

The python code looks like below

from scipy.signal import hilbert
import numpy as np
 def Hilbert(self,i=0):
    analytical_signal = hilbert(self.sensor["s"+str(i)])
    amplitude_envelope = np.abs(analytical_signal)
    return amplitude_envelope

The Matlab implementation looks like this

Matlab code is as follows:

figure;
plot(abs(hilbert(signal)),'r');
hold on;
plot(signal,'b');

The line data is as follows Signal.
I am wondering which one is correct?

It works fine for me:

from scipy.signal import hilbert
import numpy as np
from matplotlib.pyplot import plot

sensor = np.loadtxt('signal.txt')
plot(sensor)

analytical_signal = hilbert(sensor)
plot(analytical_signal.real)
plot(analytical_signal.imag)

amplitude_envelope = np.abs(analytical_signal)
plot(amplitude_envelope)

What are you doing differently?  Maybe you're throwing away the imaginary part of analytical_signal somehow?


Downsample Equiripple FIR filter I am designing some Downsampling filters using the remez exchange algorithm (Matlab FDAtool). I see that some times i get a spike at the first and last coefficients. does anybody knows why is that? my intuition seems to be that this can introduce some artifacts to the signal..
What will be the consequance of multiplying the filter by a window (hamming and so on) to reduce this edge effects?

The "spikes" (they are quite small in magnitude, but large compared to the nearby filter taps) at the beginning and end are part of how the equiripple property is achieved.  In his book "Multirate Signal Processing for Communication Systems" (a book that I recommend quite highly), fredric harris indicates that it is sometimes advantageous to eliminate the spikes because it causes the frequency response to continue to decrease after the cut-off frequency instead of staying constant.  The price you pay is that the filter is no longer equiripple which appears to cause (see example below and comments by JasonR) the transition band to be not as sharp.

EDIT:
I created a low-pass filter in Octave using remez.  I intentionally made the transition-band very small to get it to create a decent spike.  The filter taps are plotted below.

Here is the frequency response of the filter.

As you can see, the frequency response is pretty poor due to the short transition-band, but the point is that it is an equiripple filter.  Below I plot the exact same filter, only I have zeroed out the two spike taps.

Below is the frequency response of the "no spikes" version of the filter.

The frequency response shows the continued drop beyond the cutoff frequency.  It also clearly shows the non-equiripple response in the pass-band.  This is entirely expected.  The one thing that is unexpected is that the "no spike" filter appears to be better than the "spike" filter at every point.  I assume that this is because of some flaw in Octave's remez algorithm.


How to convert wave from real to complex and vice versa? I have wave expressed by array of real numbers (double in C++). But I want to express it as a complex.
I tried to create complex variable and assign to its real the array of my wave and to its imaginary just array of zeros. OK it's complex now. But I feel there is something wrong. It's not authentic complex. It's just surrogate. For example I can't manipulate in phase (only reverse phase). When I multiply it by for example $ i^{0.4} $ it's just changing amplitude little bit but phase is without change.

In the other hand if I have authentic complex wave how to express it as a real numbers but including information from imaginary values?

For any help great thanks in advance.

Rather than saying the questions you are asking are non-sensical, I will say they reveal a lack of understanding of complex numbers.  There is lots of material to be had with a few simple searches.  I would also recommend that you read my blog article The Exponential Nature of the Complex Unit Circle which gives an explanation of Euler's equation which is absolutely essential to understand in order to comprehend the meaning of DFT bin values.

For your first question:  Technically, a real number is a complex number with a zero imaginary component.  So setting the real part equal to your values and the imaginary part to zero is the correct method.  When you multiplied it by $ i^{0.4} $ it should have rotated each value a tenth of a cycle.  Thus the amplitude stays the same, but the phase shifts by $\pi/5$ radians.

For your second question:  A complex number consists of two values.  A real number consists of one.  Therefore converting a complex number to a real number is a non-sensical proposition.  However, you can calculate the magnitude, which is a real number.  You can also ask what the real part is, which is also a real number.  The imaginary part is a real number multiplied by $i$.

A real tone is a sinusoidal, flat as a pancake.  A complex tone is a spiral, like a slinky.  The definitions I use in my blog articles are:

$$ S_n = A \cos( \alpha n + \phi ) $$

For a real tone, and:

$$ S_n = A e^{i (\alpha n + \phi)} $$

For a complex one.

Notice that a real valued tone is just the average of two complex tones:

$$ \frac{ A e^{i (\alpha n + \phi)} + A e^{-i (\alpha n + \phi)} }{2} = A \frac{ e^{i (\alpha n + \phi)} + e^{-i (\alpha n + \phi)} }{2} = A \cos( \alpha n + \phi )$$

Because

$$ \cos( \theta ) = \frac{ e^{i\theta} + e^{-i\theta} }{2} $$

Hope this helps.

Ced


Tricks for plotting the magnitude of a DTFT? To preface, this isn't a homework question but rather a self-study question to help me to understand the basics of finding the DTFT and magnitude of the DTFT based on a discrete time signal sampled from a continuous time signal.

I was wondering if anyone has any methodology to approaching questions where you're asked to plot the magnitude of the DTFT of a discrete-time signal, when
$x(t)$ is sampled according to $x[n] = x(nT_s)$`, with $F_s$ equal to a given value.

For example if $$1/T_{s}=8W$$ and :
$$x(t) = T_{s}\frac{1}{4W}\frac{\sin(2\pi Wt)}{\pi t}\frac{\sin(2\pi 2Wt)}{\pi t}$$

What are the first questions that you would ask yourself that will help you to plot the magnitude of the DTFT.  I am having trouble translating the DTFT to the plots themselves.

So the example that you've listed appears to be somewhat odd. Assuming Ts is the time step, you're saying that the continuous signal x(t) is dependent on the time step you've chosen. That just seems odd.

It also seems that you've formatted your question in the sense that you understand the Fourier Transform (FT), but you're just having a difficult time relating the FT to the associated Discrete Time Fourier Transform (DTFT) after sampling the continuous signal. I'm going to speak at a relatively high level to explain how to go from the FT to the DTFT.

The first thing you should realize is that the DTFT will be related to the FT of the continuous signal. But how so? The sampling theorem (which is the most fundamental concept in DSP) tells us that the bandwidth of a discrete signal (i.e. sampled signal) is equal to that of the sample rate$^1$. We typically write this as -Fs/2 to Fs/2 where Fs is your sample rate. Furthermore, we know that analog frequencies that are integer multiples of the sample rate will be interpreted as equivalent frequencies. This is aliasing.

For example, say Fs = 1kHz. Let's randomly select a frequency of 0.5kHz. All frequencies that are integer multiples of the sample rate at this frequency will be seen as 'the same thing' when sampled. 0.5kHz, 1.5kHz, 2.5kHz, etc. etc. will all appear to be 'the same signal' when sampled. This is also true of other frequencies such as 0.25kHz. In this case, 0.25kHz, 1.25kHz, 2.25kHz, etc. will all appear to be the same. Note that this phenomenon happens at EVERY frequency.

This is what we refer to as aliasing: different frequencies be interpreted as the same signal due to sampling. A common misperception is that aliasing only happens in one direction: that being that high frequencies are aliased onto lower frequencies. However this isn't true. Looking at our above examples, you can see that aliasing only comments on the fact that different signals are interpreted the same. It doesn't say anything about direction. Thus you could also state that lower frequency signals are aliased onto high frequency signals. Because of this aliasing of high frequency signals onto low frequency signals, and low frequency signals onto high frequency signals, the frequency spectrum of sampled signal becomes periodic. So let's put this together into a simple list of steps on how to get the DTFT from the FT.

1) Obtain a copy of the FT of the continuous signal. You will need to do this analytically.

2) Slice up the FT along the frequency spectrum at intervals of Fs. Let's say your first interval would range from -Fs/2 to Fs/2 (this is centered at zero).

3) After you've 'cut up' your FT into intervals, you can take all of those and add them together at a single interval (i.e. do the super position of them). By doing this step, we are essentially introducing aliasing into one interval.

4) Now repeat step 3 for every interval. By doing this, you're accounting for the fact that aliasing doesn't just occur from high frequency to low frequency, but it essentially occurs everywhere!

5) Look at what you've done. By doing the super position at every interval, you'll see that your frequency spectrum is now periodic (and still continuous!).

This is essentially what the DTFT is doing. Visualizing these steps is very useful, I might add illustrations if I get time. Now assuming that you've got a plot of the FT, you now have a plot of the DTFT.

[1] It would be more accurate to say that the relevant bandwidth is that of the sample rate.


How to deal with low fundamental when using AMDF for pitch extraction? I am using the Average Magnitude Difference Function to estimate the fundamental frequency of a quasi-periodic audio signal. The AMDF is defined as

$$ D_n = \frac{1}{N-n}\sum_{k=n}^{N-1}|S_k - S_{k-n}| $$

where $N$ is the length of the signal. This function exhibits a minimum when the signal is shifted by an amount equal to its period.

This is the code I am using to extract the pitch (in Matlab):

 a = amdf(f);
 a = a/max(a);
 [p l] = findpeaks(-a, 'minpeakprominence', 0.6);
 pitch = round(sample_freq/l(1);

However, I am dealing with an audio signal where the fundamental frequency is very low:

As a consequence, a pitch doubling problem arises: the detected minimum corresponds to half the period of the signal (i.e. the second harmonic):

I tried to extract the largest peak and not just the first, but sometimes this problem remains. How can I improve my code and/or the AMDF function in order to deal with low fundamental?

This is what we call in the pitch-detection biz, the "octave problem".
First of all, I would change the AMDF to ASDF.  And I would not reduce the window size as the lag increases.  (Also, I am changing notation to what I consider to be more conventional.  "$x[n]$" is a discrete-time signal.)
The Average Squared Difference Function (ASDF) of $x[n]$ in the neighborhood of sample $x[n_0]$ is:
$$ Q_x[k, n_0] \triangleq \frac{1}{N} \sum\limits_{n=0}^{N-1} \left(x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor] \ - \ x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor + k] \right)^2 $$
$\left\lfloor \cdot \right\rfloor$ is the floor() function and, if $k$ is even then $ \left\lfloor \frac{k}{2}\right\rfloor = \left\lfloor \frac{k+1}{2}\right\rfloor = \frac{k}{2} $.
Now, expand the square and consider what the summations look like as $N \to \infty$ (not that $N$ is going to infinity, but to give you an idea if $N$ is large).  The ASDF is directly related to the autocorrelation.  It is essentially the autocorrelation turned upside down.  These steps I will leave to you. take a look at this answer.
So now consider this finite-length "autocorrelation" (in the neighborhood of sample $x[n_0]$) defined from the ASDF:
$$ R_x[k,n_0] = R_x[0,n_0] - \tfrac12 Q_x[k, n_0]  $$
where
$$ R_x[0, n_0] \triangleq \frac{1}{N} \sum\limits_{n=0}^{N-1} \Big(x[n+n_0-\left\lfloor \tfrac{N}{2}\right\rfloor]\Big)^2 $$
This value $R_x[0,n_0]$ is a measure of the mean power of the signal $x[n]$ in the neighborhood of $n \approx n_0$. Since $Q_x[0,n_0]=0$ and $Q_x[k,n_0] \ge 0$ for all lags $k$, that means that $ R_x[k,n_0] \le R_x[0,n_0] $ for all lags $k$.
Another useful way to look at this autocorrelation taking place in the neighborhood centered at sample $x[n_0]$ is to normalize $R_x[k, n_0]$ with $R_x[0, n_0]$:
$$ r_x[k,n_0] \triangleq \frac{R_x[k,n_0]}{R_x[0,n_0]} = 1 - \frac{Q_x[k,n_0]}{2 R_x[0,n_0]}  $$
This normalized autocorrelation has $r_x[0,n_0]=1$ and $r_x[k,n_0] \le 1$ for all other $k$.
Suppose for a minute that $x[n]$ is periodic with period $P$ (and $P$ happens to be an integer), then
$$ x[n+P] = x[n] \quad \forall n $$
and  $Q_x[mP, n_0] = 0$ and $R_x[mP, n_0] = R_x[0, n_0] \ge R_x[k, n_0]$ for any integer number of periods ($m$ is an integer).  So you get a peak at $k=0$ and at $k$ equal to any other multiple of $P$ if $x[n]$ is periodic.  If $x[n]$ is not perfectly periodic, what we might expect is the biggest peak at $k=0$, another peak (but slightly smaller) at $k=P$ (the period we are looking for) and progressively smaller peaks for larger multiples of $P$.
We can then expect that the value of the normalized autocorrelation, $r_x[k,n_0]$ evaluated at a lag of $k=P$ or other multiples of $P$ should be pretty close to 1.  That value $r_x[P,n_0]$ can be thought of as a measure of the degree of periodicity (sometimes called the pitch confidence) of the estimated period $P$ for the quasiperiodic $x[n]$ in the neighborhood of $n \approx n_0$.  If $r_x[P,n_0]=1$, we can say that $x[n]$ is perfectly periodic with period $P$.  If the best $r_x[k,n_0]$ you can get (with $k$ that's not close to $k=0$) is very small, then $x[n]$ shows no periodicity and your pitch confidence is low.
So the octave problem comes about because of a couple of reasons.  First of all, $P$ is not necessarily an integer.  That is an interpolation problem, not a big deal.
The second reason and more difficult problem is that of subharmonics.  Consider that you're listening to a nice periodic tone at exactly A-440 Hz and it sounds like an A that is 9 semitones above middle C.  Now suppose someone adds to that tone a very tiny-amplitude (like down 60 dB) A-220?  What will it sound like and mathematically what is the "true" period?

Choosing the "right" peak for the period.
Let's say you run your note through a DC-blocking filter, so that the mean of $x[n]$ is zero.  It turns out that causes the mean of the autocorrelation $R_x[k, n_0]$ for every $n_0$ to also be zero (or close to it if $N$ is large).  That means $R_x[k, n_0]$ must sum (over $k$) to be about zero which means there is as much area above zero as below.
Okay, so $R_x[0, n_0]$ represents the power of $x[n]$ in the vicinity around $n=n_0$ and must be non-negative.  $R_x[k, n_0]$ never exceeds $R_x[0, n_0]$ but can get as large as it when $x[n]$ is periodic.  $R_x[P, n_0] = R_x[0, n_0]$ if $x[n+P]=x[n]$.  So if $x[n]$ is periodic with period $P$ and you have a bunch of peaks spaced apart by $P$ and you have an idea for how high those peaks should be.  And if the DC component of $R_x[k, n_0]$ is zero, that means in-between the peaks, it must have negative values.
If $x[n]$ was "quasi-periodic", one cycle of $x[n]$ will look a lot like an adjacent cycle, but not so much like a cycle of $x[n]$ farther down the signal in time.  That means the first peak $R_x[P, n_0]$ will be higher than the second at $R_x[2P, n_0]$ or the third $R_x[3P, n_0]$.  One could use the rule to always pick the highest peak and expect the highest peak to always be the first one.  But, because of inaudible subharmonics, sometimes that is not the case.  sometimes the second or possibly the third peak is oh-so-slightly higher.  Also, because the period $P$ is likely not an integer number of samples but $k$ in $R_x[k, n_0]$ is always an integer, so the true peak will likely be in-between integer values of $k$.  Even if you were to interpolate where the smooth peak is (which I recommend and quadratic interpolation is good enough), and how high it really is between integer $k$, your interpolation alg could make a peak slightly higher or slightly lower than it really is.  So choosing the absolutely highest peak can result in spuriously picking the second over the first peak (or vise versa) when you really wanted the other.
So somehow you have to handicap the peaks at increasing $k$ so that the first peak has a slight advantage over the second, and the second over the fourth (the next octave down), etc.  How do you do that?
You do that by multiplying $R_x[k, n_0]$ with a decreasing function of $k$ so that the peak at $k=2P$ is reduced by some factor, relative to an identical peak at $k=P$.  Turns out that the power function (not the exponential) does that.  so compute
$$ k^{-\alpha} \ R_x[k, n_0] $$
So, if $x[n]$ were perfectly periodic with period $P$, and ignoring interpolation issues for non-integer $P$, then
$$ R_x[2P, n_0] = R_x[P, n_0] $$
but
$$\begin{align}
 (2P)^{-\alpha} R_x[2P, n_0] & = \\
 (2P)^{-\alpha} R_x[P, n_0] & < P^{-\alpha} R_x[P, n_0] \\
 \end{align}$$
The factor by which the peak for a pitch of one octave lower is reduced is the ratio
$$ \frac{(2P)^{-\alpha} R_x[2P, n_0]}{P^{-\alpha} R_x[P, n_0]} = \frac{(2P)^{-\alpha}}{P^{-\alpha}} = 2^{-\alpha} $$
So if you want to give your first peak a 1% boost over the second peak, which means you will not choose the pitch to be the sub-harmonic pitch, unless the sub-harmonic pitch autocorrelation is at least 1% more than the first peak, you would solve for $\alpha$ from
$$ 2^{-\alpha} = 0.99 $$
That is the consistent way to weight or de-emphasize or handicap the peak corresponding to the subharmonic pitch one octave below.
It still leaves you with a thresholding issue.  You have to choose $\alpha$ well. But this is a consistent way emphasize the first peak over the second, which is an octave lower, but not so much that if the note really is an octave lower, but the energy in all of the even harmonics was strong, compared to the odd harmonics, this will still leave a possibility for the second peak being chosen.


Difference between librosa.stft and scipy.signal.stft I read the source code of librosa.stft and sicpy.signal.stft, and notice that the calculation results of STFT(short-time fourier transform) in these two libraries are quite different:
In scipy.signal.stft, the stft result is scaled by 1.0/win.sum(), while in librosa.stft no scaling or normalization procedure is done. So why scipy.signal.stft do the additional scaling procedure? And is there any other difference in the calculation of STFT in these two libraries?

Fft scaling (normalization) isn't strictly a part of an fft algorithm although some implementations (like the one you mention) include it.
Since scaling depends on the actual fft algorithm, it makes sense, in my opinion,  to include it in the implementation because it may not be obvious for programmers to decide what scaling factor(s) to use (the usual scaling factors are $1/N$, $2/N$ (an efficient fft algorithm that is using only half the spectrum so we need to multiply by 2 to compensate for the missing energy), $1/sum(wnd)$, $2/sum(wnd)$
As for your second question, I haven't gone over the code (you haven't provided any links to the source), but some (usually minor) calculation differences might exist
given that an fft algorithm can be implemented in a number of ways ( e.g real fft vs complex fft, decimation in frequency vs decimation in time,
radix 2 vs mixed radix etc and combinations thereof)


sequence to perform filtering? which filter first? Hi all Is there a generally accepted rule as to which order to perform filtering in?

for example if I have three iir filters, a low pass filter a notch and a high pass filter, which sequence would it be best to perform these in? or does it not make any difference at all?

cheers, Graham.

Theoretically (with infinite precision arithmetic) order of LTI filters does not matter. If any of the filters is not LTI then order will be important.

Plus, for finite precision arithmetic one should also take care of numerical error build up and potential instability issues especially with IIR filters where pole locations may become critical


Matching delayed signals I have pairs of 1d digitised waveform signals which are almost identical - except there are sections which are the same in both but one is delayed slightly.
I need to find sections and the delay with the corresponding section in the other channel.

I can set a maximum delay BUT I need to measure a delay to less than the sample interval - so some sort of convolution/correlation approach rather than a simple feature id.

Any suggestions where to start?

This may be a good candidate for a two step approach:

Step 1 would be a pretty coarse running cross correlation with a threshold detector to identify what parts in the signals are matching.

Step 2 would than determined the actual delay. There are various ways to get sub sample resolution:

Upsample to desired resolution and then cross correlate
Short term fourier transform and match a linear phase difference with a weighted least squares approach
Delayed lock loop with a fractional delay filter
Match the the two signals in with an adaptive filter. Then calculate the fourier transform of the filter impulse response and calculate the delay through a linear phase match


why did convolution coding achieve BER less than Reed-Solomon? I tried convolutional coding and Reed-Solomon (RS) (block coding technique) for wireless channel where the input data is bitstream.

Convolutional coding achieved lower BER comparing with RS. If the SNR of channel decrease around 15dB, convolution code's BER becomes more than that of the system without channel coding, where RS achieved a small improvement in the BER.

I wonder if any one has an explaination for what happened ?

Reed-Solomon codes are good for correcting bursts of errors (such as those that often occur in computer memory systems such as RAM, or on less volatile media such as CDs and DVDs) or at the output of block decoders -- either a received word is completely correct, or it has multiple errors all within the codeword length -- or at the output of convolutional decoders -- if a decoder goes off the correct path, there is a burst of errors until the decoder gets back on the correct path. Reed-Solomon codes are not all that good for correcting random errors such as those that occur in AWGN channels.  It has been common practice for many students (as well as their faculty member advisors, who really ought to know better), not to mention engineers in the communications industry, to put up graphs comparing how well convolutional codes do compared to Reed-Solomon codes in preparation for touting the advantage of their thesis work or product over the competition.
As a matter of fact, systems needing to meet the most stringent performance criteria often use a concatenated coding system in which the outer code is a Reed-Solomon code and the inner code (the one whose codewords are transmitted over the channel) is a convolutional or turbo or LDPC code, all of which are well-matched to the characteristics of the channel.


Intuition for sidelobes in FFT I was wondering if there's a intuitive way to understand why sidelobes appear when performing an FFT on a signal of fixed length?

I have two explanations with the intent of providing additional intuitive insight beyond the concise mathematical explanation; first, an explanation from the Fourier Series Expansion along with the idea of repeating the truncated time domain sequence which shows the implied discontinuities that result, requiring more frequency components to reconstruct than would really exist if the waveform was not truncated. And second, an explanation from looking at the DFT as a bank of not-so-great filters.

First Explanation: Fourier Series Expansion  and Time Periodicty

The Fourier Transform of a time limited sequence is identical to the Fourier Transform of a sequence that is periodic for all time.

This is true for the Fourier Transform as well as the Discrete Fourier Transform and is known as the Periodicity Property of the DFT:

Periodicity Property: Given the N-point DFT vector X[k], with an Inverse DFT samples x[n], with both k and n ranging from 0 to N-1; if n lies outside of the range of 0,1,2...,N-1, then

$x[n] = x[modulo(n,N)]$

Likewise, given a N-point time sequence x[n], with a DFT X[k], if k lies outside of the range of 0, 1,2,..., N-1, then

$X[k] = X[modulo(k,N)]$

To help gain an intuitive understanding, a key take-way is that anything that is sampled in one domain, becomes periodic in the other domain. Similarly, anything that is periodic in one domain, becomes sampled (discrete-valued) in the other domain. Here specifically, "sampled" means that the signal will only exist as non-zero values at discrete locations in the domain (stream of impulses).

Sampling in One Domain -> Periodicity in the Other Domain: Consider the spectrum of a 3 Hz analog cosine wave before and after sampling with an A/D converter. The digital spectrum can be viewed as periodic; a cylinder view of the digital spectrum is also a valid view to explain the periodicity, but I find this extension to the analog frequency domain (by having the frequency extend to +/- infinity) helps some gain an intuitive view of the signal processing involved.

Periodicity in one Domain -> Sampling in the other Domain: Consider the Fourier Series Expansion as a simple example demonstrating this property.  The Fourier Series Expansion is done over a finite time domain interval from 0 to T. When decomposed into separate frequency components, the only frequencies used are DC, the fundamental frequency 1/T, and integer multiples of 1/T (harmonics). In effect, since the frequencies can only exist at multiples of 1/T (and DC),  the frequency domain has been sampled.

Further, if we reconstruct the time domain waveform, by summing together the individual frequency components, we can also view the implied periodicity in the time domain if we allow the frequency components to extend beyond the interval 0 to T. It is because of this periodicity that frequency components cannot exist at any other frequencies than multiples of 1/T (because of the opposite condition: if they did exist, then they would not begin and end consistently in the 0 to T time interval, and therefore periodicity cannot exist).

Understanding the above will hopefully help provide an intuitive explanation of spectral leakage.  So now I will repeat a main point:

The Fourier Transform of a time limited sequence is identical to the Fourier Transform of a sequence that is periodic for all time.

Spectral Leakage with the "Fourier Series Expansion" View:

Consider two sinusoidal waveforms in the time interval from 0 to T, the first with an integer number of cycles over the time interval, and the second case with a non-integer number of cycles.

Clearly in Case 1 we can see with the repetition view that even repeating our pure sinusoid remains a pure sinusoid, but in Case 2 our sinusoid now suffers from abrupt transitions, and using the Fourier Series Expansion view of reconstruction, would requred several frequency components to reconstruct such a waveform in the time domain.

Second Explanation: Filter Bank View of the DFT

Another intuitive explanation for spectral leakage (and helps significantly with understanding the DFT in general) is what I call the filter bank view of the DFT.
To see this, consider a simple 4 pt DFT as shown in the figure below and observe for each bin, we are effectively rotating the signal, and then passing the rotated values through a 4 tap unity gain FIR filter. For the first bin, which corresponds to DC, there is no rotation, so we just sum the four samples, and for the other bins we progressively rotate at higher frequencies as we move through the DFT bins:

(Side note- If we performed a streaming DFT, where we calculated a new 4 point DFT on a 4 point sequence as we scanned through a waveform, it would exactly be such a filter bank, but regardless if we do that or not, this view gives great insight into spectral leakage in addition to the usual convolving of a sinc function in the frequency domain that the mathematical explanation reveals)

Now consider the frequency response for each equivalent FIR filter, using the coefficients given in the DFT (For example use freqz([coeff]) in Matlab or Python) as shown in the figure below. (Note that the filter coefficients are in reverse order from that shown in the DFT matrix. This is because the coefficients of the filter are reversed in time; the earliest samples in time are at the far end of the FIR filter state register, so reversed in time):

Here is the main point: Since each Filter in the construction of the DFT is fundamentally a unity gain FIR filter, the shape of this filter in frequency approaches a sinc function as the length of the DFT gets longer (and is an aliased sinc function for small N). So we will call these sinc filters and note that a sinc filter has relatively high sidelobes and has an envelope that rolls off very slowly with frequency (at 1/f for a pure sinc). With the phase rotators in the DFT, we just move the main lobe of this sinc filter to each bin of interest, but the side lobes that exist for each bin allow frequencies at other locations to cause energy to appear in this bin. The amount of leakage is completely predicted by these filters.

Using this view, consider an input single-toned signal with a frequency that is somewhere between two frequency bins as shown in the figure below (an input that exists exactly on any bin will have an integer number of cycles in the time domain interval and therefore NO spectral leakage as we showed earlier). The top filter shows the amplitude at this frequency location that will "leak" to the first bin. The second filter shows the amplitude (slightly higher), the third bin (that our frequency in closest to) will have the highest response and the fourth bin will be lower.

SUMMARY

I have presented two explanations with the intent of providing additional intuitive insight beyond the concise mathematical explanation of multiplication by a rectangular window in the time domain is convolution in the frequency domain (and the leakage we see therefore is the result of a sinc function convolving in frequency with our waveform of interest that we truncated in time); first, an explanation from the Fourier Series Expansion along with the idea of repeating the truncated time domain sequence which shows the implied discontinuities that result, requiring more frequency components to reconstruct than would really exist if the waveform was not truncated. And second, an explanation from looking at the DFT as a bank of filters, and poor filters at that (specifically unity gain filters that approach a sinc function frequency response as N increases). Such filters provide the best frequency resolution but poorest dynamic range.


FFT requirements of real time signal from sensors I have an electrical current sensor which generates measurements every second. I want to calculate the FFT of this signal. Since my real-time values are received every second, I think that the fastest I can go with the FFT sampling period is also 1 sec. In that case, what I am basically doing is that I collect X number of samples in a buffer and I calculate the FFT. I then store the new values as they arrive at the end of this buffer and repeat the FFT calculation.
Q1: is this a correct setup?
Q2: is there a minimum requirement for the buffer (FFT window?). I have now set it up to process 10 samples.
Q3. what is the maximum number of frequencies that can be detected in this scenario?

Q1: is this a correct setup?

I'm guessing "no" but that really depends on your application and goal. The process will produce numbers but if these numbers are useful or not depends on what you want to do with them.

Q2: is there a minimum requirement for the buffer (FFT window?).

Mathematically speaking : no. You can run a FFT on a single sample if you want to. t doesn't do anything useful but it can be done. The real requirements come from your application.

I have now set it up to process 10 samples. Q3. what is the maximum number of frequencies that can be detected in this scenario?

An FFT setup like this will you give  you ten numbers which are the complex amplitudes from -0.5Hz to +0.5Hz in 0.1 Hz steps. That's not the only frequencies you "detect" but that's a more complicated topic. If you feed a 0.15Hz sine wave in the process all 10 numbers will be non-zero.
You may also get significant aliasing with this setup, if you don't use a decent anti-aliasing filter before the current sensor.


time in Sine Wave equation According to Wikipedia,

  The sine wave or sinusoid is a mathematical function that describes a
  smooth repetitive oscillation.

The formula for the Sine wave is,

A = Amplitude of the Wave
ω = the angular frequency, specifies how many oscillations occur in a unit time interval, in radians per second
φ, the phase,
t = ?

Here ω, is the angular frequency i.e ,

It defines how many cycles of the oscillations are there.

Now let see the frequency,

  Frequency is the number of occurrences of a repeating event per unit
  time

For example, if 100 events occur within 15 seconds the frequency is:

Events = 100;
Time = 15;
FREQ = Events/Time;

Which means there is total 100 cycles of wave in 15 seconds , am i correct ?

Events = 71;
Time = 15;
FREQ = Events/Time;
W = 2*PI*FREQ;   % 2*PI for circular path
A = 2 ;

Now what is the use of small 't' ? in the Sine wave equation , what time is it ?

In the context that you described, $t$ is a variable indicating time. That is, the wave will take on a value of $y(t)$ at the time instant $t$. If this is measured in seconds, then $\omega$ specifies the number of radians that the wave passes through in one second. And, as you noted, if the wave has a frequency of $\frac{100}{15}\ \text{Hz}$, then it will pass through 100 periods in 15 seconds.


Does waveform-based ASR with lowpass filtering preprocessing still require the assumption of FFT? I read some recent papers about time-domain speech recognition, and wondering if it steers away from the stationary assumption of FFT. However, some preprocessing operations like pre-emphasis and low-pass filtering still employ FFT. Does it mean waveform-based ASR with these preprocessing still has the FFT assumption of stationarity?

Many thanks!

Pre-emphasis and low-pass filtering can be done completely in the time domain (using IIR or FIR filtering, for instance), thus without an FFT (or with, by fast convolutions filtering methods).


Common Spatial Pattern for EEG Signals I am having difficulty in understanding the use of CSP for EEG signal feature extraction and subsequently. Since I am using two classes, this query will be restricted to it.

In essence CSP requires 2 classes of trial data, i.e C1, the covariance matrix for class 1, and C2, the covariance matrix for the second class, for calculating W, the spatial filter matrix, using eigenvectors/eigenvalues. Following the procedure the log variance is calculated.
This is useful during the training phase.

However when I conduct a trial, it will belong to either C1 or C2. I cannot calculate the CSP for this particular trial.

How are features of the EEG Signal extracted then?

I am not too sure if this is the correct forum for asking this question. I apologise in advance. Please suggest the correct forum for the same.

I am not sure I understand your question, but for the particular question

  How are features of the EEG Signal extracted then?

CSP gives you a matrix $W$ that represents a spatial transformation of your EEG $X$. When using CSP, $W^TX$ will have the property of maximal variance for one class while minimal variance on the other class. Therefore, when using CSP for classification, the features are the variance of $W^TX$. In practice, the log-variance has a better distribution.

An interesting and brief explanation of CSP can be found on doi:10.1109/MMSP.2010.5662067 (until section II, then it goes on Riemannian geometry which might be out of your scope).


Relationship between $E_b/N_0$ and $\rm SNR$ with channel coding I would like to know what would be the relationship between $E_b/N_0$ and $\rm SNR$ when there is a coding rate of $R_c=K/N$. Normally when there is no channel coding according to this link.

$$\textrm{SNR(dB)} = \frac{E_b}{N_0}(\textrm{dB}) + 10\log10\left(\frac{R_b}{B}\right)$$

However what would happen if I have a linear block code or other coding schemes with code rate $R_c =1/2$ or other values, how to calculate the relationship between $\rm SNR$ and $E_b/N_0$ in this case ?

The relationship between SNR and $E_b/N_0$ is independent of the code rate. Note that $E_b$ is the energy per data bit (not the coded bits), and $R_b$ is the (uncoded) data bit rate. As long as you keep using these values, you can use the formula given in your question. Of course, when going from an uncoded system to a coded system, the values of $E_b$ and $R_b$ will usually change, so the resulting SNR will generally be different.


The current state of the art in audio signal interpolation Three questions:

What are all the metrics one can use to measure audio interpolation quality, objectively? (but also in terms of psychoacoustics if possible)
By those metrics, what is the current state of the art in audio interpolation?
Suppose I were to render two files from a sequence of notes from virtual instruments in two resolutions and then compare an upsampling of one file with the high-frequency rendered version, what software could one use to compare these objectively? - ideally using the before-mentioned metrics

So far, from I've been able to gather, these resamplers provide some of the best quality

http://www.mega-nerd.com/SRC/
http://sox.sourceforge.net/SoX/
http://www.izotope.com/tech/src/

One of the problems that these resamplers seem to have is pre- and post ringing.

I should note that of key interest is signal reconstruction (insofar as that term makes sense), so upsampling more than downsampling.

Edit: Interpolation time efficiency is irrelevant in this context.

Best regards,
The curious :-)

There is this pretty thorough comparison of resampling algorithms: http://src.infinitewave.ca/

You can see the tests they used there. Aliasing is a big one, and is easy to visualize with a spectrogram of a sine sweep. There's also high frequency response -- SRC can roll off high frequencies in addition to allowing them to cause aliasing. You can visualize phase response with the impulse response graph, or with a plot of the phase response.


Halving a floating point time series and it's effect on noise levels Why does halving a floating point time series (audio) cause the noise shown in the spectogram to change?

From having a look, all the steps from getting the first half of X, stft() and amplitude_to_db() output the same values except in the last two frames closest to the halving point (which is expected but not expected to change specshow?).

So why does the graph show an increase in noise throughout the data?

Here is my code to create the two graphs:

audio_file = './nonvoice_snippets_stepvoices/t15_14steps.wav'
X, sample_rate = librosa.load(audio_file)

import matplotlib.pyplot as plt
import librosa.display
plt.figure(figsize=(15, 10))
D = librosa.amplitude_to_db(librosa.stft(X), ref=np.max)
plt.subplot(4, 2, 1)
librosa.display.specshow(D, y_axis='linear')
plt.colorbar(format='%+2.0f dB')
plt.title(audio_file)

#lets get the first half of that.
half = len(X)/2
half_X = X[:half]

plt.figure(figsize=(15, 10))
D = librosa.amplitude_to_db(librosa.stft(half_X), ref=np.max)
plt.subplot(4, 2, 1)
librosa.display.specshow(D, y_axis='linear')
plt.colorbar(format='%+2.0f dB')
plt.title(audio_file)

And the reason why I'm asking this is because I plan to get the Sum of absolute differences (SAD) for use in a Weiner filter to remove the noise you can see in both files. I figured that the first half of the file was a better representation of the noise in the file to get my SAD from but I need to know that this increase in noise won't cause unexpected effects.

Thanks community for being so helpful.

...why does the graph show an increase in noise throughout the data?

It does not. The graph is a depiction. You should focus on the numbers themselves.

Your silence segment includes at last three relatively loud spikes.

Spikes are sharp and brief transitions of the amplitude in the time domain. The bandwidth of a spike is inversely proportional to its duration. That is, the shorter the pulse, the more sinusoids it takes to synthesize it. A quick way to show this is via the Discrete Fourier Transform of the sinc pulse.

Spikes and more generally impulse-like waveforms, manifest themselves in a spectrogram as vertical "line" structures. There are three, relatively clear vertical "lines" in your spectrogram towards the end of the first plot.

The point is that the spike is a sharp transition that stands out from its surrounding samples (in the time domain) and causes relatively large sums in the frequency domain.

This is what is picked up by your normalisation by np.max and is made worse by your choice of y_axis='linear'.

If you try to depict this sequence of numbers x=[8,3,6,4,7,6,120,5,3,6,5,3,2] then the visualisation algorithm will assign the colour at the top of the color map to 120. Effectively, you would get one really bright value which would "swallow" all the rest. If you were to depict the same sequence of numbers but instead of 120 you substituted with 6.5, then you would observe a much more "colourful" output.

This is what is happening in this case. When you plot the whole file, you include the pulses that "skew" your visualisation (but not the underlying numbers) but this is not happening when you simply use the first half of the recording.

To deal with this, you can either use a simple logarithmic transform like log(1+x) or simply switch your scale to logarithmic (i.e. Decibel).

Both will have the same effect: Large values will remain large but small values will be "boosted". This does not distort the number. It simply changes the way the number is presented.

This looks like the beginning of a vinyl recording which includes a bit of hiss along with the occasional pop, because of depositions on the record (?)

Hope this helps.


Difference between DC component and zero frequency component of signal We know that Fourier Transform of a signal exists if it is absolutely integrable
and it exists for periodic signals if impulse functions are allowed.

If we consider the fourier transform of $\text{rect}(t)$ , we get $\text{sinc}(f)$ in frequency domain. Which has got zero frequency component as $1$. But we all know, DC value of $\text{rect}(t)$ is zero.

My question is:

If a signal has got a zero frequency component  in frequency domain ,There must be DC value in time domain. But why there is no DC value in case of $\text{rect}(t)$ in time domain?
What is the difference between DC component and zero frequency component?

Let's first have a look at the rectangular signal given as an example in your question. If you have a rectangle $s(t)$ in the time domain which is $1$ in the interval $[-T/2,T/2]$ and zero elsewhere, its Fourier transform is $S(f)=T\text{sinc}(Tf)$, where I use $\text{sinc}(x)=\sin(\pi x)/(\pi x)$. The value of its Fourier transform at $f=0$ equals $S(0)=T$, which corresponds to

$$\int_{-\infty}^{\infty}s(t)dt=T\tag{1}$$

Its time average (or mean, or DC value) is given by

$$\bar{s}=\lim_{T_0\rightarrow\infty}\frac{1}{T_0}\int_{-T_0/2}^{T_0/2}s(t)dt=0\tag{2}$$

It is clear that any function for which the integral in (1) is finite, must have a DC-value of zero. The integral in (1) is the value of the Fourier transform of the signal at DC, and this is probably what confuses you. The DC value of a signal, and the value of its Fourier transform at DC are not the same. Any signal with a finite Fourier transform at DC has a DC value of zero, i.e. $\bar{s}=0$. Any signal with a non-zero DC value $\bar{s}\neq 0$ has a Dirac delta impulse component in its Fourier transform at DC.

If you write a signal as

$$s(t)=\bar{s}+\tilde{s}(t)$$

where $\bar{s}$ is the DC component as computed from (2), and, consequently, $\tilde{s}(t)$ has a DC component of zero, then its Fourier transform is

$$S(f)=\bar{s}\delta(f)+\tilde{S}(f)$$

where $\tilde{S}(0)$ is finite.

EDIT:
Also note that when the Fourier transform of a signal $s(t)$ has a certain non-zero value at a frequency $f_0$, then this does not entail that the signal has a pure sinusoidal component at that frequency. The same is true for DC. If the Fourier transform has a finite value at DC, the time-domain signal has no DC component, otherwise there would be a Dirac impulse at $f=0$, just as there would be a Dirac impulse at $f_0$ if the signal contained a sinusoid at the frequency.


Does the Wiener–Khinchin-Einstein theorem hold for non-Gaussian processes? If so are there any assumptions? The Wiener–Khinchin-Einstein theorem states that the auto-correlation $(r_{xx}(\tau))$ and spectral density $(S(f))$ are Fourier duals, i.e.

$$r_{xx}(\tau) = \int^{+\infty}_{-\infty} S(f) \exp\left( 2\pi i \tau f \right)\:df$$

This relationship has several assumptions:

-The process must be stationary, i.e the spatial correlation depend only upon the the distance between two points ($\tau$) and not the orientation ($\vec{\tau}$). For Gaussian processes this requires the first and second order moments to be stationary.

My question is:

Does this relationship hold for non-Gaussian processes?
And if so are there any assumptions or restrictions in place regarding its use? For if it holds for non-Gaussian processes, is it suffient for the process Weak or wide-sense stationarity, i.e. only the first and second order moment need to be stationary and higher order moments can be non-stationary, or are there more strict conditions? My gut feeling is the former since the equation requires no-higher order information, however I'm not sure if that is the case.

first sentence of Wikipedia article:

  In applied mathematics, the Wiener–Khinchin theorem, also known as the Wiener–Khintchine theorem and sometimes as the Wiener–Khinchin–Einstein theorem or the Khinchin–Kolmogorov theorem, states that the autocorrelation function of a wide-sense-stationary random process has a spectral decomposition given by the power spectrum of that process.

(emphasis by me)

So this has nothing to do with Gaussian or not, this has to do with stationarity of your process – a process that is Gaussian for every single point in time isn't necessarily wide-sense stationary.

For example, if you have a random signal with an amplitude that is normal distributed for every instant, but whose variance decreases exponentially with time, you can't apply the Wiener-Khinchin theorem.

If you, on the other hand, have a process that actually only takes discrete values $-1$ and $1$, but does so with constant probability, then you can apply the Wiener-Khinchin theorem.

So:

  Does this relationship hold for non-Gaussian processes?

Yes.

  And if so are there any assumptions or restrictions in place regarding its use?

The actual restriction is wide-sense stationarity, as you've noticed:

  For if it holds for non-Gaussian processes, is it suffient for the process Weak or wide-sense stationarity, i.e. only the first and second order moment need to be stationary and higher order moments can be non-stationary, or are there more strict conditions?

Um, that's not exactly the definition of WSS that I have in mind:

Weak sense stationary means that the autocovariance only depends on the time difference, not on the absolute times. This also means that the autocovariance of zero shift, ie. $E[x(t)x(t)]-\mu(t)\mu(t)$ is constant, which can, by considering that these are two expressions for energy, only happen if the first moment $\mu$ is constant.


The relation of sub-carrier and FDMA In Mengyu Liu, Yuan Liu, "Charge-then-Forward: Wireless Powered Communication for Multiuser Relay Networks", it is said that

  while information signals are conveyed by using FDMA over $N$
  sub-carriers (SCs) in the next phase of WIT(wireless information
  transmission). For information transmission, we deﬁne a binary SC
  allocation variable $x_{k,n}$ with $x_{k,n} = 1 $ representing that SC
  $n$ is allocated to pair $k$ for WIT and $x_{k,n} = 0 $ otherwise. Each
  SC is allocated to at most one pair at slot 1 for WIT to avoid
  interference

Does it mean that if the source transmits information to relay, it needs a carrier to help it send the message, so if I have $K$ source, then I have $K$ subcarrier?

If so, why does it still need two different variables? And I can't understand that if my thinking is right, why does it still have two situations, $x_{k,n}=0$ and $x_{k,n}=1$, for that?

Does it mean that if the source transmits information to relay, it needs a carrier to help it send the message, so if I have K source, then I have K subcarrier?

There's $K$ pairs of nodes, if you read the text closely. That means you can have, for example, one source with four information sinks, or four sources and three sinks, or…


Do high frequency oscillations use more than one bandwidth Please excuse me for asking such a basic question. I was having a look at how AM radio works. So there is a carrier frequency and a higher frequency inside that that carries the data.

Does this mean that there would be two frequencies, one high and one low?

What about for stuff like 3G/GSM when stuff is sent over EM would it create noise at other frequencies, say out of their licensed spectrum?

If not how do they squeeze it all into such a tiny band of frequencies?

A modulated radio signal has non-zero bandwidth.  The modulation signal that carries the information has a much lower max frequency compared with the carrier frequency, keeping the total bandwidth (carrier plus both modulation sidebands) of an AM signal within FCC legal limits.

Other modulation schemes are used for GSM/3G, such as OFDM, that create more widely spaced sidebands, but the carrier frequency is so much higher (GHz range) than the modulation sideband width, that the resulting signal bandwidth is only a tiny fractional slice of the total RF spectrum.


Decimation Question about using a low pass filter before downsampling I have a bunch of audio files that have a sample rate of 380kHz, however I would like to reduce this sampling rate to 44.1kHz. From doing some research online, it seems like I will have to use a low pass filter in order to avoid aliasing before I downsample.
I am using python and thinking of using the butter_lowpass_filter. However, I am confused as to what my cutting frequency would need to be and also how much I will need to decrease my sampling rate in the downsampling step.
Any help would be appreciated!

Sample rate conversation is easy in theory but tricky in practice.
Assuming you want to convert to the standard rate of 44.1 kHz (not 44 kHz), you have an awkward conversion ratio. $3800 =2^3 \cdot 5^2 \cdot 19$ and $441 = 3^2 \cdot 7^2$  are mutually prime that means that rational sample rate conversion is impractical,so you need irrational sample rate conversion. If you do not need 44.1 kHz exactly you can use $69/8  \rightarrow 44058Hz$ or $112/13 \rightarrow 44107Hz $
In either case you need a low pass filter. The choice of the filter depends A LOT on your specific application. Specifically

How much energy is above 20 kHz. There really shouldn't be a lot but you don't know before you look. If there is a lot above 40 kHz you can perhaps pre-filter this out before doing the actual conversion
What amount of aliasing and signal to noise can you tolerate.
Do you care more about preserving transients in the time domain or an accurate phase spectrum
Up to which frequency do you have to meet your requirements. Designing a converter that's good up to 20 kHz is a LOT harder than designing one that's good up to 18 kHz. Don't even try to go up to 22 kHz.
Do you care about latency, memory consumption and/or execution time.

The "art" of building a good sample rate converter is to really understand your requirements and trade-offs and then optimize the filter design for your specific case.
A bunch of "generic" hints for audio

Polyphase FIR filter is typically a good choice and a good trade off between time smearing and phase preservation
32 taps, 64 phases gets you around 90-100 dB SNR at 18 kHz and much better below.
If you need to do irrational sample rate conversion, do linear interpolation between the phases.
You have a conversation factor that is "far away" from unity. In this case, it's probably best to convert this in multiple steps and not in one go. Maybe from 380->95 (down by 4) and then from 95->44.1


MFCC extraction for speaker verification Regarding extraction of MFCC from speech signal...

I have checked for the detailed procedure of extraction of MFCC from http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/

For each frame they have extracted the features..
If the signal consists of 100 frames..
then each frame will give 39 values as such..
there by we will get a 39 X 100 matrix for a single speech signal

My question is how to get 39 MFCC coefficients for the whole single signal?
People are saying that they have used 39 MFCCs for a signal, How?

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

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

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

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


Can I assume system is LTI when given by DTFT of impulse response I'm having hard time to grasp it probably because i don't fully understand it.
I understand that when a system is given by $h(t)$ (in general $h(t-\tau)$) i can assume that it is a LTI system.
So i guess my gap of knowledge is that im not sure if it is possible to FT or DTFT impulse response of a non LTI system $h(t;\tau)$. I assume not because the transform depends on the time variable but here i have a 2 variables function and im not really sure how the transforms behave under these conditions.
So like the title says, if a system is given by $H(e^{j\theta})$ can i assume that it is a LTI system? and why is that?
Thanks in advance.

A 1D LTI system is completely characterized by the function $h(t)=T\{\delta(t)\}$ which is denoted as the impulse response of the system. Given an LTI system with impulse response $h(t)$ you can then find its frequency response as $H(j\omega)=\mathcal{F}\{h(t)\}$ where $\mathcal{F}$ stands for continuous time Fourier transform.

When a system is not LTI but, for example, linear time-varying, then it's characterization is possible with the 2D function $h(t,\tau)=T\{\delta(t-\tau)\}$ which is the response of the system to an impulse applied at time $t=\tau$

When such is the case, the Fourier transform can be aplied in the following sense, treating $\tau$ as fixed (or as a parameter) the Fourier transform $\mathcal{F}\{h(t)_{\tau}\} = H_{\tau}(j\omega)$. So in effect you have to consider 1D CTFT for each $\tau$. Whether this is a useful representation or not, or there are other more useful representations for such systems, requires an indepth analysis of time-varying systems.


FIR filter coefficients from frequency responses I have frequency responses of FIR filter of order 13 from obtained 50 frequency samples. How can I get FIR filter coefficients of order 13?

For example, in MATLAB a lot of possible ways I tried ;

I have a frequency samples : w=[0 1] (1x50)

I have frequency responses of FIR filter of order 13: hd (1x50)

b1=fir1(12,wc) % b1:FIR filter ideal coefficients of order 13

hd=freqz(b1,1,w) %hd :FIR filter frequency responses with respect to b1

When I apply inverse fourier transform to FIR filter frequency responses ( ifft(hd,13) ), I don't reach b1. I am supposed to reach b1 coefficients in theory. Am I right?
I apply inverse fft shift for moving to true domain before apply Fourier Transform( ifft( ifftshift(hdg),13) ) but also I don't reach b1. So what else can I do?

There is incompatible things(may be between frequencies and transforms ) I don't understand.

The frequency argument for freqz() should be in radians/second, not normalized frequency in [0,1). It has a form that makes it easy for you:

% length of frequency response vector
N = 50;
% cutoff frequency
wc = 0.25;
% design 12th-order lowpass filter
b1 = fir1(12, wc);
% calculate its frequency response at `N` points around unit circle
hd = freqz(b1, 1, N, 'whole');
% inverse transform the frequency response
b2 = ifft(hd);
% `b1` and `b2` should be the same.
plot(1:length(b1), b1, 1:length(b2), b2);

This yields the plot:

You only see one trace because they are right on top of one another, as they should be. b2 is longer in duration because it has the same length as the frequency response vector, which was arbitrarily chosen to have length 50.


FIR Filters - Type 3 So I understand that a type 3 filter is not suitable for a highpass filter design, but is there any reason why it isnt suitable for a lowpass filter?

So ultimately, can a type 3 linear-phase FIR filter be used to design a lowpass filter? Why or why not?

A type 3 FIR filter has odd symmetry and an odd number of taps. For this reason it has frequency response zeros at $\omega=0$ (DC) and $\omega=\pi$ (Nyquist), corresponding to transfer function zeros at $z=1$ and $z=-1$:

$$H(1)=\sum_{n=-M}^{M}h[n]=0\\
H(-1)=\sum_{n=-M}^{M}(-1)^nh[n]=0$$

where $N=2M+1$ is the filter length. So it can neither be used as a high pass nor as a low pass filter.

Apart from that, the phase shift of $\pi/2$ caused by the odd symmetry is usually undesirable for frequency-selective filters. You could use such a filter for implementing Hilbert transformers or differentiators.


How to Generate White Gaussian Noise with Known PSD in MATLAB I want to generate white Gaussian noise with known PSD. Using wgn I can generate noise with specific power (dBW) but I want to generate a white noise with constant or know PSD.

Can anyone give me some idea.
Thanks
Franklin

The specific PSD is basically white noise which was filtered with a filter which its Magnitude is the same as the PSD.

This is a result of the Wiener Khinchin Theorem.

So if you have the shape of the PSD do the following:

Create a filter which in the Frequency domain has the magnitude (More accurately the squared magnitude $ {\left| H \left( \omega \right) \right|}^{2} $) of the PSD you want (Namely, you don't care about the phase).
Create a White Noise.
Filter the white noise with the filter created in 1.

I would note that in the case of a White Noise the PSD is constant.
Hence in this case all you need is to generate white noise with normalized variance and then to multiply it by a constant (Changing its variance) to the value you're after.

Enjoy...


Convolution of two rectangular pulses intuition From what I understand the convolution of the impulse response of a system with the input to that system gives the output.

Now if the impulse response is a rectangular function and the input is also a rectangular function we get a triangular function as the output.

What I dont understand is: How could this ever be the case?
As the rectangular function of the impulse response has only two levels as does the input. How can the system ever produce a value that is between these levels let alone a ramp of values between these levels.

It's obvious I am missing something but I cant quite put my finger on what that is.

A convolution integral is an overlap integral, i.e., for any given shift of the two aperiodic functions being convolved, the convolution integral is simply the overlap area. McGillem and Cooper [1, p. 58] defined the convolution integral of $x_1$ and $x_2$ as

$$\mathrm {x_3 =x_1*x_2 =\int_{-\infty}^{\infty}x_1(\lambda)x_2(t-\lambda)\,\mathrm d\lambda \tag{1}}$$

As a simple graphical illustration of the defining integral, they considered the following two rectangular pulses:

With $x_1$ and $x_2$ as shown in the above figure, their convolution is shown in the figure below:

This figure is redrawn from [1, p. 59]. The shaded areas are the overlap areas as a function of the shift, $t$, and the resulting convolution has a trapezoidal shape. If the rectangular pulses had had equal width, then the convolution would havec simplified to an isosceles triangular shape.

1  C.D. McGillem, G.R. Cooper, "Continuous and Discrete Signal and System Analysis", 2nd Ed., Holt, Rinehart and Winston, NY, ©1984, pp. 58-59.


Establishing Phrase Frequency in Audio Files Is there any Application, Programming Library or Theory that could be used to identify human language phrases in a set of audio files, so as to calculate the frequency of occurrence of each such phrase. For example, how often does the phrase "Thank You" occur in these files?

One solution might be to convert the audio to an International Phonetic Alphabet form but I cant find such a tool.

I'm basically trying to find patterns in speech and the frequency at which they occur. Of course the pattern for "Thank You" is slightly different in different dialects, but I've no problem in having different matches in such cases.

This is not a simple problem.  Almost certainly any "home-brewed" system you create looking for "frequency patterns" is going to have horrible performance, unless it develops complexity comparable to full-fledged speech recognizer.  The reason for this is that recognizing speech, while it seems like a simple task to our ears, is in fact quite complex, and requires large amounts of context to have good accuracy.  What I'm trying to say is that if you just try to recognize the phonemes from the audio locally, you're going to get gibberish. You need to take into account the neighboring phonemes, the vocabulary of the language being spoken, and the large scale structure and content of the surrounding sentences in order to get anywhere near human level performance.

At the end of the day if you try to write something which "recognizes frequency patterns of speech" from scratch, you're going to be reinventing the wheel, and a pretty dang big wheel at that.

For this reason, the best thing for you to do would be to take an existing speech recognition system, and modify its output to suit your needs.  I would suggest using CMU Sphinx to decode the audio into an N-Best list, and then search this list for the text you are looking for.  CMU Sphinx has good recognition performance and is highly configurable.


Computational burden of EMD/Huang-Hilbert vs wavelet I am working on an online application of signal processing and pattern recognition. It involves sampling the signal at 2 MS/s, extracting features and classifying.

My classifier has pretty good accuracy so far using the DWT (discrete wavelet transform) with 3 decomposition levels. However, I have been studying the EMD (Empirical Mode Decomposition) and trying to apply it to the problem since my data is highly non-stationary. The algorithms that I got are waaaaaay slower than the wavelet and PSD (power spectral density) estimation by Welch's method (another one that I was using.)

Is this normal/expected behavior? For wideband data is possible to implement it to work online?

Sorry about the broad question, it seems that it is more a curiosity for discussion that anything.

There are relationships between fixed wavelets and EMD, for certain classes of signals, cf. Empirical Mode Decompositions as data-driven wavelet-like wavelet-like expansions. A related relation between a fixed basis and an adaptive one can be found between Fourier transforms and Karhunen-Loève/PCA type decompositions.

Of course, when the vectors used for representations are learnt instead of being fixed beforehand, they are likely to require more computations. And are likely to be more difficult to implement online. Plus, DWT are often critically sampled, while EMD are redundant (hence greedier). Yet, some papers (see below) evaluate it proportionally to the burden of the FFT. Roughly to that of a redundant stationary wavelet transform.

Then, for non constrained versions, EMD can be quite sensitive to border effect, local spike instability, etc. As a result, from one signal to a similar other, a different number of IMFs can be necessary, which is not the case of a fixed-level DWT (but you can find that with adaptive wavelet packets).

Further, each can be sensitive to the features you can extract.

So yes, it is not so surprising to me, with experience on seismic data and combustion signals with sporadic instabilities. Yet, it can work on some data.

Finally, some online implementations exist:

On-line empirical mode decomposition biomedical microprocessor for Hilbert Huang transform, 2011

  On-chip implementation of Hilbert-Huang transform (HHT) has great
  impact to analyze the non-linear and non-stationary biomedical signals
  on wearable or implantable sensors for the real-time applications.
  Empirical mode decomposition (EMD) is the key component for the HHT
  processor. In tradition, EMD is usually performed after the collection
  of a large window of signals, and the long latency may not be feasible
  for the real-time applications. In this work, the architecture of
  on-line EMD for biomedical signals is proposed. The on-line
  interpolation method with data reuse as well as component and
  iteration loop decomposition is applied to obtain low latency and low
  hardware cost. The first chip of EMD processor is fabricated in UMC
  90nm LL process and consumes 57.3µW.

Online Empirical Mode Decomposition, 2017

  The success of Empirical Mode Decomposition (EMD) resides in its
  practical approach to dissect non-stationary data. EMD repetitively
  goes through the entire data span to iteratively extract Intrinsic
  Mode Functions (IMFs). This approach, however, is not suitable for
  data stream as the entire data set has to be reconsidered every time a
  new point is added. To overcome this, we propose Online EMD, an
  algorithm that extracts IMFs on the fly. The two key elements of
  Online EMD are a sliding window to compute local IMFs, and a stitching
  procedure to gradually append local IMFs to the final result. Using
  synthetic data we show that the decomposition quality of Online EMD is
  similar to classical EMD. We also present results obtained with a real
  data set to expose the practical advantages of Online EMD when dealing
  with data stream or large data set.

On the computational complexity of the empirical mode decomposition algorithm, 2014

  It has been claimed that the empirical mode decomposition (EMD) and
  its improved version the ensemble EMD (EEMD) are computation
  intensive. In this study we will prove that the time complexity of the
  EMD/EEMD, which has never been analyzed before, is actually equivalent
  to that of the Fourier Transform. Numerical examples are presented to
  verify that EMD/EEMD is, in fact, a computationally efficient method.


Understanding noise removal method using wavelets I am trying to understand how wavelet transform can be used to denoise a time series or signal and how to plot the scalogram image. My signal has a lot of fluctuations and as such I am finding it difficult to denoise. Morevoer, to plot the scalogram I need to know the frequency. I don't know what is the frequency for this particular kind of time series obtained from a dynamical system of the form: Logistic Map given by:
$$x[n] = 4\big(x[n-1]\big)\big(1-x[n-1]\big)$$
Systems similar to this type of dynamical systems are the Lorenz, Mackey-Glass. Can somebody please help:

How to properly denoise the signal? As observed, from the plot the denoised signal denoised does not look exactly the same as the clean signal x (black dotted line), so what are other parameters or wavelet types that I could use and how to decide which ones to use. Is there a rule of thumb?

What is the sampling and nyquist frequency for this kind of signal and

how to plot the scalogram image: I used wt() to obtain the wavelet coefficients. After that how to plot the image of scalogram so that X axis is time and Y axis is Frequency?
 x(1) = 0.1; % initial condition (can be anything from 0 to 1)
 M = 50; %number of data points (length of the time series)
 for n = 2:M, % iterate
     x(n) = 4*x(n-1)*(1-x(n-1));
 end

 %add noise
 x_noise = awgn(x,10,'measured');

 %denoise using wavelet
 denoised = wdenoise(x_noise, 3,'Wavelet','db3',...
     'DenoisingMethod','Bayes',...
     'ThresholdRule','Median',...
     'NoiseEstimate','LevelIndependent');
 figure
 plot(x_noise)
 axis tight
 hold on
 plot(denoised,'r')

 fb = cwtfilterbank('SignalLength',M);
 [cfs,frq] = wt(fb,denoised);

Your signal (with initial par x0 =0.1) is already noise like and high frequency. It will be hard to distinguish it from the added white noise... One thing you can do is to interpolate (resample) the time series by a large enough factor and then later add the white noise. This will artifically help to separate the noise spectrum and your signal spectrum but the signal lengths will also be increased. whether it is what you have to do is up to you !
The following modification apparently improves the noise removal, but fundamentlaly it's separating the noise spectrum from the signal. So whether this is a viable option is upto your applications:
 M = 50;        % number of data points (length of the time series)

 x(1) = 0.5;    % initial condition (can be anything from 0 to 1)
 for n = 2:M,   % iterate
     x(n) = 4*x(n-1)*(1-x(n-1));
 end

 U = 10;       % interpolation factor
 xU = resample(x,U,1);  % just interpolate the obtained sequence

 % add noise onto the interpolated sequence xU
 x_noise = awgn(xU , 10 , 'measured');

 %denoise using wavelet
 denoised = wdenoise(x_noise, 3,'Wavelet','db3',...
     'DenoisingMethod','Bayes',...
     'ThresholdRule','Median',...
     'NoiseEstimate','LevelIndependent');

 denoised = resample(denoised,1,U);   % downsample de-noised sequence back

 figure
 plot(x_noise(1:10:end))   % down-sample noisy seqeunce on the fly for displaying
 axis tight
 hold on
 plot(denoised,'r')
 plot(x,'c--');
 legend('noisy','denoised','clean');

 fb = cwtfilterbank('SignalLength',M);
 [cfs,frq] = wt(fb,denoised);

The result looks like :


How do I scale Gaussians for Edge detection Say I have an $m\times n$ image and I want to use DOG for edge detection. I can see in this answer he is using:

gaussian1 = fspecial('Gaussian', 21, 15);
gaussian2 = fspecial('Gaussian', 21, 20);

I was working with this image for reference, and I tried to move the size from 1 to 30. I cant make sense of which output is better. I think the sharpest edges I received were for hsize=6, but how can I tell?

After following the advice on the comments and reading the wiki entery's abstract on the subject I am still puzzled. It is my understanding that the band pass parameters is dependent in the image and there is no efficient automatic way to derive the sigma and size parameters of the filters from the image.Does this means I have to manually adjust the filter sizes for any new image?

Is there an automatic way to do so? (I am sure there is since this algorithms work independently on the market)
say I want to generate response for different filters as an experiment:

Which iterations should I go through?
When I am looking at the result image, what should I look for to ssay it is the desired response?

Does this means I have to manually adjust the filter sizes for any new image?

yes! You need to know the kind of edges you want to detect. And that "kind of edge" is a "width of the edge in pixels" in disguise. So, you need more "intelligence" than just a DOG if you want to detect your favourite kind of edge in pictures of any zoom / resolution.

  Is there an automatic way to do so?

Plenty! You could, for example, simply try out some parametrizations for your DOG and automatedly select the parametrization that leads to an edge image that fits your signal model best – for example, if you want to find square boxes in images, looking for the "cleanest" straight lines.

Again, all depends on what you want to detect, and how you model your signals. There's no "universally" best edge detector. (or else, that would be the only one everyone uses, right?)