Connecting edges detected by an edge detector I have a binary image obtained from a canny edge detector. The edges are not nicely detected in the center and I need to join them. The connection of edges is orientation and neighborhood dependent. I need to connect the vertices if they are below a threshold (say less than 5 pixel distance). If my almost linear features are nicely oriented, this threshold should be relaxed. (This is the best case.)

I tried the Hough transform, but it didn't work for me because I don't have straight lines. Dilation and erosion are also not good; they make the images messier.

The approach I am trying is to first detect vertices and nodes (in MATLAB with bwmorph) and then make leaves as a separate feature. This is done as follows:

Scanning in a 3x3 window, look for neighbors.
Traverse the complete connected object.
Try to fit a line (or maybe a quadratic polynomial).
Check feature-by-feature if it's worth a connection or not.

The implementation is not straightforward since the decision-making part where vertices need to be connected is tricky.

This might not be complete solution, but will give you good direction.

Basically, what is the key criteria of to say that edges match? That "locally" the gradient of the edge matches and to some extent the distances are reasonable against how long the edge is continuous.

If you have geometric edges, like long straight lines, Hough will do very seamless job right away. But this doesn't work when edges are arbitrary curve. In this case, you can still think of curve as roughly piece wise constant segments (good enough in your case), hence, you should take Hough locally. i.e. you can take a small portion of image (say a block) compute the Hough and identify some peaks. Based on this, you can identify that de-houghing is creating meaningful gap, if so, keep it or move on.

Once smaller gaps are filled up, you can extend the same to take up on a larger size. the peaks will be more, but you may select fewer from it.


Fourier Transform Identities We know the below,

$$
\mathscr{F}\big\{x(t)\big\}=X(f) \tag{1}
$$
$$
\mathscr{F}\big\{x(-t)\big\}=X(-f) \tag{2}
$$
$$
\mathscr{F}\big\{x^*(t)\big\}=X^*(-f) \tag{3}
$$

Now, if for some signal

$$
x(-t)=x^*(t) \tag{4}
$$

Then, is it safe to assume the following?

$$
X(-f)=X^*(-f) \tag{5}
$$

or does it depend on the type of signal?

You are correct. Your last equation is simply a fancy way of saying that $X(f)$ is real valued.

In general: if it's real in one domain, it's conjugate symmetric in the other.


Need for Oversampling in OFDM and symbol time determination I am trying to implement OFDM on gnuradio. I have two questions.

The sampling frequency which is slightly higher than available bandwidth  satisfies Nyqvist criterion for complex baseband signals. So why do we need to add oversampling factor again? Suppose say oversampling factor is 2. How does that help us in recovery of input signal?
I have seen in an OFDM implementation that

$$
t_{\rm symbol}=t_{\rm sampl} *(N_c+N_G)*{\rm oversamplingrate }
$$

where

$t_{\rm symbol}$ = Symbol duration,

$t_{\rm sampl}$=sampling duration,

$N_c$ = No. of subcarriers consisting of data and pilot

$N_G$ = No. of guard band taps.

How does this formula come about? I also have a doubt as to whether each subcarrier corresponds to a single symbol or you can have multiple symbols per subcarrier. What is parameter that decides this?

As the questioner refers to my code at github in the comments, I'd like to clarify.

There is no need to apply oversampling in OFDM, nor is oversampling required to recover input signals. However, in practical implementations you have to consider some additional constraints regarding the signal's bandwidth. The WARP v2 radio board is clocked with 40 MHz, i.e., the Nyquist frequency is 20 MHz. However, the WARP v2 architecture limits the generated baseband signal's bandwidth to ~10 MHz by TX and RX lowpass filters. This is a meaningful measure to avoid unpleasant aliasing effects. Hence, the OFDM signal's bandwidth should be limited to 10 MHz in the baseband. You can easily achieve this by padding all subcarriers above 10 MHz with zeros. This will work perfectly well.
An alternative (but probably rather uncommon) way to reduce the bandwidth is to apply 2x oversampling while using the majority of all subcarriers. I implemented oversampling in my generic (not-OFDM-related) transmission procedure and finally enabled it in my OFDM code as I could achieve a slightly higher SNR in our office environment with infrequent interference and some narrowband interference. However, you can safely disable oversampling in my code and apply zero-padding instead.
An OFDM symbol is the time-domain representation of the symbols concurrently transmitted over all subcarriers. It consists of just as many samples as there are OFDM subcarriers since it is calculated as the IDFT of all subcarriers' symbols. When adding the cyclic prefix, you add $N_G$ samples. The duration in seconds is the number of samples times the time per sample. When oversampling, the symbol duration extends according to the oversampling factor.

Kind regards,
Robin


Discrete signal testing for periodicity How would one go about determining if the following discrete time signal x[n] is periodic, and if it is, determine its fundamental period?

I understand that the period for the second exponential term is 6, but apart from that I am unable to further my calculation.

The answer to the question above is:

but I fail to understand how the equation highlighted in yellow comes to be. Any help with this matter would be appreciated.

The first term is just a constant, so it's not relevant to periodicity. The sequence as written is periodic with a fundamental frequency of 1/6.

Your cited solution makes no sense to me.


best implementation of a real-time, fixed-point iir filter with constant coefficients I have implemented IIR filter in direct form, parallel form, and cascade form. The input in all the cases is a kronecker delta function. How to judge which implementation is better is in real time and why? What are the factors to be considered while selecting the best implementation? Please help.

what is fraction saving? can you write a code.so that i can understand more clearly?

Let's call the quantizer operator $\operatorname{Quant}\{\cdot\}$ .  So the output of the quantizer, with $v[n]$ going in, is

$$y[n]  = \operatorname{Quant}\{ v[n] \}$$

which we shall model as an additive error source:

$$y[n]  =  v[n] + q[n]$$

No matter how the quantizer works, we can always derive $q[n]$ in the run-time system with

$$q[n]  \triangleq  y[n] - v[n]  =  \operatorname{Quant}\{ v[n] \}  -  v[n]$$

Even though $q[n]$ is, strictly speaking, a deterministic function of $v[n]$, often, if the signal swings are much larger than the step size, $\Delta$, of the quantizer, we model $q[n]$ as some kinda random number and a "noise source".  It's when the signal gets very small, on par with $\Delta$, that this assumption is no good, and that's when you need to consider adding dither before the quantization operation.  I doubt you will need dither for your situation, but "ya never know".

The simplest quantizer to implement in a fixed-point 2's complement system is simply dropping the less significant bits that are to the right of the rounding point or "binary point":

$$ \operatorname{Quant}\{ v[n] \}  \triangleq  \Delta \left\lfloor \frac{v[n]}{\Delta} \right\rfloor $$

$\lfloor \cdot \rfloor$ is the floor function, the same as the what the C function floor(): returning the most-positive integer that does not exceed the argument to the function.  $\lfloor u \rfloor$ is the sole integer that satisfies:

$$  u - 1 < \lfloor u \rfloor  \le  u  $$

or

$$ \lfloor u \rfloor  \le  u  <  \lfloor u \rfloor + 1 $$

If $u \ge 0$, often we say that $\lfloor u \rfloor$ is the "integer part" of the real value $u$.  The "fractional part" of $u$ is $u - \lfloor u \rfloor$.  It is clear that the fractional part of $u$ is

$$ 0  \le  u - \lfloor u \rfloor  <  1 $$

So the quantization error is the fractional part remaining after quantization is

$$ q[n]  =  \Delta \left\lfloor \frac{v[n]}{\Delta} \right\rfloor -  v[n] $$

That is the additive noise signal in terms of the input to the quantizer.

$$ -1  <  \frac{q[n]}{\Delta}  =  \left\lfloor \frac{v[n]}{\Delta} \right\rfloor - \frac{v[n]}{\Delta}  \le  0 $$

So we know that the additive quantization error, $q[n]$, has a range of

$$ -\Delta < q[n] \le 0 $$

Now, we assume (that's the only "rule-of-thumb") that essentially the value of $q[n]$ is equally likely in that range so then $q[n]$ is viewed as a random number with uniform p.d.f.:

$$ p_q(\alpha) = \frac{1}{\Delta} \operatorname{rect}\left(\frac{\alpha}{\Delta} + \frac{1}{2}\right) $$

where $\operatorname{rect}(\cdot)$ is the standard unit rectangular function:

$$\operatorname{rect}(u) \triangleq \begin{cases}
1, & \text{if } \ -\frac{1}{2} < u \le \frac{1}{2} \\
0, & \text{otherwise}
\end{cases}$$

The expectation value or mean of $q[n]$ is $-\frac{\Delta}{2}$ (a half-step bias toward the negative direction) and the variance is $\frac{\Delta^2}{12}$.  This is well known about uniform quantization.

In the frequency domain:

$$Y(z)  =  V(z)  +  Q(z)$$

(for frequency response, we substitute $z = e^{j 2 \pi f T}$  where $T$ is the sampling period and in reciprocal units as $f$. Normalized frequency (using the DTFT or "discrete-time Fourier transform") is $\omega \triangleq 2 \pi f T$ and ranges $-\pi < \omega < \pi$.  In the DTFT, if the quantization error really is random, $Q(e^{j \omega})$ would have a Dirac impulse at $\omega=0$ of magnitude $-\frac{\Delta}{2}$ (to represent the DC bias of $-\frac{1}{2}$ LSB because we're always rounding down) and otherwise a flat spectrum from $-\pi$ to $\pi$, that would, because of Parseval's theorem, have area of $\frac{\Delta^2}{12}$.  Maybe there's a $2\pi$ factor in that (which is why i don't usually like angular frequency over "ordinary frequency").

we know that:

$$ \int\limits_{-\infty}^{+\infty} p_q(\alpha) \ d\alpha \ = \ 1 $$

The mean of quantization error (DC component):

$$\begin{align}
\mu_q & \triangleq \overline{q[n]} \\
& \triangleq \lim_{N \to +\infty}\frac{1}{2N+1} \sum\limits_{n=-N}^{N} q[n] \\
& = \int\limits_{-\infty}^{+\infty}  p_q(\alpha) \ \alpha \ d\alpha \\
& = Q(e^{j 0}) \\
& = -\frac{\Delta}{2}
\end{align}$$

The mean square of quantization error:

$$\begin{align}
\overline{\left| q[n] \right|^2} & \triangleq \lim_{N \to +\infty}\frac{1}{2N+1} \sum\limits_{n=-N}^{N} |q[n]|^2 \\
& = \int\limits_{-\infty}^{+\infty}  p_q(\alpha) \ \alpha^2 \ d\alpha \\
& = \frac{1}{2\pi}\int\limits_{-\pi}^{+\pi} \left| Q(e^{j \omega}) \right|^2 d\omega \\
& = \frac{\Delta^2}{3}
\end{align}$$

The variance of quantization error (AC power):

$$\begin{align}
\sigma_q^2 & \triangleq \overline{\left| q[n] - \overline{q[n]} \right|^2} \\
& \triangleq \lim_{N \to +\infty}\frac{1}{2N+1} \sum\limits_{n=-N}^{N} |q[n] - \overline{q[n]}|^2 \\
& = \int\limits_{-\infty}^{+\infty}  p_q(\alpha) \ (\alpha-\mu_q)^2 \ d\alpha \\
& = \overline{\left| q[n] \right|^2} - \left|\overline{q[n]}\right|^2 \\
& = \frac{\Delta^2}{12}
\end{align}$$

Now (to get back to "fraction saving") suppose $v[n]$ is derived from your input, $x[n]$, as so:

$$v[n] = x[n] - q[n-1]$$

So whatever error we added to $v[n]$ by quantizing it, we will subtract, before quantization, from the quantized sample immediately following.

(n.b.: Remember that $q[n] \le 0$ and when we "add" $q[n]$ in the quantization process, we are effectively subtracting whatever bits are to the right of the binary point by zeroing those bits, which is always rounding down.  So when we "subtract" $q[n-1]$ from the quantizer input, $x[n]$, before quantization, we are adding those bits we dropped, extended on the left with zeroed bits, not sign extended.)

$$\begin{align}
y[n] & = \operatorname{Quant}\{ v[n] \} \\
     & = \operatorname{Quant}\{ x[n] - q[n-1] \} \\
     & = \quad x[n] - q[n-1] + q[n] \\
     & = \quad x[n] + (q[n] - q[n-1])  \\
\end{align}$$

or, in the frequency domain:

$$\begin{align}
Y(z) & = X(z) \quad + \quad Q(z) - Q(z) z^{-1} \\
     & = X(z) \quad + \quad (1 - z^{-1}) Q(z) \\
     & =  H_x(z) X(z) \quad + \quad H_q(z) Q(z) \\
\end{align}$$

So the transfer function from input $x[n]$ to output $y[n]$ is still just a wire:

$$ H_x(z)  =  1 $$

But the transfer function from the additive quantization noise $q[n]$ to the output $y[n]$ is

$$ H_q(z) =  1 - z^{-1}  =  \frac{z-1}{z} $$

So, like FIR filters, there is a pole at the origin $z=0$ (big hairy deel), but there is a zero right at the DC point, $z=1$.  This kills any contribution of quantization error to the output at the frequency of 0.  Even that $-\frac{\Delta}{2}$ bias, from always rounding down, is killed.

The magnitude of the noise to output frequency response is

$$\begin{align}
|H_q(e^{j \omega})| & = \frac{|e^{j \omega} - 1|}{|e^{j \omega}|} \\
     & = |e^{j \omega/2}| \frac{|e^{j \omega/2} - e^{-j \omega/2}|}{|e^{j \omega}|} \\
     & = 2 \left| \sin\left(\frac{\omega}{2} \right) \right| \\
\end{align}$$

The transfer function, $H_q(\cdot)$, is 0 (or $-\infty$ dB) and when you're up at Nyquist (or $\omega=\pi$) the gain is 2 (or 6 dB).  overall, there is no decrease in quantization noise energy, in fact, because energy is magnitude-squared, there is an increase in total quantization energy, but at frequencies below $\frac{\pi}{2}$, the quantization noise is reduced.  For us audio guys, we like that because none of us hears too well above 11 kHz anyway.

But one thing it helps is that it destroys that limit cycle of IIR filters that gets them stuck on a non-zero DC value, even when the input truly goes to a dead zero.

Code for a biquad (Direct Form I, not Direct Form II which you don't wanna use, especially for fixed-point) that might look like that might be as below. To be clear in sign convention, the transfer function from input to output is:

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}} {1 - a_1 z^{-1} - a_2 z^{-2}} $$

This corresponds to a DF1 difference equation as:

$$ y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] + a_1 y[n-1] + a_2 y[n-2] $$

With quantization it's:

$$ y[n] = \operatorname{Quant}\{ b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] + a_1 y[n-1] + a_2 y[n-2] \} $$

I am sorta assuming that short is 16-bit and long is 32-bit. because, in terms of actual mathematical values, any stable IIR will have  $-2 < a_1 < +2$ and $-1 < a_2 < +1$ ($b_0$, $b_1$, and $b_2$ are relative to each other so we worry less about their fixed range), so let's say that all our coefficients are all scaled up by a factor of $2^{14} = 16384$ (so that the equivalent integer value of the fixed-point number ranges from $-2\cdot2^{14} = -32768 \ <$ a1 $= 2^{14} a_1 < \ +2\cdot2^{14} = +32768$).  so the quantization or rounding point is at bit 14 of the 32 bit number.  The coefficients are also stored as 16-bit short, but we are casting them to long upon loading so that the cast is not continuously necessary in the inner loop.  We want coefficient times state to be a long, even if both multiplier and multiplicand going in are short.

Version 1:

long accumulator;

long a1, a2, b0, b1, b2;   /* 2nd-order IIR coefficients */

short state_x1, state_x2, state_y1, state_y2, state_error;

short x, y;

short get_input(void);

void put_output(short y);

short i, num_samples;
         /* num_samples is the number of samples per block */

/*  load up your coefficients, cast all to long */

/*  load up your states from whereever */

accumulator = (long)state_error;

for (i=0; i<num_samples; i++)
    {
    x = get_input();

    accumulator += b0*x;
    accumulator += b1*state_x1;
    accumulator += b2*state_x2;
    accumulator += a1*state_y1;
    accumulator += a2*state_y2;

    if (accumulator > 0x1FFFFFFF)
       {
       accumulator = 0x1FFFFFFF;     /* clip value */
       }
    if (accumulator < -0x20000000)
       {
       accumulator = -0x20000000;    /* clip value */
       }

    y =  (short)(accumulator>>14);   /* always rounding down */

    state_x2 = state_x1;             /* bump the states over */
    state_x1 = x;
    state_y2 = state_y1;
    state_y1 = y;

    accumulator = accumulator & 0x00003FFF;
    /* keep the fractional bits that you dropped for the next sample
       otherwise clear the accumulator */

    put_output(y);
    }

state_error = (short)accumulator;

/*  save your states back to whereever  */

Now another trick, to speed things up, is to adjust the a1 coefficient by subtracting the scaled fixed-point version of 1 (which, in this case, is $2^{14}$).  Then the code might look like:

Version 2:

long accumulator;

long a1, a2, b0, b1, b2;   /* 2nd-order IIR coefficients */

short state_x1, state_x2, state_y1, state_y2, state_error;

short x, y;

short get_input(void);

void put_output(short y);

short i, num_samples;
         /* num_samples is the number of samples per block */

/*  load up your coefficients, cast all to long */

/*  load up your states from whereever */

a1 = a1 - 0x00004000;    /* adjust a1 coefficient */
/* 0x00004000 is "1" in our scaled-up fixed-point coefficients */

accumulator = (long)state_error + (((long)state_y1)<<14);

for (i=0; i<num_samples; i++)
    {
    x = get_input();

    accumulator += b0*x;
    accumulator += b1*state_x1;
    accumulator += b2*state_x2;
    accumulator += a1*state_y1;
    accumulator += a2*state_y2;

    if (accumulator > 0x1FFFFFFF)
       {
       accumulator = 0x1FFFFFFF;     /* clip value */
       }
    if (accumulator < -0x20000000)
       {
       accumulator = -0x20000000;    /* clip value */
       }

    y = (short)(accumulator>>14);    /* always rounding down */

    state_x2 = state_x1;             /* bump the states over */
    state_x1 = x;
    state_y2 = state_y1;
    state_y1 = y;

    put_output(y);
    }

state_error = (short)(accumulator & 0x00003FFF);

/*  save your states back to whereever  */

Note the missing bit-masking instruction at the bottom of the loop.  So you have the left bits having an extra value of $y[n-1]$ in there.  That's why we have to fix a1 by subtracting 1 (again scaled by $2^{14}$) from it.  That does affect the range of a1 by bumping it over.  Now it's $-2 < a_1 - 1 < +2$ or $-1 < a_1 < +3$ but the range of all of the other coefficients are unchanged.


How to compute the energy of two signals? Say I have two signals, one being $g_1(t)$ and the other one being $g_2(t)$, and say that I am being asked to compute the energy of $g_3(t)$, where $$g_3(t) = g_1(t) + g_2(t)$$

Is it correct to compute the energies of $g_1(t)$ and $g_2(t)$ separately and simply add the two values to find $g_3(t)$ or is it correct to compute them as one value?  I.e. $$\int_ {-T}^{T} \big[g_1(t)+g_2(t)\big]^2 \,dt $$

You simply have to apply the definition of energy. Assuming all signals involved are real, the energy of $g_3(t)$ is given by

\begin{align}
E &= \int_{-\infty}^\infty g_3^2(t) \, dt \\
&= \int_{-\infty}^\infty \big(g_1(t) + g_2(t)\big)^2 \, dt \\
&= \int_{-\infty}^\infty g_1^2(t) \, dt + 2\int_{-\infty}^\infty g_1(t)g_2(t) \, dt + \int_{-\infty}^\infty g_2^2(t) \, dt.
\end{align}

One interesting case is where $g_1(t)$ and $g_2(t)$ are orthogonal, in which case the middle integral is zero. In that case, the energy of $g_3(t)$ is the energy of $g_1(t)$ plus the energy of $g_2(t)$. One simple example is two rectangular pulses that do not overlap.


Equation visualization of complex sinusoid function I'm trying to learn Fourier Transform & Signal

I can visualize how this expression $ A \sin(2 \pi f t)$ could turn into that curvy sinusoidal signal.
But i couldn't visualize how $F(\nu)e^{2 \pi i f t}$ would look like as a signal....

Why the natural exponent suddenly popped up? what does it do?
i understand that $F(\nu)$ gives information of the amplitude and the phase shift, but what sort of information i could get from $e$ ?
Does all signals always decaying/growing by time?

I saw this a lot in Laplace Transform: $e^{-st}$ but don't understand how it would look like as a signal.

The exponent function is directly related to the sine and cosine via Euler equation:

$$e^{i \theta} = \cos(\theta) + i \sin(\theta)$$

For example for function $z=e^{i \pi t}$:

$$z=\cos(\pi t) + i \sin(\pi t) $$

You could plot it in 3D, where x axis is the time, and y & z axes are real and imaginary part of your complex signal, respectively. In 3D that is going to be a helix:

But if you look at projections onto the y=0 and z=0 planes, then you will notice sine:

and cosine waves:

On the other hand, if you look at the x=0 projection, then you will observe a (unit) circle:

You can find the applet with this visualization in Wolfram Demonstration Project.


What is the meaning of spectral contamination in image processing? I resize the Multispectral image to perform IHS-to-RGB transformation. I used the Nearest Neigbour Interpolation method when resizing. While researching this issue, I saw that this is the best option to prevent spectral contamination, but I do not know exactly what the definition of spectral contamination means. I searched, but could not find detailed information on this topic.
What is the meaning of spectral contamination in image processing?

"Spectral contamination" is used for designating a few phenomena connected with spectral analysis in optics, acoustics, analytical chemistry. I have never heard of these words to be used in connection with resizing of digital color images.
The gray-color images can acquire colored pixels in the course of resampling, one user of Adobe Photoshop calls this effect a "color variations".
In computer graphics, the colored edges appear in screen font glyphs rendered with LCD and OLED displays as the result of subpixel rendering. This technique is also used to enhance  the resolution of all image types on layouts which are specifically designed to be compatible with subpixel rendering.
You can find online articles where research is done to detect the spectral contamination with the help of multispectral image cameras, but this is the spectal contamination connected with spectral analysis in optics, acoustics, analytical chemistry and it is not related to the color variations resulting from resampling the digital images.


Why is the output of the LPF half that of the real/imaginary part of the complex envelope? The passband signal $\ u_p(t) =
I_{[−1,1]}(t) \cos 400πt $ is passed through an LPF as $\ u_p(t) \cos(401\pi t).$ The question is to find the output. The example does the following:

If we let $\ u_1 = u_{c1} +ju_{s1}  $ denote the complex envelope with respect to the reference $\ e^{j401t} $,  the output of the LPF is $\ u_{c1}/2.$ Why is it half?

The passband signal can be rewritten as,
$\ u_p(t) = \Re(I_{[−1,1]}(t)e^{−jt}e^{j401t})$ to find $\ u_{c1}$
but how is this justified if we want to find the output of $\ u_p(t) \cos(401\pi t) $ and not just $\ u_p $ that is just rearranged.  Shouldn't the translation induced by the cosine be taken into account. I'd appreciate any help.

The idea of this exercise is to show what happens if you have a frequency offset at the receiver, i.e. if the actual carrier frequency is different from what the receiver thinks it should be.

You can write the passband signal as

$$u_p(t)=\text{Re}\left(I_{[-1,1]}e^{j400\pi t}\right)\tag{1}$$

Since you're demodulating with a frequency of $\omega_0=401\pi$, it is advantageous to rewrite the passband signal with that reference frequency:

$$u_p(t)=\text{Re}\left(I_{[-1,1]}e^{-j\pi t}e^{j401\pi t}\right)\tag{2}$$

From $(2)$ you can read off the new complex envelop with respect to the actual demodulation frequency $\omega_0=401\pi$:

$$u(t)=u_c(t)+ju_s(t)=I_{[-1,1]}e^{-j\pi t}=I_{[-1,1]}\cos(\pi t)-jI_{[-1,1]}\sin(\pi t)$$

Demodulation with $\cos(401\pi t)$ will give you the I component $u_c(t)$. The factor $1/2$ simply occurs because by multiplying with a cosine you get $u_c(t)\cos^2(\omega_0t)=u_c(t)(\frac12+\frac12\cos(2\omega_0t))$, and the LPF removes the component at $2\omega_0$, which leaves you with $\frac12 u_c(t)$.


Finding when a part of audio starts and ends in another audio I have two audio files in which a sentence is read (like singing a song) by two different people. So they have different lengths. They are just vocal, no instrument in it.

A1: Audio File 1
A2: Audio File 2
Sample sentence :  "Lorem ipsum dolor sit amet, ..."

I know the time (in decimal seconds) every word starts and ends in A1. And I need to find automatically that what time every word starts and ends in A2. So, how to find sound of a word (in A1) in another audio (A2) that has different duration (of word) and different voice?

You may consider using Dynamic Time Warping (DTW) for this.
One approach which generally would work is as follows:

Divide both the audios into smaller (analysis) windows with appropriate overlap as per the needed temporal resolution
Transform the audio from time domain to cepstral. You may compute  MFCCs.
Copute DTW to analyse the warping


How to filter OFDM Signal with an analog filter I am using OFDM system in Matlab, I want to use an analog filter after the DAC. Because we don't a simulation for the DAC, the generated OFDM analog signal does not have a specific frequency (they only have normalized frequency). When I design the analog filter, I have to specify the frequency (not normalized frequency as in digital filters). So I am not able to filter my OFDM signal with the analog filter. Any ideas on how to do that?
I have been trying to use slim Matlab command to do that, I managed to apply an analog filter on sinusoidal signals as follows:
clc
clear all
close all
f1 = 1000;                  % signal frequency
f2 = 2000;                  % signal frequency
f3 = 3000;                  % signal frequency
fs= 10000;                  % DAC sampling Rate
Ts= 1/fs;                   % DAC sampling time
N= 4000;                    % Number of points
n= 0:N-1;                   % sampling index
x= cos(2*pi*f1*n*Ts)+ cos(2*pi*f2*n*Ts)+cos(2*pi*f3*n*Ts);
df=-0.5:1/N:0.5-(1/N);      %frequency bins
X=fftshift(fft(x,N));
figure;plot(fs*df,abs(X)/N);    %original signal

%%
%filter design
order=14    % Filter order
Rs=30   %stopband attenuation
fc=1500; % filter cutoff frequency
[b,a] = cheby2(order,Rs,[2*pi*fc],'low','s');
figure;freqs(b,a);
% [z2,p2,k2] = cheby2(order,Rs,2*pi*fc,'low','s');
% [b2,a2] = zp2tf(z2,p2,k2);
% figure;freqs(b2,a2);
% [zeros, poles, gain] = butter(order,[2*pi*fc],'low','s'); %generate LPF zeros, poles, and filter gain.
% [b2,a2] = zp2tf(zeros, poles, gain); %convert to transfer function form of b, a coefficients
% figure;freqs(b2,a2);

%%
%apply the filter
% create system
myFilter = tf(b,a);
% apply filter to time domain signal
[y, time] = lsim(myFilter,x,n*Ts);
Y=fftshift(fft(y));
figure;plot(fs*df,abs(Y)/N);

same as your previous RRC question: you simply apply a digital filter that approximates the analog filter.
There's no analog signals in your PC, and OFDM is a digital signal anyway, so only the digital equivalent of the analog system would matter to the receiver, anyway.
Look at it this way: You're supposed to simulate a digital system. Therefore, the receiver is digital. Therefore, the analog effects that would happen in a real-world transceiver system are only relevant in their digital form.
And if Nyquist has taught you anything, any band-limited analog system can be represented through a digital one.
So, stop worrying about building analog filters, but describe the thing you want to simulate and then implement that digitally.
(I mean, really! Look at the purpose of OFDM: it's a purely digital thing, meant to do the inverse to a purely analog multipath channel, which is an analog filter. It works. You don't have to model anything in analog, unless you're not really understanding what OFDM is about.)


3 dB bandwidth of a filter I am solving one of the question from signal processing about calculating the 3dB bandwidth of the filter. But I am stuck once I got the magnitude response of the filter. Can someone help me with this?
Thanks

From the given transfer function you can see that the filter has a low pass character. The gain at DC is given by $H(1)=1$. This is also the maximum value of $|H(e^{j\omega})|$. So you can try to solve
$$\big|H(e^{j\omega_c})\big|^2=\frac{\max_{\omega}|H(e^{j\omega})\big|}{2}=\frac12\tag{1}$$
because this equation defines the frequency $\omega_c$ at which the filter gain equals $1/\sqrt{2}$, which is the definition of the $3$ dB cut-off frequency.
Your calculation of the magnitude of $H(e^{j\omega})$ is correct, so you can just continue by solving $(1)$.
As a final note, there are values of $a$ for which the $3$ dB cut-off frequency doesn't exist because the filter never reaches an attenuation of $3$ dB. You can try to find that limit of $a$ as a bonus exercise.


Outlier Detection after Detrending a Time Series With Missing Values or NaN Goal

Substitute outliers in a time series by most recent valid data

Problem

The time series (end-of-day stock prices) has several 'uncomfortable' properties:

It is non-stationary and can have components of low and high frequency (trends and sudden price moves)
There can be missing data for some days, typically single occurrences, but sometimes a whole consecutive segment is missing
Some values are erroneous (implausibly high or low i.e. due to misplaced decimal sign). These are the real outliers and don't belong to the series.

There is an additional caveat: The substitution algorithm must not look into the future of the time series since it will be used in a machine learning scenario.

Strategy

What I thought of so far:

Detrending the series. Maybe rolling linear detrending by the last n values. Will still be a problem with very short term moves.
Replacing the missing values by extrapolating the trend identified in (1)
Calculating the z-score and replacing values with z > 6 by extrapolating the trend identified in (1)

As the tags suggest, I will realise this in python / scipy. However I don't have much experience with time series analysis. I have no idea if this is an appropriate approach or if I'm missing something important. Any help is appreciated. Thanks.

Problem

  The time series (end-of-day stock prices) has several 'uncomfortable' properties:

  It is non-stationary and can have components of low and high frequency (trends and sudden price moves)

These are not outliers though. These are valid signal behaviour.

  There can be missing data for some days, typically single occurrences, but sometimes a whole consecutive segment is missing
  Some values are erroneous (implausibly high or low i.e. due to misplaced decimal sign). These are the real outliers and don't belong to the series.

To reduce the effect of implausibly high or low values you can use a rolling median filter at a suitable time scale (e.g. depending on your sampling frequency, you might want to get the median of a day's transactions).

Missing data is a more challenging problem. A few missing values might be "patched" with some interpolation scheme but if you have almost a day's worth transactions missing, it is probably better to break the analysis into segments too. Otherwise, you would be assuming a lot for a part of the signal that you don't really know what happened in reality.

Hope this helps.


How can I distinguish between two similar EEG signals? I have two EEG signals that are very similar. The difference is only in amplitudes. However, they are coming from two different cognitive processes. What are some methods, beside FFT, for distinguishing between two very similar signals in terms of statistical information? Please let me know. Thank you

BSS.

yet another thing to try could be wavelet transform. Wavelets are powerful tools for analysis of signal. You could take a DWT of your signal and compare may be the most significant coefficients. Anyway but it would be helpful if you could post some sample data.


How to detect the presence one- and two-dimensional patterns in distributions of spatial points? I have many point patterns, and I want to detect those patterns that have points arranged in lines that may be curved, intersected or broken. For example, on the following picture, the upper left pattern is not interesting but other three patterns are those that I want to detect.

Not all points in the interesting patterns are aligned. Some points in those patterns are still random. The lines do not have any parametric form, so I cannot use the Hough transform. I searched some literature on spatial statistics, for example this one and only found metrics for assessing
spatial uniformity and detecting clusters such as the index of dispersion. Which metrics exist for detecting the presence of not only clusters but curved lines?
I also need to do the same analysis for three-dimensional points. How to detect those 3D patterns that have points roughly arranged in lines and surfaces that can be wiggly?

Perhaps you could perform blob detection to find the individual dots and then map their location to a 2d space, that preserves their spatial structure, through the use of some measure for example the center of mass of the blob.
Once you have mapped each blob to its coordinates in the 2d space you could perform clustering or fit a regressor since the patterns somewhat resemble piecewise polynomials, you could model them as second or third degree polynomials, and build your algorithm upon RANSAC, it is a probabilistic algorithm that fits a model while discarding the outliers, you could also introduce some extra hyperparameters:

the percent threshold of points that minimize the cost function (you could use the sum of squared errors as a base cost function to build upon) and thus belong in the pattern.
the number of piecewise relations if you opt in the piecewise polynomial approach. Then you could exclude valid points that fit one branch of the piece wise polynomials from next iterations until all such polynomials are found.

Another classic algorithm that can cluster points based on their geodesic structure rather than just the distance would be DBSCAN, choosing a suitable $\epsilon$ for your dataset as the points that belong to the pattern are rather close to each other. A more complex classical machine learning algorithm is ISOMAP, it can also cluster points that based on geodesic properties


How to get some law parameters from its FFT figure? Given an hidden law $$I(t)=\frac{1}{2}[A+B.cos(2\pi\nu_0t+\phi_0)]$$ I observe its signal with 256 samples and compute its FFT over this interval, and  represent the FFT result from 0 to 1 (normalised frequency).

I get a figure whose pic is at 0.83.

How can I retreive A,B and $\phi_0$ from this graphics ?

A second question is : is it normal that this figure is not symmetric ?

A final question is : doing an analytical computation on classical Fourier transform gives a result. How much does what I observe from FFT and the analytical result differ ?

If you feed an FFT a finite set of samples from a function of theoretically infinite range, you have implicitly windowed your function with a rectangular window before the FFT.  Thus, the FFT result will be the convolution of the spectrum of your function and the transform of a rectangular window (a periodic Sinc).  This convolution will be an identity operation only if the samples happen to span exactly an integer multiple of periods of a purely periodic function.  So how much your FFT result will change or differ will depend on the offset and length of your FFT aperture/window with respect to some reference point on and the scale of your function.

In the case of a single sinusoid in zero noise, I reference 2 papers by Clay Turner on solving for the 3 unknowns using only 3 or 4 non-aliased samples, somewhere on my web page: http://www.nicholson.com/rhn/dsp.html   No FFT needed.


Upsampling signal using convolution-based interpolation filters I am currently reading this paper which discusses several image interpolation methods, such as nearest neighbor and linear interpolation, using convolution filters. I first want to do this in 1D with discrete-time signals (e.g. [0 0 1 0 0]), however, I have some trouble understanding on how to derive the convolution filters from the mathematical definitions:

$$h_{NN}(x)=\cases{
1, &if$\quad-\frac{1}{2} \leqslant x < \frac{1}{2},$\\
0, &otherwise,
}$$
$$h_{Lin}(x)=\cases{
1-|x|, &if$\quad 0 \leqslant |x| < 1,$\\
0, &if$\quad 1 \leqslant |x|.$
}$$

My main question is, how do I calculate the filter weights as well as the filter size? Furthermore, I am aware of the fact that I need to convolve the signal multiple times to upsample it, which the following equation in the paper also shows.

$$I(\mathbf{x}) = (\dots((I_s(\mathbf{p})*h(x_1))*h(x_2))*\dots)*h(x_N),$$

What are the several x in there? In case of images, they are the positions of the original image (before it is reconstructed from $I_s(p)$), but how do I apply this concept to one-dimensional discrete-time signals?

I would be thankful for any help!

Starting from the equation with the convolutions and assuming unit distance between sampling grid points, for a one-dimensional signal you'll have just:

$$I(\mathbf{x}) = I((x_1))= I_s((p_1))*h(x_1),$$

where $\mathbf{x}$ is the vector of continuous coordinates of which there is only one: $x_1 \in \mathbb{R}.$ Likewise there is just a single discrete coordinate $p_1 \in \mathbb{Z}.$ For simplicity, you could drop the vector notation and subscripts altogether and have just:

$$I(x) = I_s(p)*h(x).$$

Without having fully read the paper, I don't think the authors dare to venture into defining what is meant by convolution between discrete-time $I_s(p)$ and continuous-time $h(x)$ which results in continuous-time $I(x)$. Such a convolution can be defined as:

$$I(x) = \sum_{p=-\infty}^{\infty}I_s(p)h(x-p)$$

For example, if you want to obtain $I(x)$ within the inter-sample interval $0 \le x < 1$ and you use the linear interpolation kernel, then the above simplifies as:

$$\begin{eqnarray}I(x) &=& \sum_{p=-\infty}^{\infty}I_s(p)h_{Lin}(x-p)\\
&=& \sum_{p=-\infty}^{\infty}I_s(p)\times\cases{
1-|x-p|, &if$\quad 0 \leqslant |x-p| < 1,$\\
0, &if$\quad 1 \leqslant |x-p|$
}\\
&=& \sum_{p=0}^{1}I_s(p)(1-|x-p|)\\
&=& I_s(0)(1-x) + I_s(1)(1-(1-x))\\
&=& I_s(0)(1-x) + I_s(1)x.\end{eqnarray}$$

Or likewise for any inter-sample interval $x_0 \le x < x_0 + 1:$

$$I(x) = I_s(x_0)(1-(x - x_0)) + I_s(x_0 + 1)(x - x_0).$$

Another way to understand such convolutions is that $I_s(p)$ is the product of the continuous-time signal multiplied by a Dirac comb, and to use a continuous-time convolution integral.

As you can see from the above, the interpolation filter is not a regular discrete-time filter, but one that converts a discrete-time signal into a continuous-time signal. The convolution filter definitions cited in the question are directly applicable, without further "sizing", when the sampling interval is 1.


Is it possible to apply a frequency filter a time domain signal? Is it possible to frequency filter a time domain signal, without converting the signal to a frequency domain..

My question comes in regards how MFCC according to wiki and other pages describe the extracting step.

Is this necessary though?

When you want to explicitly apply a frequency domain operation on a given signal $x[n]$ then you need to obtain those frequency domain samples $X[k]$ of $x[n]$,

However, it's possible that you can also apply an equivalent operation on the time domain samples $x[n]$ which will implicitly realize the proposed frequency domain operation on $X[k]$ as a result. However in order to accomplish that, you shuld be able to convert the frequency domain operation into a time domain equivalent, usually by applying Fourier transform properties.

If you cannot find such a (simple) time domain equivalent of frequency domain operation, then, no, you cannot apply it on time domain samples.

As an example if you would multiply the DTFT $X(e^{j\omega})$ by a rectengular window of cotoff frequency $\omega_c$; $Y(e^{j \omega}) = W(e^{j \omega}) X(e^{j \omega})$ then you could also perform it on the time domain as a convolution: $y[n] = w[n] \star x[n]$ , but if you would do the following operation on the frequency samples $Y(e^{j \omega}) = log( |X(e^{j \omega})| )$ then it's hard to find an simple time domain equivalent case...

(Note that I've used the theoretical DTFT, $X(e^{j\omega})$ to show the expressions while referring to the frequency samples which actually belong to DFT $X[k]$, please keep in mind the circularity issue when finding the exact expression.)

For the mel-scale operations that are included in the box such as the $\log$ operation etc, the degree of nonlinearity will probably disallow an (simple) equivalent time domain expression to be formulated, and I don't know one such, but if you can find one, then why not.


Difference Between Linear Convolution and Circular Convolution If I understood correctly (and this page should confirm: http://www.cs.ioc.ee/~khoros2/linear/convolution-teo/front-page.html)

if I convolve linearly (the usual point-to-point multiplication and addition without tricks for each element of the image matrix) an image with a kernel I obtain RESULT1
if I calculate the DFT of the image and of the padded kernel (to have same size as the image) and perform their multiplication in frequency and then the inverse DFT, the result is RESULT2 and should be like applying the circular convolution in time space

RESULT1 slightly differs from RESULT2 as in the example

Is this correct?

The difference applies only to the borders of the image.

In the linear convolution you assume the values of pixels beyond the border (examples being mirror of the image pixels, or 50% grey).

In the circular convolution (or DFT, product, IDFT), the pixels beyond the border are the pixels on the other side of the image, just as if you had a repeated tiling of the image.


Z-domain transfer function question For the difference equation below:

$y(k)=1/2{x(k)+x(k-1)}$

How i can find he z-domain transfer function?

Thanks.

I assume this is homework, but transforming a difference equation to the $z$-domain is simple; just recall the time-shifting property of the transform.

$$
x[n] \Leftrightarrow X(z) \rightarrow x[n-k] \Leftrightarrow z^{-k}X(z)
$$

So then we have:

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

$$
Y(z) = \frac{1}{2}X(z) + z^{-1}X(z)
$$

The transfer function can be written as:

$$
H(z) = \frac{Y(z)}{X(z)} = \frac{1}{2} + z^{-1}
$$


How can I calculate the cyclic (periodic) convolution? I'd like to understand how to calculate the cyclic convolution as well as understand what that means exactly. How should I go about finding the output for various periods for a system?

I have an example:

$$
x(n) =
\begin{cases}
n & \textrm{for} \quad 1 \leq n \leq 3\\ 0 &\textrm{otherwise}\end{cases}
$$

$$
h(n) =
\begin{cases}
n & \textrm{for} \quad 1 \leq n \leq 2\\ 0 & \textrm{otherwise}
\end{cases}
$$

If I perform the convolution, then I get the following values for $y(n)$:

$$y(2) = 1;\quad y(3) = 4;\quad y(4) = 7;\quad y(5) = 6$$

Now, if I want a period = 3, then:

$$x(n) = x(n+3k) \quad\textrm{and}\quad h(n) = h(n+3k)$$

At this point, I'm unsure of what to do to get values that correspond to a period.

People generally define the cyclic convolution of periodic sequences
$x$ and $h$ of period $N$ as
$$y[n] = \sum_{m=0}^{N-1}x[m]h[n-m], n = 0, 1, \ldots, N-1.\tag{1}$$
Note that the above expression consists of $N$ different sums that
you have to compute, and if while computing any particular sum,
the value of $n-m$ is not in the range of numbers for which you
are given $h[\cdot]$, then you use the periodicity ($h[n-m] = h[n-m+N]$
or $h[n-m] = h[n-m-N]$) to get the argument into the range for which
you know the value of $h[\cdot]$.
Also, note that $(1)$ holds for all integers $n$, but we don't
need to calculate more than $N$ sums like $(1)$ because
$y[n]$ is also a periodic sequence of period $N$ and so we have
for any integer
$M$ that  $y[M] = y[M \bmod N]$ where, of course, $0 \leq M \bmod N \leq N-1$.

Exercise: write out the above formula explicitly, meaning
no summations, for $n = 0, 1, 2$ and proceed from there.  Go on;
you can do it. There are only three sums of three terms each.


Steady-State Output from Transfer Function The progress I have made is as follows:

$\sin(t)$ is our signal therefore $\omega = 1 = 2\pi f$ and $f$ = $\frac{1}{2 \pi}$

Also, $f_s$ = 10Hz therefore T = $\frac{1}{f_s}$ = 0.1s

$H(z) = \frac{z}{z-0.7}$ so let $z = e^{j\omega T}$ where $\omega T = 1 \times0.1 = 0.1$

$H(\omega) = \frac{e^{j0.1}}{e^{j0.1} - 0.7}$

Taking the denominator:

$\cos(0.1) + j\sin(0.1) = 0.995 + 0.099j$

therefore,

$\cos(0.1) -0.7 + j\sin(0.1) = 0.311∠0.324$

$H(\omega) = \frac{1∠0.1}{0.311∠0.324}$

$H(\omega) = 3.21∠0.224$

From here I am out of ideas on how to continue.  Any advice appreciated.

You may know that one important property of linear time-invariant (LTI) systems is that the complex exponential $e^{j\omega_0}$ is an eigenfunction, and the corresponding eigenvalue is given by the system's frequency response evaluated at $\omega_0$.

So the response to an input signal $$x[n]=e^{jn\omega_0}$$ is given by $$y[n]=H(e^{j\omega_0})e^{jn\omega_0}$$

Now try to show that from this it follows that the response to a sinusoidal input signal $$x[n]=\sin(\omega_0n)$$ is given by $$y[n]=\big|H(e^{j\omega_0})\big|\sin\left[\omega_0n+\angle H(e^{j\omega_0})\right]$$ (if the system is real-valued). Now should you have all the information to compute the reconstructed output signal.


Doppler Shifting A Chirp Signal So I have a digital linear chirp signal:

$$ x[n] = \sqrt{\frac{2}{N}} \cdot \cos\Big( 2 \pi \frac{f_{0}}{f_{s}}n + \pi \frac{B}{N} n^{2}\Big) $$

Where $f_{0}$ is the initial frequency, $f_{s}$ is the sampling frequency, $B$ is the bandwidth and $N$ is the number of samples.

I want to simulate a Doppler shift. I am not sure exactly how to do this.

I know I need to multiply by initial frequency $f_{0}$ by the Doppler shift $d$, and do the same with the bandwidth. This will give me:

$$ x_{d}[n] = \sqrt{\frac{2}{N}} \cdot \cos\Big( 2 \pi \frac{d \cdot f_{0}}{f_{s}}n + \pi \frac{d \cdot B}{N} n^{2}\Big) $$

I know Doppler compresses or expands my waveform in time, and I know how to simulate the Doppler shift in the case of a simple sinusoid. However, I am a bit confused here.

Do I need to change the number of samples $N$ as well? If so, by how much exactly?

Thanks

You need to resample the signal to simulate the Doppler induced dilation. The resampling factor is I=Td/Ts
where
Td = Duration of the signal after dilation
Ts = Actual duration of the transmitted signal

In Matlab you can use the function resample(), but you need to find the resampling parameters P and Q from I.

Also, I=1+v/c for signal expansion and I=1-v/c for signal compression. v=speed of movement in m/s and c=speed of wave propagation in the medium. For instance sound underwater c=1500 m/s.

In Matlab, you can find P and Q as [P,Q]=rat(I);


FFT to spectrum in decibel Here is a 10 seconds-long 440hz sine wave
 normalized at $0\textrm{ dBFS}$.

When computing the STFT (with the code below) of this audio file, I noticed that max(abs(STFT)) is around 248.33.  (more generally, it seems to be approximately fftsize/4 for this particular file).

Even after doing magdB = 20 * math.log10(abs(STFT)) to have decibel values, we get a max of 47.9 dB.

But this $47.9\textrm{ dB}$ doesn't mean anything! We would like to have ~ $0\textrm{ dB}$ instead because this audio only contains 1 sine wave component (1 "harmonic") at $0\textrm{ dB}$ volume.

Audio tools display it properly, here is the spectrum display of this sine wave:

Question: how to have an absolute, canonical $\rm dB$ values from a FFT, that makes that a pure $0\textrm{ dBFS}$ sinewave has a peak of $0\textrm{ dB}$ in the spectrum display ?
More generally, is there a canonical way to go from FFT values to $\textrm{dB}$ in order to display a spectrum analysis?

NOTE: I know that $\textrm{dB}$ are for ratio between 2 things, so all I mentioned here is referred to $\textrm{dBFS}$ (full-scale)

import scipy, numpy as np
import scipy.io.wavfile as wavfile

def stft(x, fftsize=1024, overlap=4):
    hop = fftsize / overlap
    w = scipy.hanning(fftsize)
    return np.array([np.fft.rfft(w*x[i:i+fftsize]) for i in range(0, len(x)-fftsize, hop)])

Definitely you will have to calibrate your system. You need to know what is the relationship between dBFS (Decibel Full-Scale) and dB scale you want to measure.

In case of digital microphones, you will find sensitivity given in dBFS. This corresponds to dBFS level, given 94 dB SPL (Sound Pressure Level). For example this microphone for input $94 \;\mathrm{dB SPL}$ will produce signal at $-26 \;\mathrm{dBFS}$. Therefore the calibration factor for spectrum will be equal to $94+26=120\;\mathrm{dB}$.

Secondly, keep in mind scaling of the spectrum while doing windowing and DFT itself. More specifically, given amplitude spectrum (abs(sp)):

Total energy of the signal is spread over frequencies below and above Nyquist frequency. Naturally we are interested only in half of the spectrum. That is why, it important to multiply it by 2 and ignore everything above the Nyquist frequency.
Whenever doing windowing, it is necessary to compensate for loss of energy due to multiplication by that window. This is defined as division by sum of window samples (sum(win)). In case of rectangular window (or now window), it is as simple as division by N, where N is the DFT length.

Here is some example source code in Python. I am sure that you can take it from here.

#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as wf

plt.close('all')

def dbfft(x, fs, win=None, ref=32768):
    """
    Calculate spectrum in dB scale
    Args:
        x: input signal
        fs: sampling frequency
        win: vector containing window samples (same length as x).
             If not provided, then rectangular window is used by default.
        ref: reference value used for dBFS scale. 32768 for int16 and 1 for float

    Returns:
        freq: frequency vector
        s_db: spectrum in dB scale
    """

    N = len(x)  # Length of input sequence

    if win is None:
        win = np.ones(1, N)
    if len(x) != len(win):
            raise ValueError('Signal and window must be of the same length')
    x = x * win

    # Calculate real FFT and frequency vector
    sp = np.fft.rfft(x)
    freq = np.arange((N / 2) + 1) / (float(N) / fs)

    # Scale the magnitude of FFT by window and factor of 2,
    # because we are using half of FFT spectrum.
    s_mag = np.abs(sp) * 2 / np.sum(win)

    # Convert to dBFS
    s_dbfs = 20 * np.log10(s_mag/ref)

    return freq, s_dbfs

def main():
    # Load the file
    fs, signal = wf.read('Sine_440hz_0dB_10seconds_44.1khz_16bit_mono.wav')

    # Take slice
    N = 8192
    win = np.hamming(N)
    freq, s_dbfs = dbfft(signal[0:N], fs, win)

    # Scale from dBFS to dB
    K = 120
    s_db = s_dbfs + K

    plt.plot(freq, s_db)
    plt.grid(True)
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude [dB]')
    plt.show()

if __name__ == "__main__":
    main()

Lastly, there is no better source on this type of spectrum scaling than brilliant publication by G. Heinzel et al. Just keep in mind that if you want to proper RMS scaling, then full-scale sinosuidal signal will be $-3 \;\mathrm{dBFS}$.

You might also find my previous answer being helpful.


General questions regarding OFDM modulation In Orthogonal frequency division multiplexing (OFDM), data symbols are modulated using an N pt-IFFT operation to produce ONE output OFDM

Assume we have total of $N$ QAM symbols which are input to the IDFT operation.

1) Is the output time domain samples $X_k$ also of length $N$ in length?

2) If ONE data QAM symbol has symbol duration $T$, is ONE OFDM symbol of symbol duration $T$ or is one SAMPLE of the OFDM symbol of duration $T$?

Thanks looking forward for your answers

Update

One example of OFDM modulator is the one used in LTE, in particular the subcarrier bandwidth is 15 KHz, and the OFDM symbol (and not one sample of the OFDM symbol) is 1/15k = 66.67 $\mu $ sec which is equal to the transmission time of one symbol (for example 1 QAM) over one subarrier. If my understanding is correct than the answer provided below is wrong. Can somoeone please double check?

Thanks

The output time domain samples are also of same length $N$. Your figure shows a part of a OFDM receiver. At the OFDM transmitter, IFFT is used.
One OFDM symbol duration is $NT$ and obviously one sample of the OFDM symbol has a duration $T$.


MATLAB's $\tt dsp.ParametricEQFilter$ bandwidth question Consider the following code:

Fs = 48000;
h = dsp.ParametricEQFilter('CenterFrequency', 1000, ...
                           'Bandwidth', 1400, ...
                           'PeakGaindB', 10, ...
                           'SampleRate', Fs);

fvtool(h,'Fs',Fs,'FrequencyScale','Log');

Let the two bandwidth $-3\rm dB$ points $f_1$ and $f_2$.
Therefore: $f_1 = 300\ \mathrm{Hz}, \ f_2 = 1700\ \textrm{Hz}$

However, in the graph I see the $-3\ \mathrm{dB}$ point (i.e. the $+7\ \rm dB$ gain) at $487\ \rm Hz$. What am I doing wrong?

Apparently, after some research, I have found that the bandwidth "bandedges" is a bit of a tricky subject: (see Bandwidth Re-visited).
Some references suggest indeed that the bandwidth is measured at the gain which is equal to half of the peak gain. Some others say that you should measure the bandwidth at $ +3\ \textrm{dB}$ for boost gain, $ -3\ \textrm{dB}$ for cut. (which eventually runs into the issue of where it should be measured if the gain is less than $ 3\ \textrm{dB}$.

Hilmar is right about the geometric mean of the corner frequencies. However, the cutoff points (i.e. the bandedges) do not exist at half the peak gain in MATLAB.   They exist at $G^{2}_B = \frac{1+G^2}{2}$, where $G$ is the parametric equalizer linear gain, and $G_B$ is the linear bandwidth gain as described here.


Bluetooth Codec (SBC) - Uses a FFT? In the A2DP specification, the SBC code uses a "polyphase filterbank" which has the formula
$$x[n] = X[m]\cos((m+0.5)(n+M/2)\pi/M).$$
I am confused to what this is? It seems to take in $N$ points and output $2N$ points in some sort of frequency domain. Is this a special type of DFT and can it be optimized with a radix-2 algorithm? I tried searching but could not find any websites that could help me.

Edit: Links - https://www.bluetooth.org/docman/handlers/downloaddoc.ashx?doc_id=8236 and bottom of Page 64

I figured out the answer. The polyphase filterbank in the SBC codec is a variation on the discrete cosine transform type 4 that also divides it into frequency bands at the same time. If you look in the spec, the double for loop is used to both to multiple across all basis and sum up the frequency bands.


Calculating PSD, Cross PSD by using FFT outputs I am having the FFT results of two distinct voice signals. Let's say the FFT results of these two signals are of the form (a+ib) and (c+id). Now, I want to calculate the PSD of each signal and also the Cross PSD between them. I'm calculating PSD and Cross PSD as follows.

PSD of signal 1 = a * a + b * b;

PSD of signal 2 = c * c + d * d;

Cross PSD between signal1 and signal2 = (ac - bd) + (ad + bc). I also want to know whether this equation represents cross-correlation or cross PSD?

And also, I have ambiguity about the importance of cross PSD. I know that cross PSD is obtained by taking the Fourier transform of the cross-correlation results. I know that from PSD, we can get an idea of the power for different frequencies present in the signal. How does cross PSD differ from PSD? I found out that by using cross PSD, we can find the power shared by a given frequency for the two signals. How do we know how much power is shared between the two signals from a single constant power value per frequency? Please help me out in understanding what the cross PSD means and what is its importance?

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

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


Nyquist noiseless channel capacity; how can bit-rate be two times the bandwidth? I'm confused by the Nyquist channel capacity formula.  How can channel maximum capacity approach double the bandwidth.

$C = 2\times BW \times log_{2}(L)$ bits/sec

The way it was explained to me way back when when I took a telecommunications course was as follows.

C = Maximum channel capacity in bits/second

BW = (highest frequency - lowest frequency) that the channel can accommodate.

L = signal levels So in the figure below, where BW = f; the frequency of the sine wave; and L = 2 for -1 and + 1; the channel capacity is C
= 2B because we can consider the positive part of the sine wave as = to a 1 bit for half the period, and as a 0 bit for the negative half of the period.

Now i'm not sure that this is correct, as that doesn't seem like "information", being that it doesn't seem to decrease entropy, or uncertainty?

so if that's not correct, what is an example of a signal with capacity twice the bandwidth?
Or am I just not understanding the entire concept?

I think you're confusing two different (but related) terms.

Nyquist says that in a channel of bandwidth $B$ you can transmit up to $2B$ orthogonal pulses per second. So, $R_p \leq 2B$, where $R_p$ is the pulse rate.

To achieve $R_p = 2B$, the pulses need to be sinc-shaped. Other, more practical pulses achieve slightly less than that. For example, raised cosine pulses with 50% excess bandwidth achieve $R_p = 1.5B$.

If you send pulses with $L=2^k$ amplitude levels, then the bit rate is $R_b = kR_p$.

So far for Nyquist, who was only interested in maintaining orthogonality. Shannon wanted to find a limit for the bit rate of reliable communication, with probability of error $P_b \rightarrow 0$. This limit is the capacity $C = B \log_2(1+\text{SNR})$, assuming an AWGN channel and orthogonal pulses. Note that $L$ is not mentioned explicitly.

  what is an example of a signal with capacity twice the bandwidth?

The capacity is a property of the channel, not the signal. A signal has a bandwidth, and implicitly a pulse rate and a bit rate. Regardless of that, the channel imposes a limit on the bit rate (its capacity), if you want reliable communication.


Single channel speech enhancement What is the difference between Magnitude spectral subtraction and Power spectral subtraction in terms of performance?

If by "performance" you mean the quality of the processed (speech) signal, then this question cannot be answered. The reason is that the quality of the output of a spectral subtraction algorithm depends on so many factors, and the precise subtraction rule is just one detail. What is much more important is the noise estimation algorithm, certain smoothing operations, etc.

What you can say is that the dynamic range of algorithms working with magnitudes is much smaller than when powers are used. This can be a great advantage when such systems need to be implemented in fixed-point arithmetic.


Design single pass band FIR filter I'm new in DSP. I want to design single pass band FIR filter using LS method. I have this information :

one pass band $[-0.1\pi , 0.3\pi]$
two stop bands $[-\pi , -0.18\pi]$, $[0.38\pi, \pi]$
weighting constant in :

pass band : $a = 1$ and stop band $b_1=b_2=2$.

In this case, can I use Matlab's firls function? Because we have negative frequency and f is a vector of pairs of frequency points, specified in the range between $0$ and $1$.

As mentioned in a comment, you can't use firls directly because this function only designs real-valued filters with (conjugate) symmetrical frequency responses. Your specification can only be met by a complex-valued filter because the frequency domain specification is not symmetrical.

In principle, you could use firls to design the real part and the imaginary part of the impulse response separately (by using the even and odd parts, respectively, of the frequency domain specification as the desired responses), but in that case you would get discontinuous desired responses, which would result in large approximation errors.

A better and simpler approach is to set up an overdetermined system of complex linear equations and use Matlab/Octave to solve them in a least squares sense. This is very straightforward if you use matrix commands:

N = 51;     % desired filter length
% frequency grid, desired frequency response, weighting function
f = [linspace(-1,-.18,164),linspace(-.1,.3,80),linspace(.38,1,124)];
d = [zeros(1,164),ones(1,80),zeros(1,124)].*exp(-1i*pi*f*(N-1)/2);
w = [2*ones(1,164),ones(1,80),2*ones(1,124)];
f = f(:); d = d(:); w = w(:);

% set up and solve overdetermined linear system
A = w(:,ones(1,N)) .* exp(-1i*pi*f*(0:N-1) );
h = A \ (w.*d);

The resulting filter is complex-valued and has a linear phase response, i.e., the real part of the impulse response is even, and the imaginary part is odd (see figure):

You can further decrease the approximation error by chosing a larger filter length.

I also have a corresponding Matlab/Octave function on GitHub: cfirls.m.


Compute PSD from a CSV file I'm sorry if my question might be simple, but I don't have a strong background on signal processing. I have a csv file of a 8 MHz bandwidth, sampled with a resolution bandwidth of 100 KHz. What I want to do is to compute the PSD over the 8 MHz bandwidth, but I'm not sure if I am doing this correctly.

For a "stable" signal, I multiply the received power over the 100 KHz sample, and then average it at the end over 8 MHz. Is that correct? Does this give me a rough estimation of the PSD?

Are there tools like octave or even python programming that can help in this process? I looked on the internet but didn't had luck.

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

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

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

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


Analog signal can be discrete time? I am new to this site and field too.I always thought analog is continuous and digital signal is discrete. I read this today and got confused about analog signal

  Analog can be
  continuous time(CT) or discrete time(DT). A discrete time signal is not a digital signal; a
  discrete time signal is one in which the amplitude is a continuum, and not discretized. On
  the other hand, if the time is discretized, this signal is still analog. Only when a discrete
  time signal is passed through an A to D converter, it becomes a digital signal.

Will anybody explain me what actually it means?

This sounds like a confusion in terminology. See http://en.wikipedia.org/wiki/Digitizing

Digitization basically involves two steps:

Discretization: Sampling the signal at discrete times
Quantization: Turning the samples from (in theory) infinite resolution to finite resolution

Look at all possible combinations you can have 4 different types of signal:

continuous time, continous amplitude
discrete time, continous amplitude
continuous time, discrete amplitude
discrete time, discrete amplitude

Typically we refer to #1 as "analog" and #4 as "digital". #2 and #3 are of little practical relevance, other than getting the mathematical models right. Quantization and Discretization are really quite different mathematically so it doesn't make a lot of sense to analyze them in combination. In other words: the theory of time discretization applies equally to continuous amplitude and discrete amplitude signals.


Real-Time Human Pitch Detection I'm trying to implement a singing game that will analise raw mic input and tell the player how good is he singing. That needs to be done in real-time.

I've come across a lot of threads asking the same question but I'm still quite not done with it, probably due to my lack of experience in the field and shallow math background. I've implemented an algorithm based on the DSPDimension website pitch shift's article: http://www.dspdimension.com/admin/pitch-shifting-using-the-ft/

I extract the true frequency and magnitude just like the article explains, but I don't know find the fundamental frequency with this. I've tried to get the bin with greatest magnitude but that only give me right results for higher pitch signals, it doesn't matter which oversampling factor I use I still get bad data for low freq signals. Is this approach completely wrong or am I in the right track but just missing something?

Thanks in Advance,

EDIT: I forgot to mention that I'm only interested in the pitch class, so it is ok if the fundamental is missing but I have a strong overtone in the sample.

EDIT2: Thanks to everyone, I just finished a version of the algorithm that's working like a charm. The low pitch estimation problem was due to my input test. When I sung the note it matched correctly. Also, I'm considering all harmonics now, not just the highest peak.

I've tried to get the bin with greatest magnitude but that only give me right results for higher pitch signals, it doesn't matter which oversampling factor I use I still get bad data for low freq signals.

That's because the harmonics are larger than the fundamental.  Plot your spectrum and you'll see.  A better method to find the true fundamental is autocorrelation.  Then you're "sliding" the waveform past itself and finding delays at which the wave shape lines up with itself.

http://www.phon.ucl.ac.uk/courses/spsci/matlab/lect10.html

Do you really want them to sing the exact note, or is it ok if they sing an octave above or below depending on their voice register?


Spectral images - How do they appear during ADC? Why is it that spectral images appear? I know they appear for $F_s$, but why those replicas appear, I don't know.

E.g., sampling $x(t) = sin(1000 * 2\pi t)$ with a sample frequency of 230 Hz.
I know it will alias, due to spectral replicas, but why do those appear at all?

One way to represent AD conversion mathematically is to think of it as a multiplication of the signal by an impulse train, which is a sequence of dirac delta distributions spaced by the sampling period $T_s$.

With this in mind and remembering that multiplicating in the time domain is equivalent to doing convolution in the frequency domain, it's possible to see that the frequency representation of the digitalized signal will be the convolution of the frequency representation of the sampled analog signal with the frequency representation of the impulse train.

As it turns out, the frequency representation of an impulse train spaced by $T_s$ is an impulse train spaced by the sampling frequency $f_s$. Also, the convolution of a dirac delta with any function $g$ is just $g$ itself. By displacing the delta function by $f_s$, we also displace the result of the convolution by $f_s$.

So, the spectrum of the digitalized signal will be an infinite sum of copies of the spectrum of the analog signal distant of each other by $f_s$. Since the spectrum is bilateral, if your analog signal has frequency components above $f_s/2$, those components will end up summed from two different replicas in the spectrum of the digital signal, and this is in fact the aliasing effect.


Impulse response of the system image processing In case of the signal, when we plot an impulse response we have our x-axis - samples (or time), y-axis - amplitude, but in case I work with image processing?
x-axis is pixels? and y is still amplitude?

You imagine several representations. A first one is the standard image representation: flat in the $(x_1,x_2)$ and a color scale  of the intensity value (left plot). A second one is a 3D vision as a topographical or elevation map, with homogeneous color (middle plot). The visualization might be a bit more complicated. A third one combine the two preceding aspects (right plot). The latter helps understanding some mathematical morphology operators (like the watershed).

L = membrane(1,100);
f = figure;
ax = axes;
subplot(1,3,1)
imagesc(L);
subplot(1,3,3)
s = surface(L);
s.EdgeColor = 'none';
view(3);axis tight
subplot(1,3,2)
s = surface(L);
s.EdgeColor = 'none';
s.FaceColor = 'b';
view(3);axis tight


Transient response of system with single pole $0 \le p < 1$ $G(z) = \frac{1-p}{z-p}$

  If the value of p satisfies $ 0 \leq p < 1$ there are no oscillations
  in the transient response.

Question: Why is that $\uparrow$  true? I know roughly what a transient response is but how is the relation between the position of the pole and the nature of the transient response?

Edit: In another paper by the same authors:

  Enforcing the stability of the controlled system means ensuring that
  the pole p is non-negative and less than 1.

Now I'm completely out, stability does not depend upon the position of the pole as long as they are inside the unit circle of the z-plane I have thought?

Paper 2:

http://www.iste.uni-stuttgart.de/fileadmin/user_upload/iste/zss/publications/supplementaryMaterial/2014-ICSE-Control-PaperReviewCopy01.pdf

Look at the sentence short above eq. (10) on page 5

Paper 1 is not available on the internet as fas as I know. It's Antonio Filieri et al. "Software Engineering Meets Control Theory"

If you have a single real-valued pole, you get a term $k\,p^{(n-n_0)}u[n-n_0],n\ge 0$, with some constant $k$ and some delay $n_0\ge 0$ in the system's impulse response ($u[n]$ is the unit step). Clearly, if $|p|\ge 1$ the transient will never settle. If $-1<p<0$, the transient will alternate in sign (i.e., "oscillate"), like e.g.,

$$\left(-\frac12\right)^nu[n]=\left\{1,-\frac12,\frac14,-\frac18,\dots\right\}\tag{1}$$

So for a transient that decays (i.e., for the system to be stable) and that does not oscillate you need $0\le p<1$.

Note that if you only have one real-valued pole $p$, stability implies $$-1<p<1\tag{2}$$ If you exclude systems with oscillating impulse responses (i.e., with a negative pole), you're left with the requirement $0\le p<1$. But, as you've noticed, this requirement is not only dictated by stability (which only requires $(2)$), but also by the requirement of having no oscillations.


Basic question about trigonometric series and transforms thereof I would like to know the relation between the parameters $\{\omega_k,A_k\;|\;k\in\mathbb{Z}\}$ of a series $\sum_kA_ksin(\omega_kx)$ and a related series, for example, $\sum_kA_k^2sin^2(\omega_kx)$.

I would also like to know why a multiplicity of peaks appears in the FT when the components of a series are raised to some power? $N$ mirrored sets of peaks with characteristic spacing on the frequency axis and scaling on the amplitude axis are manifest when taking the FT of $\sum_kA_k^Nsin^N(\omega_kx)$. I am interested to know the relation between the parameters of these $N$-fold multiplicities as well. For instance, the amplitudes and the frequencies appear to scale geometrically (e.g., $10$ Hz, $5$ Hz, $2.5$ Hz, $1.25$ Hz).

There is no general relationship between the Fourier transform of $f$ and that of $g(f)$ where $g$ is an arbitrary function. The Fourier transform does have the linearity property, so if $g$ is something simple like an affine transform, then the same linear relationship applies to their transforms $F$ and $G$.

With respect to your second question, where $h = \sum_k \left(A_k \sin(\omega_k x)\right)^N$, the presence of more than $k$ peaks in the Fourier transform of $h$ is easily explained using the power-reduction trigonometric identity:

$$
\sin^n(\theta) =
\begin{cases}
\frac{2}{2^n} \sum_{k=0}^{\frac{n-1}{2}} (-1)^{(\frac{n-1}{2}-k)} \binom{n}{k} \sin{((n-2k)\theta)}, & n \text{ is odd} \\
\frac{1}{2^n} \binom{n}{\frac{n}{2}} + \frac{2}{2^n} \sum_{k=0}^{\frac{n}{2}-1} (-1)^{(\frac{n}{2}-k)} \binom{n}{k} \cos{((n-2k)\theta)}, & n \text{ is even}
\end{cases}
$$

(the above is shamelessly borrowed from Wikipedia)

So, when you raise a sinusoid to a power, the result can be expressed as a weighted sum of sinusoids at different frequencies, where the number of individual terms is related to the power. That's why you see additional peaks in the spectrum of $h$.

You can come up with a more general relationship for some cases by taking advantage of the multiplication property of the Fourier transform. That is, if $g = f \cdot e$, then its Fourier transform is $G = F * E$ (where $*$ indicates convolution). You could apply this relationship repeatedly to the sinusoid raised to a power to derive the same result as above.


Obtaining the FFT of a data set collected through an accelerometer I'm using MATLAB to analyse a very large sample of acceleration data collected using the MPU9150 IC. My aim is to plot the FFT for this large sample of data.  The daya is of three axes; $x,y,z$ and I have successfully imported the data into matlab as vectors.

As I am new to MATLAB I have taken over month to plot the FFT graph but I believe its highly dependent on the smoothing factor. I wish to first illustrate my approach and then showcase my full code for your kind reference. Sampling frequency of of the accelerometer is 10 Hz.

Approach in exact order.

Import the raw acceleration vector with $x,y,z$ axes Determine
butterworth Filter parameters
Apply butter-worth filter using filtfilt()
use smooth() function with a smooth factor. ( NOTE: smoothing factor
effects FFT shape significantly)
obtain FFT using abs(fft(acceleration_DC_shifted_Filtered_Vector))
smooth the FFT line using smooth function (this smoothing has no
effect on FFT there fore can be skipped if required.)
plot FFT

My code is below,

%importing data from workspace.
     accelerationd2 = Xaccelday2;
% Setting basic Parameters. Number of samples should be same in all day files.         number_of_samples = length (accelerationd2); sampling_frequency = 10;

%Setting filtering paramters.
          cut_off_frequency = 5;
          wn = cut_off_frequency/sampling_frequency;

           order_of_butterworth_filter = 1;
              [b1,a1]= butter(order_of_butterworth_filter,wn,'low');

%Filtering data from all ten days.
                  accelerationd2_filtered = filtfilt(b1,a1,accelerationd2);

%DC shift removal`

              accelerationd2_smoothned= smooth(accelerationd2_filtered,36.5);%acceleration 1filtered was passed into all smoothfunctions
          accelerationd2_DC_shifted = accelerationd2_filtered-accelerationd2_smoothned;

%Computing FFT
          accelerationd2_FFT = smooth(abs(fft(accelerationd2_DC_shifted)),8000); %1000 is ballpark
          bin_vals = (0 : number_of_samples-1);
          fax_Hz = bin_vals*sampling_frequency/number_of_samples;
          N_2 = ceil(number_of_samples/2);

%Detecting peaks
          figure;
          plot((fax_Hz(1:N_2)),(accelerationd2_FFT(1:N_2)),'r');
          hold on
          legend('Day-02');

%[maxtabDay1, mintabDay1] = peakdet(accelerationd1_FFT(1:N_2), 15, fax_Hz(1:N_2));
          [maxtabDay2, mintabDay2] = peakdet(accelerationd2_FFT(1:N_2), 15, fax_Hz(1:N_2));

%This is where peaks are dressed with markers. Markers must have the same color of the line. Markers can be different symbols. Types of markers that can be used : ^ , *
           plot(maxtabDay2(:,1), maxtabDay2(:,2), 'r*');
          xlabel('Frequency (Hz)');
          ylabel('Magnitude (dB)');
          title('Single-sided Magnitude Spectrum (Hz)');
          %hold on
          hold off
%axis tight

My problem is that the FFT changes highly as the smoothing factor in

 accelerationd2_smoothned= smooth(accelerationd2_filtered,36.5)

is changed. Any help would be so much appreciated.

Is it my thinking that there is only one correct FFT and the FFT should not change depending on the process I obtain it from?

Another important reason to note is that all samples have a DC shift. Which is why the smooth function has been used. If it is not used, the time domain is not brought about $y=0$-axis before inserting to the FFT function.

I would appreciate if a better procedure can be explained that can be used to obtain the FFT through MATLAB.

P.S.
It would be nice to have a method I can obtain the FFT for my data set. I collected this data about the motion of an particular animal. And I need to generate the FFT for its movement. I have many records but not a proper method. My data goes like,
$X_{\rm acceleration}, Y_{\rm acceleration}, Z_{\rm acceleration}$.
But there can be DC shifts seen on the time domain because the animal because the collar which has the logger was not perfectly in the middle. Therefore I must remove these DC components caused by gravity.

From you code it appears you do a low pass butterworth filter:

[b1,a1]= butter(order_of_butterworth_filter,wn,'low');
accelerationd2_filtered = filtfilt(b1,a1,accelerationd2);

Then another low pass moving average filter which you minus from the previous.

accelerationd2_smoothned = smooth(accelerationd2_filtered,36.5);
accelerationd2_DC_shifted = accelerationd2_filtered - accelerationd2_smoothned;

Yes, this second part will remove DC, but I think this is the first thing you should do. That is, use a HIGH pass butterworth filter (with a few more orders perhaps)

cut_off_frequency = 5;
wn = cut_off_frequency/sampling_frequency;
order_of_butterworth_filter = 4;
[b1,a1]= butter(order_of_butterworth_filter,wn,'high');
accelerationd2_filtered = filtfilt(b1,a1,accelerationd2);

Then apply the FFT. If there is two much high frequency noise, you can apply a low pass butterworth filter to the above output.

As for you analysis, the FFT uses the entire signal to compute. So any frequency magnitudes represent the average over the whole time series. That is, there is no time information in there. Since you have a large sample of animal data, I imagine you would want to see how things change over time, so I suggest instead using the STFT or some other wavelet transform.

In MATLAB you might get good results out of the box using spectrogram, no pre-filtering required.

window = 1000;
overlap = round(0.9*window);
s = spectrogram(accelerationd2,window,noverlap)


Understanding the units of wavelet time & frequency resolution While using a continuous wavelet transform for my research project, I came up with some questions.
Having a set of equations for a Morlet wavelet,

$\psi(t)=\sqrt{\frac{2}{\pi}}e^{-t^2/2}e^{i6t}$
$\hat{\psi}(\omega)=2e^{-(\omega-6)^2/2}$

I have calculated the time and frequency resolution (i.e., $\Delta t$ and $\Delta \omega$, respecitvely) such that

$\sigma_\psi^2=\frac{\int_{-\infty}^{\infty}(t-t_0)^2|\psi(t)|^2dt}{\int_{-\infty}^{\infty}|\psi(t)|^2dt}=\frac{1}{2}$ and $\sigma_\hat{\psi}^2=\frac{\int_{-\infty}^{\infty}(\omega-\omega_0)^2|\hat{\psi}(\omega)|^2d\omega}{\int_{-\infty}^{\infty}|\hat{\psi}(\omega)|^2d\omega}=\frac{1}{2}$ where $t_0=0$ and $\omega_0=6=2\pi f_0$
$\Delta t=\sigma_t=s\sigma_\psi$ and $\Delta \omega=\sigma_\omega=\frac{\sigma_\hat{\psi}}{s}$ where the scale $s$ is approximated with $\frac{f_0}{f\cdot \text{sampling period}}$

I have been trying to calculate the Heisenberg box as the intervals $[t-\frac{\sigma_t}{2},t+\frac{\sigma_t}{2}]$ and $[\omega-\frac{\sigma_\omega}{2},\omega+\frac{\sigma_\omega}{2}]$, but the computation resulted with counterintuitive values, with $\Delta t$ being extremely large and $\Delta \omega$ being extremely small for almost every $t$ and $\omega$ values I plugged in.
As I read some previous posts (including here), it seems the problem may be coming from the units of the resolutions. I read that the units of $\Delta t$ and $\Delta \omega$ are

$[\Delta t] = \frac{\text{samples}}{\text{cycles}\cdot\text{radians}}$
$[\Delta \omega] = \frac{\text{cycles}\cdot\text{radians}}{\text{samples}}$

If I convert samples to seconds, the resulting intervals for the Heisenberg box seem to make sense. However, even if I change samples to seconds, the units still include cycles in them, and I do not know how I should interpret the values qualitatively.
My main questions therefore are (1) if the equations defining $\Delta t$, $\Delta \omega$, and $s$ are mathematically correct, and (2) how I should change the units of time and frequency resolution so I can interpret them in seconds and Hz.
Thank you very much.

Confusion seems rooted in notation:

In $s = f_0 / (f \cdot \text{sampling rate})$, $f_0$ is the center frequency of mother $\psi$, and $f$ is the center frequency of the scaled $\psi$, i.e. $\psi_s$.
In $\sigma_{\hat\psi}=$ ..., $\omega_0$ is the (radian) center frequency of the wavelet of interest, or $\psi_s$, rather than $\omega_0=2\pi f_0$.
However, you can compute sigma of scaled from sigma of mother via: $\sigma_{\hat\psi_s} = \sigma_{\hat\psi} / s$ (see under 4.54).
$\text{sampling rate}$ is only relevant in discretized computations and must get canceled to yield physical units (Hz, sec, etc), so if you seek only physical use $s = f_0 / f_s$ (adjusted $f\rightarrow f_s$ to avoid confusion; recall that $f$ (or $\omega$) is what we integrate over, rather than plug in)
Do not involve $\text{sampling rate}$ in defining $\psi, \hat\psi$, scaled, or computing $\sigma$, or effectively most formulas you come across unless they explicitly include it, it's only to be used after computation, e.g. you first find $\sigma_{\hat\psi}$ without $\text{sampling rate}$, and then rescale as in $(4)$ here.
You commented "$\sigma_\omega$ Hz", but it's actually "rad-Hz"; to get Hz do
$f=\omega / (2\pi)$.

Let me know if this makes sense (if not, again with a specific example).


Understanding the $\mathcal Z$-transform I was studying $\mathcal Z$-transforms and found pretty good material on the topic, though I feel I do not have a proper understanding of the concept. Could someone help me clarify this?

I know that we apply $\mathcal Z$-transforms to sequences and convert them to the $z$ domain and find the values of $z$ for which we get zeros and poles and also find the region of convergence. What I didn't get is:

Why are we multiplying some $z$ values to the given sequence and finding for what values of $z$ the result gets pretty huge or dampens out?
What does poles and zeros have to do with the given sequence? It's not like we are applying the transform to a system and finding the system's response, it's just a sequence.
And they also say that its a frequency domain transformation. How do we realise frequency here?

Consider a liner discrete-time system. Assume we can define it in terms of an input-output relation as follows (you can assume a more general model but it is enough for our purpose):

$$a_0y[n]+a_{1}y[n-1]+\cdots+a_{N}y[n-N]=b_0x[n]+b_{1}x[n-1]+\cdots+b_{M}x[n-M]\tag{1}$$

When the coefficients $\{a_i\}$ and $\{b_i\}$ are constant, we call it a finite-order constant-coefficient ordinary difference equation.
It explains the current output $y[n]$ in terms of a weighted sum of the current and past inputs and the past outputs:

$$y[n]=\frac{-a_{1}y[n-1]-\cdots-a_{N}y[n-N]+b_0x[n]+b_{1}x[n-1]+\cdots+b_{M}x[n-M]}{a_0}$$

It is very similar to a differential equation in continuous time.

Dealing with such equation to find the output of the system can become complicated. It is a recurrence relation. It defines the values recursively and for arbitrary inputs it is not straightforward to express the output in a closed-form representation.

How to deal with it in an easier way?

Consider the following transform (assume the sum exists):
$$\mathcal{Z}(x[n])=\sum_{n=-\infty}^{+\infty}x[n]z^{-n}$$
You are right, it is just a transform that accept $x[n]$ and gives you $X(z)$.
But it has a useful property. Lets calculate the $z$-transform of $x[n-\alpha]$:
$$\begin{align}\mathcal{Z}(x[n-\alpha])&=\sum_{n=-\infty}^{+\infty}x[n-\alpha]z^{-n}\\[10pt]
&=\sum_{n'=-\infty}^{+\infty}x[n']z^{-(n'+\alpha)}\tag{2}\\[10pt]
&=\left(\sum_{n'=-\infty}^{+\infty}x[n']z^{-n'}\right)z^{-\alpha}\\
&=\left(\mathcal{Z}(x[n])\right)z^{-\alpha}
\end{align}$$
where in $(2)$ I changed the variable $n'=n-\alpha \Rightarrow n=n'+\alpha$.
Assume $\mathcal{Z}(x[n])=X(z)$. We have seen that a property of $z$-transform is: $$\boxed{\mathcal{Z}(x[n-\alpha])=z^{-\alpha}X(z)}$$

Applying this to the difference equation $(1)$ makes it (note that it is a liner transform and can be applied term by term):

$$a_0Y(z)+a_1z^{-1}Y(z)+\cdots+a_Nz^{-N}Y(z)=b_0X(z)+b_1z^{-1}X(z)+\cdots+b_Mz^{-M}X(z)$$

Note that all terms are expressed in terms of $Y(z)$ and $X(z)$. So we can say that

  $z$-transform reduces usually difficult-to-handle recurrence relations to much easier to manipulate algebraic relations.

With an algebraic relation, we can factor the terms:

$$Y(z)\left(a_0+a_1z^{-1}+\cdots+a_Nz^{-N}\right)=X(z)\left(b_0+b_1z^{-1}+\cdots+b_Mz^{-M}\right)$$
and consequently,

$$Y(z)=X(z)\frac{b_0+b_1z^{-1}+\cdots+b_Mz^{-M}}{a_0+a_1z^{-1}+\cdots+a_Nz^{-N}}$$

where $$H(z)\triangleq\frac{b_0+b_1z^{-1}+\cdots+b_Mz^{-M}}{a_0+a_1z^{-1}+\cdots+a_Nz^{-N}}$$
is called the system function. So we distilled all the system into an algebraic expression as $H(z)$. Hence, the zeros and poles of $H(z)$ have a direct impact on the output.

There are also other reasons that make $z$-transform and $H(z)$ important, and we would prefer working with the $z$-domain rather than time-domain.

Assume the input to the system is $x[n]=z_{k}^n$, where $z_k$ is just a complex number (don't mix it with the $z$ variable in the $z$-transform at this moment! Also there are some restrictions on $z_k$ as we will see next).
The output of the LTI system can be calculated by convolution of the input and the impulse response $h[n]$:

$$\begin{align}
y[n]&=h[n]*x[n]\\
&=\sum_{m=-\infty}^{+\infty} h[m] x[n-m]\\
&=\sum_{m=-\infty}^{+\infty} h[m] z_k^{n-m}\\
&=z_k^n \sum_{m=-\infty}^{+\infty} h[m] z_k^{-m}\\
&=z_k^n H(z_k)
\end{align}$$

What does it really mean?

The output of the system is an scaled version of the input. So the signal $z_k^n$ passes through the system and the only impact on it is a scaling. Note that the scaling factor is the system function when evaluated at $z_k$. Why is it important?

If we could decompose an arbitrary signal (not of the form $z_k^n$) into several components of the form $z_k^{n}$, then we can easily calculate the output of the system to each of the components, and then add them together to find the overall output (benefiting from the superposition property of LTI systems).

How can we decompose the signal $x[n]$ into components of the form $z_k^n$?
It is easy as the following (assume the sum exists for now)

$$\left(\sum_{n=-\infty}^{+\infty}x[n]z_k^{-n}\right)z_k^n$$

Let's refer to this component as $X(z_k)z_k^n$. So the output to this component is $y_k[n]=H(z_k)X(z_k)z_k^n$. We can have as many components as possible. Potentially, we have thousands of $z_k$. So we refer to all of them by a variable $z$ and the component-wise output becomes $H(z)X(z)z^n$. The overall output $y[n]$ is in form of an integral (since we need to add them up). It can be shown that $y[n]$ is given by the following contour integral:

$$y[n]=\frac{1}{2\pi j}\oint_C\frac{H(z)X(z)z^n}{z}dz$$

where $C$ is a counterclockwise path inside the ROC (where $X(z)$ and $H(z)$ exist) that encircles the origin, which is the same as the inverse $z$-transform of $X(z)H(z)=Y(z)$.

Regarding your second question, in the definition of $z$-transform (when it exists on the complex unit circle) if you choose $z=e^{j\omega}$ it becomes the definition of the Discrete-Time Fourier transform Hence is the frequency-domain interpretation.


homework impulse response of system Problem

A student is given an unknown DT system. Upon entering a test signal x_tst[n] = [-1, 1] into the system, an output y_tst[n] = [0,-1,-1,2] is observed. Determine the impulse response of the system.

note: bold digit - x[0] or y[0]

What I tried:

Something similar to this:

Choices
A.) h[n] = [0,-1,-2]

B.) h[n] = [0, 1,-2]

C.) h[n] = [0, 1, 2] -->correct answer

D.) h[n] = [0,-1, 2]

How do you solve these kind of problems?

First of all those sequences you provided must have a time index reference as well. I assume them to be causal sequences and $n=0$ to be the first sample.

Then you can use different approaches to find the impulse reponse for example using the Z-transforms provide: $$ H(z) = \frac{Y(z)}{X(z)} $$
Then you have to convert $H(z)$ back into time to find the impulse response $h[n]$.

But for this simple problem, I would instead use a time-domain approach directly and compute the samples of $h[n]$ by the following procedure:

Given $x[n] = -\delta[n] + \delta[n-1]$ of length $N_x=2$

and $y[n] = -\delta[n-1] -\delta[n-2] + 2\delta[n-2]$ of length $N_y=4$

(note that the first sample $y[0]$ of $y[n]$ is zero but here we consider $y[0]$ to be in the support-domain of $y[n]$ no matter what its value is, the result won't change for this problem)

1- Assuming an FIR $h[n]$ of length $L = N_y - N_x + 1 = 4-2+1 = 3$ samples.

2- Now propose the causal $h[n]$ to be $$h[n] = a\delta[n] + b\delta[n-1] + c\delta[n-2]$$

3- Now apply the convolution as the sums of shifted $h[n]$ rows as follows:

$$ y[n] = h[n]\star x[n] = x[0]h[n] + x[1]h[n-1] $$
$ $

$~~~~~~~~~~~~~~~ -a ~~ -b ~~ -c$

$~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~a ~~~~~~~b ~~~~~c$

$ + ~~~ ------------$

$y[n]=~ -a ~~(a-b) ~~ (b-c) ~~c$

$y[n]=~ -0 ~~~~(-1) ~~~~ (-1) ~~~~~2$

from which you can easily solve $a,b,c$ to be

$c = 2$ , $b-c=-1 => b = 1$ and $a=0$, therefore the impulse reponse is found to be $$ h[n] = 0\delta[n] + 1\delta[n-1] + 2\delta[n-2]$$


How to keep phase response of IIR filter? We know that IIR digital filter has nonlinear phase response in general.

How can I linearize the phase response of IIR filter without altering its magnitude response?

To my knowledge, you can't linearize it arbitrarily, but that being said, you can obtain approximately linear phase by processing the output through a cascade of all-pass filters.

An all-pass filter used this way is usually called a phase compensator or a phase equalizer. A single pole-zero all-pass IIR filter has the following structure
$$ H(z) = \frac{p - z^{-1}}{1 - p z^{-1}} $$
where $p$ controls the location of the pole/zero. You can optimize a cascade of filters of this type to equalize your phase response using least squares or some other filter design technique.

Note that this technique reduces the primary advantage of choosing an IIR filter to begin with. Usually, the point in choosing an IIR filter over an FIR filter is to reduce computational complexity, so performing the phase equalization by cascading all-pass sections might result in a situation where using a linear phase FIR filter in the first place would have been more economical.


What does the $n$ do in $1/(j \omega)^n$ I'm currently working through a research paper on beam forming. In this paper a magnitude compensation is introduced to compensate for frequency dependency. Due to other calculations low frequencies are lower in magnitude than high frequencies.

To compensate for the frequency dependency they multiply the signal with:
$$
 \frac{1}{(j \omega)^n}
$$

I know when $n = 1$, the formula is just a simple integrator. But I don't understand what happens when $n$ is higher, and how I should implement it in MATLAB.

What do you mean with Beam in this context.... Particle Beams in LINACs? Structural Beams in Buildings?.

$j\omega$ is a phasor term from frequency analysis.

$\frac{1}{j\omega}$ is the integration as @msm showed.

If you want to integrate two times, just use cumtrapz two times.

But as you said, this is only a "term" complementing a frequency model, which it would be really good to have here in order to help you more... You never simulate that term alone.


Call recognition using Digital Signal Processing We are looking for writing applications which would recognize Bird calls and Frog calls using Sound processing.

EDIT: I want to identify the species of bird / frog and do not want to differentiate whether the call is from a bird or a frog. There are about 1300+ bird species in India and about 300+ frog species which we are looking to identify based on the call.

I have a background of iOS application development and C / C++ programming but have never worked on Digital Signal Processing.

We do have few researchers who would collect multiple calls of birds and frogs and categorize it as per species. These calls can be used as reference purpose (?) to identify the bird or frog species by recording the call and comparing this call in user's device.

I do not know the technical feasibility of it, but I would like to know:

Are there any open source solutions which can be extended / used to achieve this?
Is this really possible to achieve with some 80% to 90% accuracy?
If open source solution does not exist, how to go ahead with the implementation?

I understand that the scope of this question is too broad and wide which should be ideally avoided in StackExchange. But would like to have your opinion.

Basically it is a classification problem. You can want to classify the incoming voice.

There are many ways to tackle it. One of them would be to extract features ( by features, I mean frequency, pitch etc) and label it accordingly(if it is of bird or frog).

Then shortlist the features, which will help you differentiate between these two classes(bird or frog) and make a dictionary.

In your dictionary , it will be like this

Class
   features
                  Pitch   Frequency etc etc

Frog        (P)  = 5     (F)=  10

Bird         (P)   =   1 (F) =       20

Whenever, your new sound comes in, you compute it features, and based on which feature it closely corresponds to (or there are many methods to do it) you choose the class.


STFT and DWT (Wavelets) STFT can be successfully used on sound data (with a .wav soundfile for example) in order to do some frequency-domain modifications (example : noise removal).
With N=441000 (i.e. 10 seconds at sampling rate fs=44100), windowsize=4096, overlap=4, STFT produces approximatively a 430x4096 array (first coordinate : time frame, second coordinate : frequency bin). Modifications can be done on this array, and reconstruction can be done with overlap-add (*).

How is it possible to do a similar thing with wavelets ? (DWT), i.e. get a similar  array  of shape a x b, with a time frames, and b frequency bins,  do some modification on this array, and at the end, recover a signal ? How ? What is the wavelet equivalent to overlap-add ?
What would be the Python functions involved here (I haven't found an easy example of audio modification with pyWavelets...) ?

(*) : Here is the STFT framework that can be used :

signal = stft.Stft(x, 4096, 4)    # x is the input
modified_signal = np.zeros(signal.shape, dtype=np.complex)

for i in xrange(signal.shape[0]):    # Process each STFT frame
    modified_signal[i, :] =  signal[i, :] * .....  # here do something in order to
                                                   # modify the signal in frequency domain !
y = stft.OverlapAdd(modified_signal, 4)   # y is the output

The goal is to find a similar framework with wavelets.

The short-time Fourier transform is generally a redundant transformation, usually implemented with the same subsampling over every frequency. If the window is well chosen, it is complete: you can invert it and recover any initial signal.

Since it is redundant and complete, it has many perfect inverses. It can be implemented and understood using more generic tools: (oversampled) complex filter banks. Given a window type and length plus the overlap provides you with an analysis filter bank for which you can compute whether it is invertible or not. If it does, you can compute a natural inverse, and optimized inverses too. The overlap-add is just one on the many potential inverses, probably the most common, which often restricts the window choice.

Standard discrete wavelet transforms are filter banks too, with the difference that the subsampling is not the same in each frequency band (or more properly scale). This turns into uneven lengths for each scales. Yet, there exists redundant wavelet implementations that yield "a rectangle array" of coefficients you can work with. The best known schemes are called under different names: shift-invariant or time-invariant wavelets, undecimated wavelets, stationary wavelet transform (SWT), and sometimes cycle-spinning. Its standard reconstruction involves steps similar to overlap-add, except they are more "embedded" due to the different sampling factors over scales. You can use them with any discrete wavelet from a library, or even by designing your own wavelet. The reason is that standard discrete wavelets were designed with non-redundancy in mind, which restricts the choice of wavelets. With redundancy, you wavelet choice increases, as the constraints to fulfill are less stringent. The "ultimate" avatar is the continuous wavelet transform, which admits "almost" every inverse synthesis wavelet. My last sentence is quite lousy, I hope you get the meaning: when a square matrix is invertible, it has only one inverse. When a "rectangular" matrix is left-invertible in a generalized way, it has several inverses.

There seems to be a python implementation of the Stationary Wavelet Transform. You can find a few references in the 2.3.4. Translation invariant wavelets chapter of the linked paper.

It is generally much more robust for detection, denoising or restoration in practical applications (geophysics, non-destructive testing, ultrasounds, vibrations).