Pitch Detection / Harmonic Product Spectrum - What's wrong? I've whipped together a small example that tries to use the Harmonic Product Spectrum Algorithm to extract the pitch of a simple sine wave. I'm not sure about my implementation and if i understood everything completely, that's why i'm hoping for some helping hand here.

The code is here: https://gist.github.com/akuehntopf/4da9bced2cb88cfa2d19
(You'll need JTransform and JavaFX to compile, though).

Please see what happens with the data during the processing using the
following graphs. Note: I'm using a sine wave with 160Hz and a sample rate of 8000. The duration of the sound data is 1 second.

First the sine wave is displayed (see first image).

Second, i'm applying a hamming window to the data (second image)

Then i calculate the power spectrum of the data as per

private float[] powerSpectrum(float[] window) {
    float[] powerSpectrum = new float[window.length];
    float[] fftBuffer = new float[window.length * 2 + 1];
    System.arraycopy(window, 0, fftBuffer, 0, window.length);
    FloatFFT_1D fft = new FloatFFT_1D(window.length);
    fft.realForward(fftBuffer);

    for (int i=0; i < fftBuffer.length / 2 - 1; i++) {
        float real = fftBuffer[2*i];
        float imag = fftBuffer[2*i+1];

        powerSpectrum[i] = (float)Math.sqrt(real*real + imag * imag);
    }

    return powerSpectrum;
}

Which gives me the spectrum in the third graph.

Observe that the peak is approximately at the right spot already, but we go for HPS, so continue (i'm not really sure about everything from here on):

Next, i'm multiplying the signal with its compressed form several times (with increasing compression)...

// 4. Compress
float[] spectrumCopy = new float[spectrum.length];
System.arraycopy(spectrum, 0, spectrumCopy, 0, spectrum.length);
for (int compression = 2; compression < 4; compression++) {
    for (int i = 1; i < spectrum.length; i++) {
        spectrum[i] = spectrum[i] * getCompressedSample(spectrumCopy, 1, compression, i);
    }
}

it uses the function

private float getCompressedSample(float[] buffer, int offset, int compression, int loc) {
    if (offset + loc * compression < buffer.length) {
        return buffer[offset + loc * compression];
    }
    return 0;
}

i found it somewhere on the interwebs. But my understanding is that we compress the power spectrum n times while taking only every 2, then only every 3 and so on samples from the power spectrum. The original spectrum is then multiplied with each of those compressed spectra.

The result i get is the fourth graph

Observe that the peaks have changed.

From my understanding the next step would be to find the bin with the highest peak, interpolate (i use quadratic interpolation) and recalculate the pitch frequency using this formula:

private float getFrequencyForIndex(int index, int size, int rate) {
    float freq = (float)index * (float)rate / (float)size;
    return freq;
}

However this gives me wrong results. I'm pretty much stuck. Any help with this topic would be very much appreciated! Thanks in advance!

For the next step i'm hoping to be able to extract pitch frequency data from a musical instrument (guitar/ukulele), but of course first things first :-)

The Harmonic Product Spectrum pitch estimation method, using pure multiplication weighting, only works well for signals that have a full set of harmonics with sufficient magnitudes.  Before using HPS, one should either validate that the signal in question meets that criteria, or add a suitably large non-zero floor to all harmonic magnitudes being multiplied together.  Otherwise a pure sinewave with zero harmonics will simply disappear from the result after being multiplied by a zero magnitude.

Guitar and ukulele sounds usually have a far richer set of harmonics than test sine waves.


Fourier transform of time division I know that Fourier transform of $t^n f(t)= i^n \frac{d}{d\omega^n} F(\omega)$.
But does this work when $n<0$?
Is there any direct relation to compute the Fourier transform of $\frac{f(t)}{t}$?

A unit step in the frequency domain $u(f)$ has the following inverse Fourier transform:
$$\mathcal{F}^{-1}\{u(f)\}=\frac12\delta(t)-\frac{1}{2\pi jt}\tag{1}$$
This follows from the well-known Fourier transform of the time-domain unit step $u(t)$, and from the basic Fourier transform property
$$\mathcal{F}\{x(t)\}=X(f)\Longrightarrow\mathcal{F}^{-1}\{x(f)\}=X(-t)\tag{2}$$
Given $x(t)$ and its Fourier transform $X(f)$, we have
$$(X\star u)(f)=\int_{-\infty}^fX(\nu)d\nu\tag{3}$$
where $\star$ denotes convolution.
Using $(1)$, the inverse Fourier transform of $(3)$ is
$$x(t)\cdot \mathcal{F}^{-1}\{u(f)\}=\frac{x(0)}{2}\delta(t)-\frac{x(t)}{2\pi jt}\tag{4}$$
Consequently,
$$\mathcal{F}\left\{\frac{x(t)}{t}\right\}=2\pi j\left[\frac{x(0)}{2}-\int_{-\infty}^fX(\nu)d\nu\right]\tag{5}$$
As an exercise, you can try to verify $(5)$ with $x(t)=\sin(2\pi f_0t)/\pi$ (resulting in an ideal low-pass with cut-off frequency $f_0$), and with $x(t)=1/\pi$ (resulting in an ideal Hilbert transformer).


Comparison of results with Direct form1 and Direct form2 implementation of IIR filters The results obtained from Direct form1 and Direct form2 implementations (in C language) of IIR (BiQuad) filters are a little different. In the implementation, i have made the b0 coefficient as 1. So, I am just using 4 coefficients (a2,a1,b2,b1) in the implementation. The maximum error between the obtained outputs from these two implementations is in the range of 0.01 for the ITU-T standard audio test signal. So, I just want to know whether this error is because of the implementation, or any other reasons? Will the Outputs of both implementations be the same, always (i.e. error is almost 0)?

In the implementation, I am using around 8 stages of filtering. (i.e. Cascaded Biquad filters)

The referred structure forms as attached below.
enter image description here

Will the Outputs of both implementations be the same, always (i.e. error is almost 0)?

That depends A LOT on the numerical precision of your data representation and mathematical operations and on the filter itself, specifically the location of the poles.

If you use double precision in a communication application, there should be little to no difference. If you are trying to do audio processing with 16-bit fixed point format, the errors induced by Direct Form II are so large that it's basically unusable.


Camera calibration/pin hole camera model and working out 3d position I have a calibrated camera and have the intrinsic parameters. I also have the extrinsic parameters relative to a point (the world origin) on a planar surface in the real world. This point I have set as the origin in the real world coordinates [0,0,0] with a normal of [0,0,1].

From these extrinsic parameters I can work out the camera position and rotation in the world plane 3d coordinates using this here: http://en.wikipedia.org/wiki/Camera_resectioning

Now I have a second point which I have extracted the image coordinates for [x, y]. How do I now get the 3d position of this point in the world coordinate system?

I think the intuition here is that I have to trace a ray that goes from the optical center of the camera (which I now have the 3D position for as described above), through the image plane [x,y] of the camera and then through my real world plane which I defined at the top.

Now I can intersect a world coordinate 3d ray with a plane as I know normal and point on that plane. What I don't get is how I find out the 3d position and direction when it leaves the image plane through a pixel. It's the transformation through different coordinate systems that is confusing me.

If you have the extrinsics then it is very easy. Having extrinsics is the same as having "camera pose" and the same as having the homography. Check this post in stackoverflow.

You have extrinsics, also called camera pose, which is described as a translation and a rotation:

$\displaystyle Pose =\begin{bmatrix}R|t \end{bmatrix} = \begin{bmatrix}R_{11} &R_{12}&R_{13}&t_x\\R_{21}&R_{22}&R_{23}&t_y\\R_{31}&R_{32}&R_{33}&t_z \end{bmatrix} $

You can get Homography from Pose this way:

$\displaystyle   H = \frac{1}{t_z}\begin{bmatrix}{R_{1x}}&{R_{2x}}&{t_x}\\{R_{1y}}&{R_{2y}}&{t_y}\\{R_{1z}}&{R_{2z}}&{t_z}\end{bmatrix}$

Then you can project your 2D points into the corresponding 3D points by multiplying the Homography by the points:

$p_{2D}=\begin{bmatrix}x &y &1\end{bmatrix}\quad$ add $\quad z=1\quad$ to make them homogeneous

$p_{3D}=H*p_{2D}  $

$p= p / p(z)\quad$   Normalize the points


Relationship between the autocorrelations of X(t) and X(nt) Defining:

$X(t)$ WSS random process with autocorrelation function $R_{X}(\tau) = \mathbb{E}[X(t)X(t+\tau)]$.

$Y[n] = X(nT)$ (sampling of $X$ at a rate $\frac1T$) with autocorrelation function $R_Y(\tau) = \mathbb{E}[Y[n]Y[n+\tau]]$.

Which is the relationship between $R_X(\tau)$ and $R_Y(\tau)$?

The answer to the OP's question is more straightforward than rb-j's comments make it out to be.

$\{X(t)\colon -\infty < t < \infty\}$ is a continuous-time WSS random process with autocorrelation function
$$R_X(\tau) = E[X[t)X(t+\tau)], -\infty < \tau < \infty.$$
On the other hand,
$\{Y(n)\colon n \in \mathbb Z\}$ where $Y[n] = X(nT)$ is a discrete-time process whose mean function
\begin{align}
\mu_Y[n] &= E\left[Y[n]\right]\\
&= E[X(nT)]\\
&= \mu_X(nT)\tag{1}\\
&= \mu_X
\end{align}
is a constant since $\{X(t)\}$ is assumed to be WSS and so its mean function $\mu_X(t)$ has fixed value $\mu_X$.  Similarly, the autocorrelation function of $\{Y[n]\}$ is
\begin{align}
R_{Y,Y}[m,n] &= E\left[Y[m]Y[n]\right]\\
&= E[X(mT)X(nT)]\\
&= R_{X,X}(mT, nT),\tag{2}\\
&= R_X((n-m)T)\end{align}
that is, the value of $R_{Y,Y}[m,n]$ depends only on the difference $n-m$ of the arguments $m$ and $n$, and not upon their individual values. We conclude that

  $\{Y(n)\colon n \in \mathbb Z\}$ is a discrete-time WSS random process with mean $\mu_X$ and autocorrelation function
  $$R_Y[m] = E\left[Y[n]Y[n+m]\right] = R_X(mT).$$

Note that just as the random variables $Y[n]$ are regularly spaced samples ($T$ seconds apart) of the $\{X(t)\}$ process, the autocorrelation function $R_Y[m]$ is just regularly spaced samples ($T$ seconds apart) of the autocorrelation function $R_X(\tau)$ of the $\{X(t)\}$ process. Indeed, Eqs. $(1)$ and $(2)$ show that if $\{X(t)\}$ were not WSS, then the mean function $\mu_Y[n]$ would be equal to $E[X(nT)] = \mu_X(nT)$, that is,  the mean function of $\{Y[n]\}$ is just regularly spaced samples of the mean function $\mu_X(t)$ of the $\{X(t)\}$ process while $$R_{Y,Y}[m,n] = E\left[Y[m]Y[n]\right] = E[X(mT)X(nT)] = R_{X,X}\left(mT, nT\right),$$
that is, the autocorrelation function $R_{Y,Y}[m,n]$ is just regularly spaced sample values of $R_{X,X}(t_1, t_2) = E[X(t_1)X(t_2)]$ where the samples are on a square grid.

Since the OP's original question and my answer have resulted in a flood of comments and queries from rb-j that are far too many to be answered in comments, I will add something here to further explain my view of the issues, and then ignore any follow-up queries or comments. If you don't like my answer, down-vote it or write your own answer explaining your view of the matter. Both actions  will be ignored by me.

The OP's question concerns sampling a wide-sense-stationary (WSS) continuous-time process $\{X(t)\colon t \in \mathbb R\}$ to create a discrete-time process $\{Y[n]\colon n \in \mathbb Z\}$ where $Y[n]$ is defined to be $X(nT)$ where $T$ is a known positive constant.  The OP then asks for the relationship between the autocorrelation functions $R_X(t)$ (continuous-time)  and $R_Y[n]$ (discrete-time) of the two processes.

My answer, which I have explained in some detail, to this question,  is that $R_Y[n]$ is just the sampled version of $R_X(t)$:  $$R_Y[n] = R_X(nT)~\text{for all} ~n \in \mathbb Z.$$
This relationship is not a two-way street; given the values of $R_Y[n]$, it is not possible to reconstruct $R_X(t)$. Nor have I claimed anywhere that $R_Y[n]$ and $R_X(t)$ are the same function.  "But, but, but,$\ldots$ what about the sampling theorem? What about ergodicity? What about $\ldots$," you splutter. Please re-read the sampling theorem, especially the hypotheses. If $R_X(t)$ is a strictly band limited function (and nowhere has it been asserted that $R_X(t)$ enjoys this property), then it can be reconstructed from its sample  values provided the sampling rate is adequately high (exactly how high is not relevant). Given only sample values, one can construct infinitely many possible functions that have the specified sample values. If the  sample values are of a function is indeed strictly band limited, then one of these reconstructions will coincide with the function that was sampled. But simply to use the Shannon-Whittaker interpolation formula on $R_Y[n]$ and insist that the result must be $R_X(t)$ is something that boggles my mind, though ymmv. As for ergodicity, the less said the better.


How to estimate the signal-to-noise ratio of a waveform? I have a signal: $f_i(t_i=i\Delta t)$, where $i = 0\ldots n-1$.

The signal seems to vary quickly around a slower varying "trend". I am assuming that the quickly varying part is noise and the slowly varying part is the real signal.

How do I estimate the signal-to-noise ratio (SNR) of the signal?

I guess that if I could decide on a treshold frequency: $\omega_t$ I could use the following expression:

$$S/N=\frac{\displaystyle\int_0^{\omega_t}|F(\omega)|^2}{\displaystyle\int_{\omega_t}^{\infty}|F(\omega)|^2}$$
where $F$ denotes the Fourier transform of $f(t)$.

If the spectrum of the signal and the spectrum of the noise do not overlap, you might be able to integrate or sum the energy in each of the two frequency ranges, and take the ratio.


Are all images 2d or also we have 3d images? Are all the images that we encounter/use in image processing are 2D or do we also have 3D images? If also 3D, then Video will be treated as a 2D image or 3D image?

The term image derive from  Latin  imago. I fully agree with wikipedia saying that it  denotes an artifact

  that depicts visual perception [...] that resembles a subject—usually
  a physical object—and thus provides a depiction of it.

As such, it is a representation of some reality, combining spatial "ordinal" cooordinates (with some natural order, in 2D and 3D for us humans) and  "cardinal" coordinates (luminance/chrominance values, spectra, categories).

The dimension attributed to images to signal/processing people is generally the one of the "ordinal" cooordinates (2D or 3D). To other practitioners who plot 2D images as 3D plots (using elevation for scalar fields in the "cardinal" dimension, non-image-processing people often call 2D image "tridimensional" (bottom-right "Arnold Schwarzenegger" in Terminator), because they look at this artifact in our natural 3D world:

.

Here, $x$ and $y$ are the ordinal dimensions, $f(x,y)$ makes the third. This visualization of images in the resulting topological/topographic map is exploited for instance in the field called "mathematical morphology",  where surfaces and regions are interpreted as geographical objects: valleys, mountains, watersheds, erosions, etc.

If you add (a 1D) time (a very ordinal dimension so far) for motion, then you can add another dimension: standard 2D+t video, 3D video, etc.

If we call for scale-space theory, where one bluntly considers that "objects are meaningful over a certain range of scales only", one can smooth a given image with smoothing kernels to get different images:

and pile them up to yield an additional "scale" dimension. Then, you can extend the discussion to images at arbitrary coordinates, like point clouds, with motion, etc.

Finally, I believe that image data (and their dimensions) can be considered as at least three objects:

how they should be in reality (if such a thing exists),
how they are modeled or described,
how we process them.

and their inherent dimension depends on that. Remember however that one can map two continuous spaces of different dimensions, albeit in a non continuous manner. For discrete spaces, like the ones for classical image/volume coordinates, mappings exist, like zigzag of Hilbert/Peano curves for image compression (JPEG, J2K), so dimensions are a matter of data structuration.


Is it appropriate to use Sobel operators to find the derivatives of an image? Strictly speaking, the $x$ direction derivative should be the difference between the left and right pixel of each pixel. So, then I should be using a $1\times 3$ filter:

\begin{bmatrix}-1&0&1\end{bmatrix}

If that is so, why is the Sobel operator a $3\times 3$?
\begin{bmatrix}-1&0&1\\-2&0&2\\-1&0&1\end{bmatrix}

Is it so that the approximation of the derivative is less susceptible to noise?

Indeed, it adds smoothing in the $y$ direction. The Sobel filter is the separable combination of the centered derivative $[−1,\;0,\;1]$ along $x$, and the $3$-point binomial smoother $[1,\;2,\;1]$ along $y$.


Beginner question on aperiodic samples and DFT behavior I have a strongly aperiodic signal and I get weird artifacts after running a fft, using a low-pass filter, and then doing an inverse fft.  I believe the artifacts are due to the way the signal is assumed to be periodic and "wrapped" at the ends. (Like I said, I'm a complete beginner here.)

I understand the usual way of handling aperiodic signals is via a window function, but applying a window function (say, Hamming) to my data looks like it squishes the data too much and like it would generate worse artifacts.

What sort of seems like it works is to

surround the signal with silence, by prepending and appending a segment of 0s,
doing the fft, applying the filter, etc., and then inverting the fft, then
clipping off the silence segments off the beginning and end.

The result of step 1 looks like applying a rectangular window function to a longer signal, so I think that's ok, but...

Is what I'm doing in any way reasonable?
If so, what are the effects? And what is the best way to determine the length of the segments of 0s? (Currently, 10% of the original signal seems like it's not bad.)
If not, is there another way to do this?

It sounds to me like you are experiencing time-domain aliasing. First, two basic facts:

Recall that if $x[n]$ is of length $N$ and $h[n]$ is of length $M$ then $x[n]\ast h[n]$ is of length $N+M-1$.

Fact 1 does not matter in what domain you perform the convolution/filtering. So, if you modify the DFT of $x[n]$ by "shaping" its spectrum, then there exists an  equivalent filter $h[n]$ in the time domain, which defines an implicit length $M$.

Taking an $N_1$-point IDFT of an $N_2$-point signal for $N_1 < N_2$ results in time-aliasing.

Fact 2 is the dual to frequency aliasing by undersampling. I'll illustrate this as follows: let $x_1[n]$ denote a signal of size $N_1$. Then, if we take an $N_1$-point DFT of $x_1$ and modify its spectrum so the size should be $N_2>N_1$, then taking an $N_1$-point IDFT is equivalent to obtaining the signal $x_2$ given by:
\begin{eqnarray}
x_2[n] = \sum_{p=0}^{\infty}x[n+pN_1]
\end{eqnarray}

This seems to be consistent with what you have mentioned. As you correctly noted, zero padding the original signal is exactly the solution! Essentially just zero-pad such that the terms in the above equation for $p>0$ are all zero-valued. Also, using the ideas above gives you a way to understand the exact amount of zero-padding necessary based on the implicit filter $h[n]$ your frequency domain shaping yields. Good luck!


How can I get the polyphase components of a signal My question, What is the polyhase components of the signal? What does mean a polyphase component for a signal?

Is it true that example :

N = 100;
h = firls(N, [0 .2 .25 1], [1  1 0  0],[.001 .0001]);

% PolyPhase Componets
%First polypase component

for i=1:length(h),
h0(i)=h(i);
end

%% Second polyphase component

for i=1:length(h),
h1(i)=h(i)*exp((-j*2*pi*((i-1)))/4);
end

%% Third polyphase component:

for i=1:length(h),
h2(i)=h(i)*exp((-j*2*pi*(2*(i-1)))/4);
end

 and so on ..

Thank you

No, that's not right:

Your polyphase components need to be ... well, phases of your signal, and hence, only take every $N$th component. What you do is something different.

To say "these are polyphase components of $h$", you need to define how many components you'll have. You're mixing different $N$s in your code, but overall, your code is just confused – it never selects the $N$th input samples.

Hence, you'd need something like

\begin{align}
p^4_0[n] &= h[4n]\\
p^4_1[n] &= h[4n+1]\\
p^4_2[n] &= h[4n+2]\\
p^4_3[n] &= h[4n+3]
\end{align}

for a $N=4$-phase polyphase decomposition of $h$ (of a specific type; "type" specifies in which order you take the $+0$,$+1$, $+2…$ components and declare them as $p_0$, $p_1$ and so on).

Feels like you need a better text book to start with! In English, the classic would be fred harris' Multirate Signal Processing for Communication Systems. It does an introduction that I'd explain differently these days, but I've worked with such systems a bit, and my initial understanding was guided by that book, so I consider it relatively good.


Criterion to choose length of an accumulator and product I am designing an IIR filter with fixed-point arithmetic and I have to select the proper length for the accumulator and product. I would like to know a standard to follow in order to choose its length. I have selected the Direct Form 1 topology to implement a 6th order digital filter (instead of a cascaded second-order block), as it has only only one point where the sum is done and thus only one point which has to be taken care of when considering quantification error and overflow (correct me if I am wrong!).

For the ACCUMULATOR, I assume that I have to consider the worst case of summing numbers with N bits, so the answer should be N+M bits to avoid overflow, where M represents the quantity of numbers which it is summed.

For the PRODUCT, multiplying two numbers of N bits, in the worst scenario gives us a N+N=2N bits result.

I don't know if it is this easy or if I am skipping an important concept.

If you have experience in this field, any documentation about designing IIR filters on microcontrollers will be welcomed.

Fixed point processing is very difficult. Floating point is A LOT easier.
The best algorithm and scaling approach depends a lot on your specific filter and the statistics of your signal. There is no "one size fits all" solution.
Cascaded second order sections are almost always the way to go. Primarily they guarantee that your coefficients values are bounded to something reasonable (2 or less). With transfer function implementation the coefficients tend to get really large or really small wich generates a lot of noise or scaling problems.
Within each section, use either Direct Form I or Transposed Form II. This helps controlling the scaling of the state variables and intermediate result
The best way to approach this to use "Q notation" https://en.wikipedia.org/wiki/Q_(number_format) to determine for EVERY SINGLE operation, which bits to use for what purpose. Doing this analysis will determine which bits to compute and which bits to keep (using multiplies, shifts, adds and masks).
Sometimes you need to treat this as a statistical problem. In many cases doing conservative scaling results in very poor signal to noise ratio, so you need to trade off the probability of clipping versus your steady state signal to noise. That's highly dependent on your signal, the specific filter and the requirements of your application.
You also need to carefully control your clipping and rounding behavior to avoid wrap-arounds and limit cycles.


Why Doesn't Compressive Sensing Work for Any Signal? My (probably naive) understanding of compressive sensing is that it is a technique that allows to efficiently reconstruct an $N$-dimensional signal $\boldsymbol x$, provided that it is sparse in some basis (without the need to know the sparsity basis), using a number $M\ll N$ of measurements.

More formally, if the signal can be written as $\boldsymbol x = \Psi\boldsymbol s$ with $\boldsymbol s$ $K$-sparse, and we are measuring in an incoherent $M$-dimensional measurement basis $\boldsymbol y = \Phi\boldsymbol x$ with $M\ll N$, then we can reconstruct the signal $\boldsymbol x$ using only $M$ measurements (provided $M\ge K$) (with the notation of Baraniuk 2007).

What I don't understand is why doesn't the above sparsity requirement apply to all possible signals.
I mean, isn't any signal sparse in its own basis?
Regardless of the choice of $\boldsymbol x$, I can always take as $\Psi$ any matrix whose first column is $\boldsymbol x$, and this will result in $\boldsymbol x$ being representable with the $1$-sparse vector $\boldsymbol s\equiv (1,0,0,...)^T$.

So why can't I efficiently reconstruct a signal which is $1$-sparse in a basis like the above?
More generally, what are the requirements on the basis in which the signal has to be sparse for a CS approach to work?

You can measure and reconstruct such a 1-sparse signal as you describe in your question. The crucial misunderstanding here is, as @MBaz points out, that you have to know the basis $\Psi$.

You can take measurements in a universal manner without knowing the dictionary $\Psi$ beforehand, if you choose the measurement matrix $\Phi$ as a sufficiently random matrix - could be a random Gaussian matrix. What universality gives you is the ability to reconstruct $x$ again with good accuracy almost no matter what $\Psi$ turns out to be. You can take such universal measurements at the encoder without knowing $\Psi$, but you still have to know some $\Psi$, that admits a sparse representation of $x$, (and $\Phi$) at the decoder.


Filtering Passband Signals using Complex Baseband Filtering So I'm given a pass band filter with specific transfer function $H_p(f)$, I want to implement this via baseband processing.

I already know how to take the input signal $u(t)$ and process it such that I get the I/Q components $u_I(t), u_Q(t)$.

I also have figured out that the following convolution relation holds for the filter output $y$'s I/Q components $$y_I=(u_I * h_I) - (u_Q*h_Q)$$ and $$y_Q = (u_I*h_Q) + (u_Q*h_I)$$ where $h_Q$ and $h_I$ are the filters I/Q components.

I have run into an issue where I can't figure out a simple way to determine $h_Q$ and $h_I$.

I know for the complex envelope of the filter I have $$h(t) = h_I(t) + j\cdot h_Q(t)$$

I can take the fourier transform of this to obtain $H(f) = H_I(f) + j\cdot H_Q(f)$.

So essentially my issue is with how to obtain $H_I$ and $H_Q$.

My book basically gives the answer but I must be missing something since I can't figure out the reasoning behind it. They state that $$H_I(f) = (H(f) + H^*(-f))/2$$ and $$j\cdot H_Q(f) = (H(f)-H^*(-f))/2$$.

I'm guessing its some property to do with the fourier that I've forgotten but I'm hoping someone can explain the reasoning for it. Thanks!

Since $h_I(t)$ is the real part of $h(t)$ you have

$$h_I(t)=\frac12[h(t)+h^*(t)]\tag{1}$$

where $h^*(t)$ is the complex conjugate of $h(t)$. For the imaginary part you have

$$h_Q(t)=\frac{1}{2j}[h(t)-h^*(t)]\tag{2}$$

Since the Fourier transform of $h^*(t)$ is $H^*(-f)$1, the Fourier transforms of (1) and (2) are

$$H_I(f)=\frac12[H(f)+H^*(-f)]\\
H_Q(f)=\frac{1}{2j}[H(f)-H^*(-f)]$$

1The fact that $H^*(-f)$ is the Fourier transform of $h^*(t)$ can be easily seen as follows:

With
$$H(f)=\int_{-\infty}^{\infty}h(t)e^{-j2\pi ft}dt$$
the Fourier transform of $h^*(t)$ is
$$\mathcal{F}\{h^*(t)\}=\int_{-\infty}^{\infty}h^*(t)e^{-j2\pi ft}dt=
\left[\int_{-\infty}^{\infty}h(t)e^{j2\pi ft}dt\right]^*=H^*(-f)$$


Division by zero in Levinson–Durbin recursion? I am doing a LPC analysis of a speech signal using the autocorrelation method. To calculate the LPCs I am using the Levinson–Durbin recursion. In my literature the error is initialized with the first element of the autocorrelation sequence. The algorithm later uses the error (from the previous pass) to calculate the reflection coefficient. However, the error is used as divisor in a division. Thus, in case the first element of the autocorrelation sequence is zero this leads to a problem in my C++ implementation of LPC analysis.

I have seen implementations setting all LPCs to zero in case the first element of the autocorrelation sequence is zero. Other implementations ensure that the first element of the autocorrelation is always at least 1.
Matlab gives me LPCs of NaN if the first element of the autocorrelation is zero in a call to its levinson function.

What is the "correct" way to cope with this divison-by-zero-problem? Or what is the most suitable for LPC analysis of a speech signal?

Matt, this is a concise snippet of code from long ago:

/**********************************************************************
*                                                                     *
* Levinson-Durbin Algorithm for LPC Coefficients                      *
*                                                                     *
*   input:                                                            *
*   n  = LPC order                                                    *
*   r -> r[0] to r[n]    autocorrelation values                       *
*                                                                     *
*   output:                                                           *
*   a -> a[0] to a[n]    LPC coefficients, a[0] = 1.0                 *
*   k -> k[0] to k[n]    reflection coefficients, k[0] = unused       *
*                                                                     *
*                                                                     *
**********************************************************************/

#define N   17

levinson_durbin(const int n, const double* r, double* a, double* k)
  {
  double a_temp[N], alpha, epsilon;         /* n <= N = constant  */
  int i, j;

  k[0]      = 0.0;                          /* unused */

  a[0]      = 1.0;
  a_temp[0] = 1.0;                          /* unnecessary but consistent */

  alpha     = r[0];

  for(i=1; i<=n; i++)
      {
      epsilon = r[i];                       /* epsilon = a[0]*r[i]; */
      for(j=1; j<i; j++)
          {
          epsilon += a[j]*r[i-j];
          }

      a[i] = k[i] = -epsilon/alpha;
      alpha = alpha*(1.0 - k[i]*k[i]);

      for(j=1; j<i; j++)
          {
          a_temp[j] = a[j] + k[i]*a[i-j];   /* update a[] array into temporary array */
          }
      for(j=1; j<i; j++)
          {
          a[j] = a_temp[j];                 /* update a[] array */
          }
      }
  }

now i see many places (well, one place, but many instances) where division by zero is possible.  any time k[i] = 1 or -1.  can't that conceivably happen with a non-zero input?


How is $\delta(at+b)=\frac{1}{|a|}\delta(t+b/a)$? This result has been used in the second last line of the pic. I don't know why it's true. Both functions are zero for $t$ not equal to $-b/a$. But at $t=-b/a$, a scaling factor $1/|a|$ has been introduced in the second function.

First of all it seems useful to establish what we mean by an equation like

$$\delta(at+b)=\frac{1}{|a|}\delta(t+b/a),\qquad a\neq 0\tag{1}$$

Since the Dirac impulse $\delta(t)$ is a distribution, Eq. $(1)$ only makes sense if interpreted as

$$\int_{-\infty}^{\infty}\delta(at+b)\phi(t)dt=\frac{1}{|a|}\int_{-\infty}^{\infty}\delta(t+b/a)\phi(t)dt\tag{2}$$

where $\phi(t)$ is a so-called test function, by which we generally mean that it has derivatives of any order, and that for $|t|\to\infty$ it tends to zero sufficiently rapidly.

In order to prove $(1)$, it is sufficient to show that

$$\delta(at)=\frac{1}{|a|}\delta(t),\qquad a\neq 0\tag{3}$$

because substituting $t+t_0$ for $t$ in $(3)$ results in

$$\delta(a(t+t_0))=\frac{1}{|a|}\delta(t+t_0),\qquad a\neq 0\tag{4}$$

which is equal to $(1)$ for $t_0=b/a$.

Proving $(3)$ means that we need to show that

$$\int_{-\infty}^{\infty}\delta(at)\phi(t)dt=\frac{1}{|a|}\int_{-\infty}^{\infty}\delta(t)\phi(t)dt,\qquad a\neq 0\tag{5}$$

Eq. $(5)$ can be shown in a straightforward manner by the variable substitution $at=\tau$, which results in $a\,dt=d\tau$:

$$\begin{align}\int_{-\infty}^{\infty}\delta(at)\phi(t)dt=\begin{cases}\displaystyle\frac{1}{a}\int_{-\infty}^{\infty}\delta(\tau)\phi(\tau/a)d\tau,&a>0\\\displaystyle\frac{1}{a}\int_{\infty}^{-\infty}\delta(\tau)\phi(\tau/a)d\tau,&a<0\end{cases}\tag{6}\end{align}$$

Since

$$\frac{1}{a}\int_{\infty}^{-\infty}\delta(\tau)\phi(\tau/a)d\tau=-\frac{1}{a}\int_{-\infty}^{\infty}\delta(\tau)\phi(\tau/a)d\tau\tag{7}$$

the result $(6)$ can be summarized as

$$\int_{-\infty}^{\infty}\delta(at)\phi(t)dt=\frac{1}{|a|}\int_{-\infty}^{\infty}\delta(t)\phi(t/a)dt,\qquad a\neq 0\tag{8}$$

which is equivalent to $(5)$ because

$$\begin{align}\frac{1}{|a|}\int_{-\infty}^{\infty}\delta(t)\phi(t/a)dt&=\frac{1}{|a|}\int_{-\infty}^{\infty}\delta(t)\phi(t)dt\\&=\frac{1}{|a|}\phi(0),\qquad a\neq 0\end{align}\tag{9}$$

Consequently, $(8)$ and $(9)$ prove $(3)$, which in turn proves $(4)$ and $(1)$.

Thanks to Laurent Duval for commenting and motivating me to rewrite this answer with a bit more (engineering) rigor.


Extracting metadata from wav files I was curious if there is a way to extract metadata from wav files. For instance, can you get bpm, key, instrument, etc. from a .wav file on upload?

As mentioned already, it is incredibly easy to read metadata off of WAV files as it is a tag based format and there are already implementations of readers for this information out there (e.g. tinytag in Python).

Having said this though, "bpm, key, instrument" are not part of typical information you would find in a WAV file. If the creator of the file put them there, you will be able to read it, but "bpm" is not essential for playing back an audio file. What you are more likely to be able to retrieve is things like Sampling Frequency, Bits per sample, number of channels and of course anything to do with any compression method used. Building an automatic "recogniser" which could be populating these fields given some WAV file, is possible but not exactly straightforward.

Bpm, key and instrument however are essential for playing back MIDI files. Perhaps these are more useful to you (?)

Hope this helps.


IQ FM Demodulation of Wow and Flutter Signal in Software Background:

I have a recording device that exhibits Wow and Flutter, or periodic frequency deviation. I use this device to record a calibrated 3.15kHz test tone from a tape film to generate a .wav file on my PC.

In an ideal world, an FFT on the data contained in the .wav file would result in a single sharp peak at 3.15kHz. In the real world, there is deviation in this frequency due to mechanical errors such as head scrape, non-concentric film rollers, motor vibration which all cause the tape film to either speed up or slow down slightly.

These deviations are called Wow and Flutter, and appear as errors frequency modulated by the 3.15kHz test tone.

What I would like to do:

Assume I have the .wav file, sampled at 48kHz. Let's say it contains an FM waveform with carrier frequency at 3.15kHz, and a message signal with bandwidth ~200Hz.

I'd like to build a python script to analyse and demodulate the signal to give me just the 200Hz signal. I've tried several models across the internet simulating PLLs, slope detectors etc. But these are all hard coded difficult to understand and I just can't seem to mould them to be used in my application.

I want to try implementing an IQ FM demodulator but I am unable to find a good reference to how FM waves are modulated via IQ signals and I definitely can't find one for demodulation.

Can someone point me in the right direction and let me know off the bat if there are any problems or roadblocks with what I'm trying to do?

Frequency is the derivative (or 1st difference) of phase.  To get an IQ signal from which you can get phase angles, first multiply your wave file samples by a cosine at 3.15kHz and a sinewave at a frequency of 3.15kHz.  That will hetrodyne your signal down to baseband IQ.  Use atan2 of the IQ array to get an array of  angles, take the 1st difference between successive angles, and note that big jumps in phase differences will likely need to be unwrapped to a smaller range.  Then low pass filter (with a 200 Hz cutoff) the result. That's a simplified FM demodulator .  See if that looks close to what you are looking for.


PSDs and Parseval's theorem OK, I have been banging my head for quite a while trying to make sense of this simple equation:

The $1/T$ is what is causing me all the grief.  I assume that the square of the Fourier transform of the signal $X$ is a finite function.  So won't the $T$ in the denominator mean that the PSD, $S_x$, will always tend toward zero?  If $X$ does increase to infinity with time, can someone explain this to me?  I am thinking that $X$ is a measure of the frequency content of the random time signal, and that shouldn't change appreciably with the time taken to measure it.

That equation is indeed simple, but it's also wrong. So your doubts are completely justified. What is probably meant is something like

$$S_X(f)=\lim_{T\rightarrow\infty}\frac{1}{T}\left|X_T(f)\right|^2\tag{1}$$

with

$$X_T(f)=\int_{-T/2}^{T/2}x(t)e^{-j2\pi f t}dt\tag{2}$$

This definition of the power spectrum is used for functions $x(t)$ with finite power (but infinite energy) for which the Fourier integral does not exist.


Unscented Kalman Filter - Multiple Consecutive Measurement Updates In trying to implement an Unscented Kalman Filter (UKF), I have come across the issue of what to do when my measurement signals come in at a different rate than my control inputs, which I use in the prediction step of the UKF. In my implementation, I've tried working around this issue by forming a prediction each time I receive a new measurement, but I've been told that it's also possible and perhaps favorable to only form a prediction when I get a new control input and form a new measurement update when I get a new measurement.

The "standard" Kalman filter algorithm in most textbooks/papers that I've seen do not delve deeply into this issue of different data rates, and they usually just show a simple predict-update loop that occurs at each time step.

I see how it is possible to form multiple consecutive predictions. However, I don't fully understand how one could do an update without a prediction right before it at the same time step, since the computation of the Kalman gain relies on the sigma points passed through the state transition function. Any help would be greatly appreciated!

For reference, I am trying to implement my UKF using the FilterPy Python library: https://filterpy.readthedocs.io/en/latest/_modules/filterpy/kalman/UKF.html#UnscentedKalmanFilter

I will use Wikipedia notations - Kalman Filter.
In most models the state transition model matrix $ F $ depends on the interval parameter $ T $. The same goes for the Process Noise Covarinace Matrix $ Q $. For instance, for the constant velocity model:
$$ F = \begin{bmatrix}
1 & T \\
0 & 1
\end{bmatrix}, \; Q = \begin{bmatrix}
\frac{ {T}^4 }{4} & \frac{ {T}^3 }{2} \\
\frac{ {T}^3 }{2} & {T}^{2}
\end{bmatrix} $$
Now, one simple policy you act by, whihc works perfectly in real life and is highly flexible, is as following:

Always remember the time stamp of the last Kalman step, be it prediction or update.
Given new input which is not a measurement do as following:

Take the input time stamp.
Define $ T = {t}_{current} - {t}_{last} $.
Build $ F $, $ Q $ and all other model matrices based on the given $ T $.
Apply prediction step.
Update the time stamp of the last Kalman step to $ {t}_{current} $.

Given new input which is a measurement do as following:

Take the input time stamp.
Define $ T = {t}_{current} - {t}_{last} $.
Build $ F $, $ Q $ and all other model matrices based on the given $ T $.
Apply prediction step.
Apply update step.
Update the time stamp of the last Kalman step to $ {t}_{current} $.

This works and it is very robust. Yet it has the overhead that you need to generate the correct matrices each step.
Also, since the UKF (Unscented Kalman filter) is different only in its Update step, it holds for it as well (For EKF as well).


2D Multi Scale Dot Enhancement Filter based on Gaussian Filter and Hessian Matrix I'm trying to implement an algorithm that enhance dot and line like structures based on Hessian matrix, the algorithm uses Gaussian filter with different scales before calculating the hessian matrix for every pixel and later the Eigen values.

I want to use the algorithm to enhance pulmonary nodules (dot like structures) to extract them for later classification.

The Gaussian filter is used to reduce noise and preserve objects with specific scales (diameters in my case).

Here is the algorithm that I'm trying to implement:

A = imread('sliceX.png');

scales = [3, 5, 7, 11, 17]; % the scales in number of pixels
% correspond to diameter of nodules in mm [1, 1.6, 2.4, 3.8, 6]
Ecircle = zeros(512, 512, length(scales)); % Ecircle will store the
% results of each enhancement scale.

for sn=1:length(scales) % sn for scale number
% Smooth the original 2D image with a 2D Gaussian function of scale
% Sigma_s.
    B = imgaussfilt(A, scales(sn));

    [gx, gy] = gradient(double(B));
    [gxx, gxy] = gradient(gx);
    [gxy, gyy] = gradient(gy); % it's normaly [gyx, gyy] but since
    % gyx = gxy it's ok

    % loop through the image to calculte the eigen values for every
    % pixel and based on that we choose the value of each scale
    %enhancement filter
    for x = 1:512
        for y = 1:512
            % construct 2*2 hessian matrix for every pixel
            h = [gxx(x, y), gxy(x, y); gxy(x, y), gyy(x, y)];

            e = eig(h); % returns a vector with the two eigen values
            % lambda1 and lambda2 in which lambda1 = e(1) and
            % lambda2 = e(2)
            lambda1 = e(1); lambda2 = e(2);

            if abs(lambda1) < abs(lambda2)
                temp = lambda1;
                lambda1 = lambda2;
                lambda2 = temp;
            end

            if lambda1 < 0 && lambda2 < 0
                Ecircle(x, y, sn) = (abs(lambda2)^2)/abs(lambda1);
            end
        end
    end
    % multiply each enhancement scale by (sigma^2) as mentioned in the
    % article
    Ecircle(:, :, sn) = Ecircle(:, :, sn) * (scales(sn)^2);
end

I = max(Ecircle, [], 3);

This implementation for enhancing dot like structures (nodules), and can enhance line structures (vessels in my case) with just changing the if condition.

The problem is I'm not getting the supposed results that I should get:

The result of the developer of the method in the paper:

And this is my result:

abvisioly the nodule in the image is greatly enhanced by there are great noise in the output image plus all the vessels joints are enhanced too (small and big)

I think the problem is with the Gaussian smoothing filter imgaussfilt(), and precisely with sigma or scales(sn) in this example.
I've read that sigma should be in the same units of x and y i.e. number of pixels so I've transformed the diameters (that I got from the experimental results of the original article) to the number of pixels using PixelSpacing attribute form the original Dicom file metadata.

The diameters are [1, 1.6, 2.4, 3.8, 6] to cover approximately all the diameters of possible nodules.

Where I did get wrong and If the problem is with the value of Sigma,How can I specify sigma? and how can I use the imgaussfilt() function correctly?

Note: here is the image used in the example:
input image

Questioner's answer...

sigma have the same units as x and y i.e. number of pixels.
In multi-scale filtering, the size of the filter must change when the sigma changes.
Obtain the number of pixels per one millimeter or the vice-versa. (I did this using the property of pixel spacing included in the DICOM metadata in Matlab you can do this as info=dicominfo('image.dcm'); and spacing=info.PixelSpacing;)

What I was doing wrong is not changing the filter size when I changed sigma in the different scales.

So, the solution to my problem is like this:
the size of the gaussian kernel should be 2 times or preferably 3 times the value of sigma on eather side of the origin as Mr. Cris suggested.
Gaussian filtering with Matlab's Image Processing Toolbox

cutoff = ceil(3*sigma);

In matlab there is two options, using fspecial() or imgaussfilt(), fspecial() is not recommended anymore in the newer versions of Matlab and latter is the recommended one.

h = fspecial('gaussian',2*cutoff+1,sigma);
B = conv2(A,h,'same');

or...

B = imgaussfilt(A, sigma, 'FilterSize', 2*cutoff+1);

Another thing that I did that corrected my results is converting my input image into type double using im2double() function.

Here is the result:


Impulse response of a system I have the following difference equation:
$y[n]=y[n-2]-0.25y[n-4]+x[n]-x[n-2]$
I need to find the first 6 values of the impulse response of the system characterized by that equation. When I do that, the values I get are:
$[1,0,0,0,-0.25,0]$, starting from index zero. I have the answer, which states that the right response is $[1,0,0,0,-0.25,1]$. I really cannot find my mistake; any help will be appreciated.

Answers in textbooks are not necessarily right.. I have got the same answer as you ;)

If you are not sure that your method is right, you can always test digital systems on your PC (matlab, scilab, python, octave, or almost in any language )


what is difference between structure image and texture image I am working in Computer Vision. One thing is confusing me, the difference between structure image and texture image. Both are looking same. but i read few papers, which showing that both are different. Someone tell me how they are different, and how can we said either this is structural or textural, or either both.

There is a growing awareness that images are not homogeneous and that it could be useful to model them as a combinaison or superposition of morphological components that have distinct properties. One can cite four main components: piecewise-smooth parts, contours, geometrical/periodic textures and a last one, barely modelled, sometimes termed noise.

The edge-cartoon/texture is a common name given to a more macroscopic version with two components: structure and texture. Structure often refers to the piecewise-smooth parts and large-scale contours, texture to geometrical and almost periodic textures, noise. Sometimes, they are called geometry and oscillatory components. Related works are mentioned in Structure-Texture Image Decomposition: See for instance: Image decomposition toolbox

When the contours of an object are not distinct, its textural aspect can help to segment it. The color aspect is often put on the "texture" side.


How channel state information is calculated from Sounding Packet What I understand is for Explicit Transmit Beamforming scenario, the Transmitter sends an NDP packet(aka HT Sounding) which consist of OFDM training symbols in packet preamble (HT-LTF) which are randomly generated +/-1 for each subcarrier(say 56).

I do not have much knowledge on this, but if anyone can help me to understand how these training symbols are calculated at the receiver end to evaluate CSI will be helpful.

According to pages 5-7 in this, the Channel State Information (CSI) contains 2 measurements per matrix element.

The CSI essentially tells the Access Point (AP) two things: where is the client (bearing) and how far away on that bearing it is. However, because of multipath propagation, the "bearing" might not correspond to the true bearing of the client but the direction at which the signal is "louder". This means that the signal could be bouncing on a wall and then finding the client, rather than traveling to it in a straight line.

So, the task we have now is to calculate these two numbers. We assume that the channel identification pseudorandom sequences in the pre-amble are transmitted at the same time from the AP:

To estimate the $A$, simply sum the "strength" of the received signal at each antenna.
To estimate the $T$, you need to estimate the time differences between the times the signal was received. At this point, you need to decide on the ordering of the elements.

You can estimate both using a matched filter whose $h$ is the pseudorandom sequence.

The matched filter will return a relatively high output when it detects the pseudorandom sequence at its input.

So, imagine a bank of matched filters behind each receiver. Each one of the antennas receives the AP signal at slightly different times and amplitudes. Each one of the matched filters produces a slightly different amplitude (the $A$ component) when it "sees" the pseudorandom sequence at its input. Consequently, because of the difference in arrival times, each one of these high outputs will occur at slightly different times. To estimate $T$, take one of your antennas as a reference element and measure the time that passes between two high outputs at each matched filter.

You might also find reading up on the fundamental of the following techniques:

The Matched Filter and its applications in Radar and Sonar
Beamforming

Which sooner or later will lead you to Adaptive Beamforming which is what you are dealing with here.

Phased Arrays

Which, again, is what you are dealing with here.

Pseudorandom binary sequences (PRBS)

Which is what your header sequence is right here

PRB Sequences with applications to channel characterisation

In brief, by examining the distortions that a PRBS undergoes when transmitted through a "channel" (whatever this might be), you get an estimate of the channel's condition or state which you can then use to adapt some aspect of transmission (Forward error control sequence, beamforming, etc).
Characterisation of an acoustic communication channel using PRBS
Transmitter identification using PRBS
And more...

Hope this helps.


Refer good books for the following topics in Image processing Thanks for reading my question. Please refer good books(links, if possible) for the following topics in Image processing. Your help is very much appreciated.

UNIT II SEGMENTATION

Fuzzy clustering, Watershed algorithm, Active contour methods, Texture feature based segmentation, Model based segmentation, Atlas based segmentation, Wavelet based Segmentation methods.

UNIT III FEATURE EXTRACTION

Texture descriptors- Autocorrelation, Co-occurrence features, Runlength features, Fractal model based features, Gabor filter, wavelet features.

UNIT IV REGISTRATION AND IMAGE FUSION

Image Fusion-Overview of image fusion, pixel fusion, Multiresolution based fusion discrete wavelet transform, Curvelet transform. Region based fusion.

for Multiresolution based fusion discrete wavelet transform you may refer the following
" A Wavelet based image fusion tutorial " by Gonzalo Pajares ( available at: www.ELSEVIERcomputer science.com)

"Wavelet based image fusion techniques - An introduction ,review and comparison by Krista Amolins

for Curvelet transform refer " Wavelet for image fusion" by shih-Gu Huang


Understanding GCC-PHAT as a Feature I'm working on a sound source localization project and am interested in using GCC-PHAT as a feature. I'm quite new to this transformation and I've been reading a few papers about it and I'm still struggling with understanding what it actually is. One of the papers I'm interested in contains the following:

The maximum possible time delay between any 2 microphones is
$τ = 0.2/340 = 0.5882ms$ where the sound speed is assume to be
$340m/s$. Suppose we are using a sampling rate of 16kHz, then the
maximum delay in samples is $n = 16000τ ≈ 21$. Hence, for each
microphone pair, only the 21 correlation coefficients near the center
contains useful information for DOA estimation.

First of all I'm not sure where $n = 21$ is coming from, seing as $0.2/340 * 16000 = 9.411$, but perhaps its simply an error on their part, but I was curious if this has something to do with GCC-PHAT. I've found several implementations and when I take a 10ms audio sample with a sample rate of 44.1kHz (441 samples) I get the following plot:

So this is where I start to have questions as to what the actual feature is. But more to the point about the quote above, what exactly are the coefficients that are being referenced? Are these simply the correlation values around the center? And by transforming my x-axis to milliseconds, it seems that the peak is around 6ms. So does this mean that the signal arrived at microphone-1 6ms after microphone 0? Also any literature on the topic would be greatly appreciated, as I've so far only really been able to find code implementations but I'm trying to get a more intuitive understanding of the thing itself. Any help is greatly appreciated.

Let there be two microphones, $x_1(t), x_2(t)$ capturing a source, $s(t)$, with some additive noise $n_{1,2}(t)$ that is uncorrelated with the source and with each other (i.e., its effect can be ignored while calculating the cross-correlation). Let $\tau_0$ (in seconds) be the time difference of arrival between the two microphones. The time domain equations of the microphones are
$$x_1(t) = s(t) + n_1(t) \\
        x_2(t) = \alpha_0 s(t - \tau_0) + n_2(t) $$
Now, taking the FFT and converting these to the frequency domain would give us,
$$
        X_1(\omega) = S(\omega) + N_1(\omega) \\
        X_2(\omega) = \alpha_0 S(\omega) \exp{(-j\omega \tau_0)} + N_2(t)
$$
The cross correlation function of the two microphone signals (or the cross-power spectrum in the frequency domain), can be written as
$$      \begin{aligned} R_{x_1,x_2}(l) &= \mathbb{E}[x_1(t) x_2(t-l)] \\
        R_{x_1, x_2}(\omega) &= X_1(\omega)X_2^*(\omega) \\
                             &= |S(\omega)|^2 \alpha_0 \exp{(j\omega \tau_0)}
        \end{aligned}
$$
The Generalized Cross-Correlation function with Phase Transform (GCC-PHAT) at lag $l$ is defined as
$$ \begin{aligned}
\tilde{R}_{x_1,x_2}(l) =& \frac{1}{2\pi} \int_{-\pi}^{\pi} \frac{R_{x_1,x_2}(\omega)}{|R_{x_1,x_2}(\omega)|}e^{j\omega l}d\omega \\
=&  \frac{1}{2\pi} \int_{-\pi}^{\pi} e^{j\omega(l + \tau_0)} d\omega \\
&= \delta(l - \tau_0)
\end{aligned} $$
where $\delta$ is the Dirac-delta function which is non-zero at $l = \tau_0$, which is the time-difference of arrival (TDOA). Typically, the lag is calculated in samples, so to convert that to seconds, you have to divide it by the sampling rate.
I hope this helps in understanding GCC-PHAT, which is useful tool for TDOA estimation, and hence, source localization.


How the noise in time domain affect the power in frequency domain after FFT? I am new to signal processing. And I am puzzled when studying FFT.
I create time series plus a random noise in Python:

import numpy as np
N=1000
x = np.linspace(-50,50,N)
noise = np.random.random(N)*0.1
y = np.sin(2*np.pi*x/10) + noise

Then perform FFT to the time series and get power spectrum.
But I find the spectrum is totally different when I execute my script each time.
I know it must be the random noise's effect, which are different in different execution. However, as you can see, the sine function is periodic, so we should see clearly a signal in the power spectrum. The puzzle is, the very strong periodic "signal" is submerged in the noise in frequency domain FOR SOMETIME.

I want to know how the noise in time domain affect the noise&signal in frequency domain.

Many thanks!

You have problems in defining the noise and sampling frequency. Your noise is not zero mean and your sampling rate is not well defined. Following simple Python code computes the DFT magnitude of the sine plus noise example. (Note: I'm not a python user so please modify the code as necessary)

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft
import numpy.random as npr

N = 1000                # Set signal sample length

t1 = -5;                # Simulation begins at t1
t2 =  5;                # Simulation  ends  at t2
x = np.linspace(t1,t2,N);# Define the sampling grid as x[n]

Ts = (t2-t1)/N;         # Compute resulting sampling period
Fs = 1/Ts;              # Compute resulting sampling rate

noise = (npr.random(N)-0.5); # Generate ZERO MEAN uniform noise

fsin = 0.2*Fs;            # Sine Wave Frequency in [0:Fs/2]
y = np.sin(2*np.pi*fsin*x) + noise    # Generate SINE + Noise

# Take FFT and display its MAGNITUDE
w = np.linspace(0,2,N)
waves = fft(y,N)
plt.plot(w,abs(waves))

Where the resulting plot is:


how to estimate the noise of a signal? please can someone help me to define a method that allows me to estimate the power of noise in any signal?
the signal is received by usrp and the problem that it is variable all the time in the FM band so the noise is variable also with frequency and time

thanks in advance for any help

Noise floor estimation is usually done after applying an FFT to windowed data segments. By tracking the noise floor in each frequency band, the frequency dependence of the noise is taken into account. If the noise is non-stationary, its time dependence can be tracked by regularly updating the noise floor estimate in each frequency band. This usually requires the noise to vary more slowly than the desired signal. Have a look at this paper to get an idea of how it can work, and also look at the references therein. Search the web for "noise floor estimation" and "noise floor tracking" to find many more papers on the subject.


length of window and overlap rate in STFT I want to use STFT to analyze my signal and am wondering what are differences between two solutions:

Use short windows (for ex. 256 samples window)
Use longer windows (to get higher resolution in freq. domain) but let the windows overlap more (For ex. I use 1024 samples window, but let 1023 samples overlap)

Thank you :)

In general, shorter windows contain less information, and will thus provide less frequency resolution.  But short windows will be more localized in time, and thus a sequence of short windows might allow greater resolution in locating or separating time domain behavior.

Longer windows will contain more information, but this information will be averaged across the entire window duration, thus potentially smearing out any transients or other non-stationary behavior that might change more rapidly than the window length or duration.  But the greater information content might allow a higher resolution (both for separating peaks and locating or estimating isolated stationary frequency peaks).

Overlapped windows will provide some partially redundant information, but perhaps better localized data regarding rate of phase change for some types of signals that are stationary across multiple frames.

Added:  Overlapping is also common when a non-rectangular window is used.  Since the windowing process (combined with finite arithmetic) can be quite lossy at the edges, overlapping helps provide some of the information that was thrown away by the windowing.


Why are 8 points needed to compute a fundamental matrix? Homograpies have 8 degrees of freedom.
Fundamental matrices have 7.
Essential matrices have 5.

To the best I could find, all those matrices are obtained by using the 8 point algorithm.

I don't see why the fundamental and essential matrices can't make do with less points.

Please enlighten me

Hum, no

A homography can be exactly fit to 4 point such that no three of them are collinear (example implementation in OpenCv).
An essential matrix can be fit to the image of 5 non-coplanar points (implementation).
A fundamental matrix can be fit to 7 points (implementation)


At what amplitude does spectral leakage start to be noticeable in the phase spectrum? I understand that discontinuities at the boundaries of a finite signal may introduce spectral leakage, that alter both the amplitude and phase spectra of the FFT, specially at low amplitudes. That is the explanation to several questions like (baffled by fft phase spectrum!"). In such questions, the examples are normally based on sine waves, where there is a dominant frequency, and the spectral leakage is noticeable in all the other frequencies.

Signals are often more complex, with a more continuous amplitude spectrum, but even in that case - specially if they are low pass-filtered, there is a point where there is a noticeable discontinuity in the amplitude and phase diagrams.

Example (in Octave):

N = 101;
x = exp(-linspace(-5,5,N).^2);
f = fftshift(fft(x));
subplot(211)
semilogy(abs(f))
subplot(212)
plot(unwrap(angle(f)))

Here the amplitude spectrum is bell-shaped as I expected at low frequencies. But when the amplitude goes below approximately 1e-11, the amplitude spectrum flattens.

At the frequencies where the amplitude spectrum is bell-shaped, the phase spectrum is a straight line with slope $\pi (N-1)/N$, which coincides with the effect of a time shift of $(N-1)/2$ from an even signal. At higher frequencies, it goes more or less linearly towards $\pi/2$ (although in this case the slope is not strictly constant).

Why there is such a harsh change in the spectra at that frequency? Is it a matter of numeric error in the computation of the FFT? The amplitude at the discontinuity is small (about 1e-11), but greater than machine error. And for other signals with narrower spectra, the amplitude at the discontinuity is higher. For instance, the same example but using a "wider" time signal:

x = exp(-0.5*linspace(-5,5*(N-1)/N).^2);

The issue is that you are not taking a large enough piece of your original x to get an accurate frequency domain representation.

Try setting

x = exp(-linspace(-20,20,N).^2);

If I do that, I get the bottom two plots below (the middle plots are for -10,10).  The top plot is a replication of your values.

In the three pairs below, the left plot is the semilogy plot of x and the right plot is the semilogy plot of its FFT f.

Derivation of the FFT of a windowed Gaussian

To answer your question in the comments, we need to look at how a windowed version of:
$$
x(t) = e^{-t^2/2}
$$
can be analyzed in the frequency domain.

The windowed version is:
$$
x_w(t) = w(t)e^{-t^2/2}
$$
where $w(t) = \mbox{rect(t/T)}$ if we use a rectangular window of width T (span from $-T/2$ to $+T/2$).

The Fourier transform will be:
$$
X_w(\omega) = \int_{-T/2}^{T/2} e^{-t^2/2}  e^{j\omega t} dt\\
= 2 \int_{0}^{T/2} e^{-t^2/2} \cos(\omega t) dt\\
= \left . \sqrt{\frac{\pi}{2}} e^{-\omega^2/2} \left[
\mbox{erf}(\frac{t - j\omega}{\sqrt{2}}) + j \mbox{erfi}(\frac{ \omega - j t}{\sqrt{2}} )
 \right] \right |^{T/2}_0\\
= \left . \sqrt{\frac{\pi}{2}} e^{-\omega^2/2} \left[
\mbox{erf}(\frac{t - j\omega}{\sqrt{2}}) + \mbox{erf}(\frac{ t + j\omega}{\sqrt{2}} )
 \right] \right |^{T/2}_0
$$
where this integral is used and the identity
$$
\mbox{erfi}(\frac{ \omega - j t}{\sqrt{2}} ) = \mbox{erf}(j \frac{ \omega - j t}{\sqrt{2}} ) / j = -j \mbox{erf}(\frac{ t + j\omega}{\sqrt{2}})
$$

Unfortunately, I'm not sure where that gets us. I'll leave it for now and see if I can doodle more later.

scilab code below (Sorry, I tend to use scilab rather than Octave).

// Q26236;
N = 101;
x1 = exp(-linspace(-5,5,N).^2);
f1 = fftshift(fft(x1));
x2 = exp(-linspace(-10,10,N).^2);
f2 = fftshift(fft(x2));
x3 = exp(-linspace(-20,20,N).^2);
f3 = fftshift(fft(x3));

clf
subplot(321)
plot2d("nl",abs(x1))
title('x1')
subplot(323)
plot2d("nl",abs(x2))
title('x2')
subplot(325)
plot2d("nl",abs(x3))
title('x3')

subplot(322)
plot2d("nl",abs(f1))
subplot(324)
plot2d("nl",abs(f2))
subplot(326)
plot2d("nl",abs(f3))


Does Circular Convolution Correspond to Periodically Expanded Linear Convolution? I have a question whether circular convolution and periodically expanded linear convolution corresponds in following case or why it does not?

Think about a signal $x[t]$ and a signal $y[t]$ both of length $N$.
We expand the singal $x[t]$ periodically to a signal $x'[t]=[x[t],x[t]]$ with length $2N$. We zeropadd the signal $y[t]$ to a signal $y'[t] = [y[t],0]$ to obtain the same length $2N$.

My question is does the circular convolution result $(x[t] \star y[t] )_{mod_N}$, correspond to the bins $[N...2N]$ in the linear convolution result $x'[t] * y'[t]$?

Best

Yes, they are the same.

Let us take the linear convolution of $x[t]$ and $y[t]$ as 2 portions $H_1$ of length $N$ corresponding to first $N$ samples, and $H_2$ of remaining $N-1$ samples. Circular convolution between $x$ and $y$ causes $H_2$ to overlap over $H_1$ because of time-aliasing. So the first $N-1$ samples of the result is $H_1 + H_2$ with only the last sample being correct (corresponding to linear convolution).

Now, for the linear convolution between $\tilde{x}$ and $\tilde{y}$, the samples from $N$ to $2N-1$ will have additional component of first $N$ samples of $\textrm{conv}(x,y)$. So these samples will correspond circular convolution of $x$ and $y$.


How can I generate offset QPSK symbols from QPSK symbols I am trying to generate Offset-QPSK modulated symbols. But I am facing some difficulties. I know how to make QPSK symbols but I am unable to figure out at what stage I have to bring in changes to make it OFFSET QPSK.

The following are the steps I did:

Generated a stream of 1s and 0 s
I have paired at 2 consecutive bits and mapped it to one among the integers 0,1,2,3.
Now I map these integers to $(1+j), \ (-1+j), \ (-1-j), \ (1-j) $ following grey code based ordering.

I know that these are the steps for QPSK. But I am not sure If I need to generate OQPSK symbols at which step should I make changes and what are the changes?

Any help is appreciated.

The idea of OQPSK is very simple:

Consider your QPSK as two BPSK, one orthogonal to the other. It helps remembering that you're free to rotate the constellation, so that the constellation points are on the I and Q axis. Instead of choosing one point of a four-points constellation based on two bits, you now choose two separate constellations points from two two-point constellations.
pulse shape I and Q independently with a pulse shaper that has more than 2 Samples per Symbol. Typically, use 5 or more.
offset (in time) either component so that it lags half a symbol time behind the other

If you already have QPSK-modulated baseband, you can still do step 3.


Multipath resistance of spread spectrum systems I have been reading about various telecommunications technologies and one thing keeps being repeated many times without a clear explanation - that spread spectrum is not only resistant to multipath interference, but that it can even use it to improve recovery of the data being transmitted. This is frequently mentioned in articles about LTE or DVB-T.

How does this work? I can imagine that something like autocorrelation could be used, but I'm also pretty sure that this is not the way it's done.

The key idea in spread-spectrum communications is to use much more bandwidth than the minimum required for the signal being sent.

Multipath effects occur when the sent signal $s(t)$ and a delayed, attenuated version of it interfere at the receiver:
$$
r(t) = s(t) + \alpha \cdot s(t-\tau)
$$

If $s(t)$ is really narrow band (a sinusoid), and $\alpha = 1$, then
$$
r(t) = \sin(\omega t) + \sin(\omega(t+\tau))
$$
and if $\omega\tau = \pi$ you can get perfect nulling:
$$
r(t) = \sin(\omega t) + \sin(\omega t + \pi) = 0
$$

If, however, $s(t)$ is not a simple sinusoid, then perfect nulling can be eliminated:
$$
r(t) = [\sin(\omega t) + \sin(\omega_2 t) ] + [\sin(\omega(t+\tau)) + \sin(\omega_2 (t+\tau)) ]
$$
even if  $\omega\tau = \pi$, if $\omega_2$ is chosen well, then $\omega_2\tau \ne \pi$ so there is still some signal to receive.


What does "A Target with Zero Doppler" Mean? I have seen several references to target at zero Doppler and some simulations that command users to place a target at zero Doppler. I know that when a target is flying tangential to the signal the Doppler is zero (is this even true?) but I have no clues as to what the target being at zero doppler means.

In Pulse Doppler Radar, the Doppler frequency that a target's return signal exhibits can be used to calculate the Radial Velocity of target.  That means, the return frequency of a pulse will be different than the frequency of the pulse that was sent out depending on how fast the target is moving toward the source of the pulse.  If the target's "Radial Velocity" is zero, zero doppler will result, i.e. the resulting doppler caused by the target's Radial Velocity will be 0 Hz.

Wikipedia shows the calculation of Doppler to be:

$$
F_{Doppler} = \frac{2 * F_{Transmit} * V_{Range} }{c}
$$

$c$ being the Speed of light.

So, if the target's Radial Velocity is zero (not moving towards or away from the observer), $F_{Doppler}$ will also be zero.

https://en.wikipedia.org/wiki/Pulse-Doppler_radar


"Perfect waveform" for cross-correlations? I am well versed in cross-correlations, as my masters thesis heavily relied on them for music classification and beat detection. So this question focuses more on the signals that I am running an xcorr on as opposed to the actual cross correlation process.

Here is a simplified version of my setup:
Suppose I have the original signal, which I time shift the signal (add a delay), and add noise/distort the signal. Then I run a cross-correlation on the two signals to find this time delay and get some value of R.

Is there some formulation to come up with a signal that will reduce ambiguity errors and optimize my value of R?
I know it will be some non-repeating, pseudo-random signal, but is there a formal formulation for this? Or should I be looking at other techniques?

Perhaps something like this would be perfect: http://tedxtalks.ted.com/video/TEDxMIAMI-Scott-Rickard-The-Wor. You see any issues with this?

I'm not talking about cross correlation processing techniques to improve the performance. I'm strictly talking about the shape of the waveform

In general, the resolution of the cross correlation is proportional to the bandwidth of the signal. This is way pseudo-random signals are common for that purpose.
So if you are looking for a signal that will produce sharp and distinct cross-correlation, you should use such signals.

However, if you are looking to improve the time-delay estimation for a given signal, you should pre-filter the signal with high-pass filter as the lower frequencies tend to produce more smooth cross correlation.


Is there any practical application for performing a double Fourier transform? ...or an inverse Fourier transform on a time-domain input? In mathematics you can take the double derivative, or double integral of a function. There are many cases where performing a double derivative models a practical real-world situation, like finding the acceleration of an object.

Since the Fourier transform takes a real or complex signal as an input, and produces a complex signal as an output, there is nothing stopping you from taking that output and applying the Fourier transform a second time... Are there any practical uses for doing this? Does it help to model some complex real-world situations?

With the same logic, nothing would stop you from taking the inverse Fourier transform of your original time-domain input signal... would this ever be useful? Why or why not?

No, taking the Fourier transform twice is equivalent to time inversion (or inversion of whatever dimension you're in). You just get $x(-t)$ times a constant which depends on the type of scaling you use for the Fourier transform.

The inverse Fourier transform applied to a time domain signal just gives the spectrum with frequency inversion. Have a look at this answer for more details.


How to amplify regions with specific color in an image? I want to amplify regions with a specified color in an image. For this purpose I have defined a simple function

$B(x,y)=s_r(x,y) \cdot w_r + s_g(x,y) \cdot w_g + s_b(x,y) \cdot w_b$

which produces a gray scale image B(x,y).

The weights I have defined with

$w_r \in [-2,2],w_g \in [-2,2], w_b \in [-2,2]$

The channels of the images are the functions $s_x$ for each channel.
I use the result of this amplification for a color thresholding. For example for yellow I have manually defined
$w_r=1$, $w_g=1$ and $w_b = -2$.
My Problem now is, to get "good" values for the weights automatically. First I thought taking the maximal values of histograms would be a good way, but the results are not satisfactorily.

How can I automatically determine the weights for the amplification?

How would you deal with white and black?

Or is there a better way to amplify the regions with a specific color?

Let us assume you are doing a binary classification with the desired color region set as 1 after applying your algorithm.
Then you can provide a set of color images (histograms) and desired outputs. The weights can be posed as the unknown quantities for this set of linear equations. Lo you have a overdetermined system & you can apply any method like Moore-penrose inversion to estimate the weights.


What is the algorithm to do a Discrete Hilbert Transform? I have a time series data. I have applied Empirical Mode Decomposition and obtained several IMFs.

Now I need to apply a Hilbert Transform to the IMFs to get the instantaneous frequencies.

I have found two papers explaining how to do that:

this one, that says this:

What Dr. Huang proposed for more accurately calculate the instantaneous frequency is to do normalization of IMF. The process is as followed:

Take absolute value of IMF.
Find extrema.
Based on these extrema, construct envelope.
Normalize IMF using the envelope. The FM part of signal becomes almost equal amplitude.
Repeat process 2-4 after the amplitude of normalized IMF retains a straight line with identical value.
Find the instantaneous frequency on the normalized IMF.

and this one that says this:

Notice that the first one says to take the absolute value of IMF and the other one doesn't.

I have tried both but repeating the loop to make all values be contained between 0 and 1 is doing nothing. Every loop the values are the same.

How do I really do that and how do I find the frequency of each IMF?

The Discrete Hilbert Transform is an "ideal" (implying that this is not the practical implementation) linear time-invariant filter with input $x[n]$ having output

$$ \mathscr{H}\Big\{ x[n] \Big\} \triangleq \hat{x}[n] = \sum\limits_{i=-\infty}^{\infty} h[i] \, x[n-i] $$

Where the impulse response of a discrete-time Hilbert transformer is:

$$ h[n] = \begin{cases}
 \frac{\big(1 - (-1)^n\big)}{\pi n} \quad & n \ne 0 \\
\\
0 & n = 0
\end{cases}$$

Because $1 - (-1)^n = 0$ for even $n$, this can be restated as

$$ h[n] = \begin{cases}
 \frac{2}{\pi n} \quad & n \text{ odd} \\
\\
0 &  n \text{ even}
\end{cases}$$

This is not a causal impulse response, nor is it finite in length.  To make it finite in length, you would need to window it with a decent window:

$$ h[n] = \begin{cases}
 \frac{2}{\pi n} w[n] \quad & n \text{ odd} \\
\\
0 &  n \text{ even}
\end{cases}$$

where $w[n]$ is some window function of width $L+1$ samples (and $L$ is an even positive integer).  If it were a Hamming Window, it would be:

$$ w[n] = \begin{cases}
  0.54 \ + \ 0.46 \cdot \cos\left(\pi \frac{n}{L/2} \right) \qquad & |n| \le L/2 \\
  0 & |n| > L/2 \\
\end{cases}$$

If it were a Kaiser window it would be

$$ w[n] = \begin{cases}
  \frac{1}{I_0(\beta)} \, I_0\left(\beta \sqrt{1 -  \left(\frac{n}{L/2}\right)^2 } \right) \qquad & |n| \le L/2 \\
  0 & |n| > L/2 \\
\end{cases}$$

where

$$ I_0(u) \triangleq \sum\limits_{k=0}^{\infty} \frac{(-1)^k \big( \tfrac{u}{2} \big)^{2k}}{(k!)^2} $$

$I_0(x)$ is the 0th-order modified Bessel function of the 1st kind. $L+1$ is the number of non-zero samples or FIR taps (the FIR filter order is $L$) and, in my centered and symmetrical case, must be even. $\beta$ is a "shape parameter", maybe around 5 or 6, i dunno.

Now, to make this causal, your impulse response has to be delayed by $\frac{L}2$ samples to be $h[n-\frac{L}2]$, but then should also all other signals that this Hilbert output is compared to, they should also be delayed by $\frac{L}2$ samples to keep the phase relationship correct.

With a windowed finite-length impulse response (which is what we call an "FIR"), the Hilbert transform output is

$$ \hat{x}[n] = \sum\limits_{i=-L/2}^{L/2} h[i] \, x[n-i] $$

and delaying the output so that the filter is causal you get

$$ \hat{x}[n-\tfrac{L}2] = \sum\limits_{i=0}^{L} h[i-\tfrac{L}2] \, x[n-i] $$

but with this delayed output $\hat{x}[n-\tfrac{L}2]$, you must compare that only to the like delayed input $x[n-\tfrac{L}2]$ in order for the two signals to have their 90° phase relationship (which is fundamentally what the Hilbert Transform is about).  Note that every even-indexed sample of $h[n]$ is zero, so the number of taps is not really $L+1$ but is $\frac{L}2 - 1$ taps with non-zero coefficients.


Time scaled version of discrete time signals I have a confusion regarding time scaling of discrete time signals.
Let u[n] be the discrete time unit step signal(or sequence)
Now if i plot u[n/3] i know there will be two zeros added between each previous samples(samples before scaling).
But when i tried plotting the same on MATLAB i got the same plots for u[n] and u[n/3] and even for u[4n].
How could this be?
I am attaching the MATLAB code :
close all;
clear vars;
t=-4:1:20;
stepd=t>=0;
subplot(2,2,1)
stem(t,stepd,'b');
xlabel('n');
ylabel('u[n]');
title('discrete unit step');
stepd1=1/3.*t>=0;
subplot(2,2,2)
stem(t,stepd1,'g');
xlabel('n');
ylabel('u[n/3]');
title('time scaled discrete unit step');
stepd2=4.*t>=0;
subplot(2,2,3)
stem(t,stepd2,'r');
xlabel('n');
ylabel('u[4n]');

I know this is a very silly doubt but please help me come out of this confusion.
I have also attached the plots which i am getting on MATLAB

Now if i plot u[n/3] i know there will be two zeros added between each previous samples(samples before scaling).

No.
You are confusing up-sampling with time scaling. The time scaled version of a unit step is a unit step. You would only insert 2 zeros if you up-sample by a factor of 3.


Recursive, non-recursive systems; FIR, IIR systems I am confused with the classifications of an LTI system as recursive or non-recursive systems and FIR or IIR systems.
I understood what the FIR and an IIR systems are, but is it correct to say that FIR system is always non-recursive?

We could express an finite accumulator up to past N inputs (FIR system) in both non-recursive and recursive forms.
Also is it correct to say that a non-recursive system is always an IIR system or vice versa?

The logical implications are the following:
"non-recursive" $\Longrightarrow$ FIR
IIR $\Longrightarrow$ "recursive"
But the opposites are not necessarily true because a FIR system can be implemented recursively (transfer function poles can be cancelled by zeros).

Of course, when referring to "recursive" or "non-recursive" we always talk about implementations with finitely many operations per output sample. Clearly, any discrete-time LTI system can be described by a generally infinite convolution sum, but that is not what we mean by "non-recursive".


Reconstruction of audio signal from Spectrogram: Implementation Is there somewhere an implementation available for the reconstruction of audio from a spectrogram? (e.g. based on this approach).

I could not find anything, but I think there should be  an implementation somewhere. Where can I find it?

Thank you

Here is an example in Matlab that I found some years ago. It has been turned to Python by one of us, apparently (STFT ISTFT Matlab Python).

function d = stft(x, f, w, h)
% D = stft(X, F, W, H)                            Short-time Fourier transform.
%   Returns some frames of short-term Fourier transform of x.  Each
%   column of the result is one F-point fft; each successive frame is
%   offset by H points until X is exhausted.  Data is hamm-windowed
%   at W pts..
%   See also 'istft.m'.
% dpwe 1994may05.  Uses built-in 'fft'
% $Header: /homes/dpwe/public_html/resources/matlab/RCS/stft.m,v 1.1 2002/02/13 16:15:55 dpwe Exp $

s = length(x);

if rem(w, 2) == 0   % force window to be odd-len
  w = w + 1;
end

halflen = (w-1)/2;
halff = f/2;   % midpoint of win
acthalflen = min(halff, halflen);

halfwin = 0.5 * ( 1 + cos( pi * (0:halflen)/halflen));
win = zeros(1, f);
win((halff+1):(halff+acthalflen)) = halfwin(1:acthalflen);
win((halff+1):-1:(halff-acthalflen+2)) = halfwin(1:acthalflen);

c = 1;

% pre-allocate output array
d = zeros((1+f/2),1+fix((s-f)/h));

for b = 0:h:(s-f)
  u = win.*x((b+1):(b+f));
  t = fft(u);
  d(:,c) = t(1:(1+f/2))';
  c = c+1;
end;

and the inversion:

function x = istft(d, ftsize, w, h)
% X = istft(D, F, W, H)                   Inverse short-time Fourier transform.
%   Performs overlap-add resynthesis from the short-time Fourier transform
%   data in D.  Each column of D is taken as the result of an F-point
%   fft; each successive frame was offset by H points. Data is
%   hamm-windowed at W pts..
% dpwe 1994may24.  Uses built-in 'ifft' etc.
% $Header: /homes/dpwe/public_html/resources/matlab/RCS/istft.m,v 1.1 2002/02/13 16:16:12 dpwe Exp $

s = size(d);
%if s(1) != (ftsize/2)+1
%  error('number of rows should be fftsize/2+1')
%end

cols = s(2);
xlen = ftsize + cols * (h);
x = zeros(1,xlen);

if rem(w, 2) == 0   % force window to be odd-len
  w = w + 1;
end

win = zeros(1, ftsize);

halff = ftsize/2;   % midpoint of win
halflen = (w-1)/2;
acthalflen = min(halff, halflen);

halfwin = 0.5 * ( 1 + cos( pi * (0:halflen)/halflen));
win((halff+1):(halff+acthalflen)) = halfwin(1:acthalflen);
win((halff+1):-1:(halff-acthalflen+2)) = halfwin(1:acthalflen);

for b = 0:h:(h*(cols-1))
  ft = d(:,1+b/h)';
  ft = [ft, conj(ft([((ftsize/2)):-1:2]))];
  px = real(ifft(ft));
  x((b+1):(b+ftsize)) = x((b+1):(b+ftsize))+px.*win;
end;

EDIT: if you want to start from STFT magnitude, then you run into the more complex problem of phase retrieval. A few paper to start from:

Signal reconstruction from STFT magnitude: a state of the art, DAFx 2011
Non-convex phase retrieval from STFT measurements
STFT phase retrieval: Uniqueness guarantees and recovery algorithms


How to apply baseline methods to MNE preprocessed raw EEG signals? I am new at EEG preprocessing and classification. I have followed Clemens Brunner's blog https://cbrnr.github.io/2018/01/29/removing-eog-ica/ as a tutorial for MNE.

This is how I have preprocessed and singled out all the individual trials for all subjects in the dataset BCI Competition IV Dataset 2a. The dataset can be found here: http://bnci-horizon-2020.eu/database/data-sets

ch_names = ['Fz', 'FC3', 'FC1', 'FCz', 'FC2', 'FC4', 'C5', 'C3', 'C1', 'Cz',
        'C2', 'C4', 'C6', 'CP3', 'CP1', 'CPz', 'CP2', 'CP4', 'P1', 'Pz',
        'P2', 'POz']

mat = loadmat('A01T.mat')

for j in range(6):

    timelist = mat['data'][0][j+3][0][0][1]    //this gives starting of each event
    timelist = list(timelist)
    timelist.append(np.array([96735], dtype='int32'))    //length of eeg for single channel is added
    timelist = np.array(timelist)

    actionlist = mat['data'][0][j+3][0][0][2]    //these are the labels

    eeg = ((mat["data"][0][j+3][0][0][0]*10e-6).T)[:22]
    raw = mne.io.RawArray(eeg, info)
    raw.set_montage("standard_1020")
    raw.filter(1, 30)

    raw_temp = raw.copy()
    ica = mne.preprocessing.ICA(method='infomax', fit_params=dict(extended=True))
    ica = ica.fit(raw_temp, picks=['eeg'])
    ica.detect_artifacts(raw_temp)
    raw = ica.apply(raw)
    trial = raw.get_data()

    label = []
    train = []

    for k in range(48):
        trial_temp = trial[:, timelist[k][0]:timelist[k+1][0]]
        train.append(trial_temp)
        label.append(actionlist[k][0] - 1)

The file is saved in the shape of 22 channels and corresponding measures for the next 7-8 seconds. Frequency is 250Hz.

How do I apply baseline methods such as SVM to this preprocessed data?

Also, how do I calculate the power of the signals using FFT or any other method?

MNE has a fantastic website with a lot of documentation and many tutorial examples. Those examples include how to access machine learning classifiers in scikit-learn (to play with SVMs) as well as the methods for getting / plotting the power, extracting frequency components with FFT. You should definitely check out the website.

Alex


DFT algorithm for FPGA without phase For my FPGA programming lab final project, I want to implement an audio spectrum analyzer. Since I won't be recreating any signals, I only need the frequency magnitude, and not the phase. I'm hoping that ignoring the phase will spare me the trouble of using complex arithmetic, and make the algorithm leaner. I found some references to the discrete hartley transform, which has "no intrinsic involvement of complex numbers" according to wikipedia. I also found the "real DFT" which splits the complex numbers into sine and consine parts. I was wondering, though, how the radix-2 fft butterfly method would be modified for the real DFT or the hartley dft. Or, maybe there's another method to get just the magnitude information that's better for hardware.

Thank you for any advice.

magnitude-only FFT can really only be done with the complex FFT followed by computing the magnitude.

the radix-2 Cooley-Tukey FFT (with the butterflies) will likely be the easiest to implement with a FPGA because of the self-similarity of the butterflies in each FFT pass.  (if you want some simple C code that you might translate into FPGA code, i can supply you with a very simple radix-2 FFT in C.  are you limited to fixed-point or can you do floating point in your implementation?)

if you're displaying spectrum magnitude in dB, then magnitude-squared (without the square root required for magnitude) is easier and just changes the scaling factor on the $\log(\cdot)$ function.

$$\begin{align} 20\log_{10}\left( \Big|X[k] \Big|\right) &= 10\log_{10}\left( \Big|X[k] \Big|^2\right) \\ \\ &= 10\log_{10}\left( \Re\{X[k]\}^2 + \Im\{X[k]\}^2 \right) \end{align} $$

so, for dB, don't bother with the square-root.

if your spectrum analyzer also wants log frequency along the frequency axis, that's another different problem to solve.


Kullback-Leibler Distance of Spectral Data I am currently reading through Music Structure and Analysis from Acoustic Signals and am having some difficulty in understanding how the modified Kullback-Leibler distance is calculated.  (I am just very recently starting to get into audio analysis and back into shape regarding more applied math, so this could just be a case of me getting confused by some symbols!)

I'm specifically having some trouble understanding the Segmentation by Clustering section.  The formula for the proposed KL distance is:

$$KL2(A, B) = KL(A; B) + KL(B;A)$$

Expanded out, this yields:

$$KL2(A,B) = \left(\frac{\mathrm{Cov}[A]}{\mathrm{Cov}[B]}\right) + \left(\frac{\mathrm{Cov}[B]}{\mathrm{Cov}[A]}\right) + (\mu_A - \mu_B) \cdot \left(\frac{1}{\mathrm{Cov}[A]} + \frac{1}{\mathrm{Cov}[B]}\right)$$

Here's where I'm running into issues:

This formula is expected to be used between to matrices of spectral data.  So, basically- how do I compute the covariance?  It seems I do it between each feature vector.  That's easy enough, but how do I combine them so that they're usable in this formula?  It shouldn't be a vector of covariances, right- how could they be divided, then?
Just so I know I'm correct, then the dot product means that the means has to be a vector, right?  If that's the case, then what happens if the segments of audio that I'm looking at happen to be different lengths?  Will the resulting vector from that subtraction just be the length of the longer vector, with the tail end just being either positive or negative entries based on if A or B had more entries?

Apologies if this is confusing!  Any help is appreciated.

The expression in your question seems to have been written for univariate Gaussians. For multivariate Gaussians,
$KL(A,B) = \frac{1}{2}\left[\log\frac{|\Sigma_B|}{|\Sigma_A|} - d + Tr(\Sigma_B^{-1}\Sigma_A) + (\mu_B - \mu_A)^T \Sigma_B^{-1}(\mu_B - \mu_A)\right]$
where $d$ is the dimensionality of the feature vectors.

(https://stats.stackexchange.com/questions/60680/kl-divergence-between-two-multivariate-gaussians)

So

\begin{eqnarray*}
KL2(A,B) &:=& KL(A,B)+KL(B,A)\\
&=& \frac{1}{2}\left[ - 2d + Tr(\Sigma_B^{-1}\Sigma_A) + Tr(\Sigma_A^{-1}\Sigma_B) + (\mu_B - \mu_A)^T (\Sigma_A^{-1}+\Sigma_B^{-1})(\mu_B - \mu_A)\right].
\end{eqnarray*}

I don't think you compute this for the spectrograms themselves but for feature vectors (some of the features may be derived from the spectrogram).


How to append two bandlimited signals and make the result bandlimited without modifying the first signal I need to append (not add) a bandlimited signal F(b) to a bandlimited signal F(a) and keep the result bandlimited without modifying the part that corresponds to F(a).

Both are bandlimited to the same frequency (for example to 1/4 the Nyquist frequency). If F(a) is 100 samples long and F(b) is 100 samples long, the result F(c) is 200 samples in length and the first 100 samples of F(c) are identical to F(a).

An obvious "solution" is to filter the boundary so that the discontinuity can be made bandlimited, but in this case, I am restricted from modifying the samples in F(a).

How can I filter (or manipulate) F(b) such that the result of appending the two signals is also bandlimited to the same frequency and preserves the integrity of F(b) as much as possible?

Edit: For context, here's the problem I'm trying to solve. I am performing procedural terrain generation using 2D heightmaps and need to add a tile of terrain to a pre-existing landscape that may or may not have been procedurally generated. I'm looking for an algorithm that will "smooth" the borders of the tiles so that the tile boundaries are not noticeable.

That's a tricky problem. First you will need to define how to design "integrity". What do you value more: time domain shape, spectrum, phase, envelope, short term Fourier spectrum envelope, etc. You will probably have to trade between those, so prioritizing would be helpful.

In general this is done bu doing "cross fades", i.e. you overlap the two signals for a small time period and then you fade out F(a) while at the same time fading in F(b). The time and shape of the fade can be controlled to satisfy the limited bandwidth requirement.

However, it seems you can't do this, since you can't modify F(a) at all. One potential work around would be to either extrapolate F(a) or use a  circular repetition of F(a) to create enough data for a cross fade. This depends on the exact properties of F(a)


PSD of complex white gaussian noise It may be a really simple question, but I'm not sure about this one:

Given a complex white Gaussian noise process with iid real and imaginary parts and a double sided power spectral density of $N_0/2$. What is the PSD of this random process?

Is it just:
$$S_n(f) = \frac{N_0}2 + j \frac{N_0} 2\qquad ?$$
Or am I wrong?

With reference to $N_o$ this usually is the symbol for the power spectral density of thermal noise, where $N_o = kT$, where k is Boltzmann's Constant and T is the temperature in Kelvin. The thermal noise signal is a complex, white Gaussian distributed noise, with half of the power in the real component and half the power in the imaginary component (but power in itself is not a complex quantity). Therefore if the total power density is $N_o$, the power density of the real component would be $No/2$ and the power density of the imaginary component would also be $N_o/2$.

Considering the complex signal we can represent the spectral density as "Single-Sided" or "Double-Sided". With a Single-Sided power spectral density, the frequency axis starts at 0 and extends to infinity, making use that the the negative frequency components have the same density as the positive. We do not need to show this redundant information, but we do need to double the quantity that is shown to account for the power that extends into the negative frequencies. If we include the negative frequency axis, as done in a Double-Sided power spectral density plot, then the density is given as $N_o/2$.

In the end in either case the power density is $N_o$ Watts/Hz. In the Single Sided case the power in 1 Hz of bandwidth is completely shown on the positive axis, while in the Double Side case, half the power is in the positive axis and half is in the negative axis.

Thus in this we have made two considerations in how the power may be distributed; either in the positive and negative frequency axis if we present Single or Double-sided spectrums, and in the real and imaginary axis if we consider how the power in the full complex signal ($N_o$) is distributed between real and imaginary components. For example, if we were only to consider the real component, the power spectral density would be $N_o/2$ and then if we viewed that as a Double Sided PSD, the PSD would be $N_o/4$.

Single and Double Sided power spectral densities are frequency encountered when working with phase noise, which may provide a further intuitive way of seeing this: The double sided PSD for phase noise is given as $\scr{L}_\phi(f)$ and represents what we would read directly off of a spectrum analyzer, where we can observe both the upper and lower sidebands of the spectrum (although we typically only report one side on the plot which causes confusion in if it should be called single or double sided; one side is shown but it is a double sided power spectral density quantity--- we just don't need to show the other side since it is identical). The single sided power spectral density $S_\phi(f)$ represents the power from both sidebands as given on a plot with a positive frequency axis only (since the negative sideband is equivalent to the positive, we can put it all on the positive frequency axis simply by doubling in power what we read off of the spectrum analyzer).  So $S_\phi(f)$ represents all the power and is used in related calculations where we integrate over a band of frequency offsets, while $\scr{L}_\phi(f)$ is convenient in that we can read it directly off of test equipment (spectrum analyzers).


Question on FIR filter I've been trying to do this problem for the last few hours. Could someone show me the workings on how it is done? Not just the answer.

The figure shows the structure of a FIR filter where $A = 1.5$, $B = 0.69$ and $C = 3.0$.

If $x[n] = u[n]$, what is the value of $y[n]$ for $n = 20$?

[EDITED] From left to right, try to find the impulse response. On the bottom branch, delay the unit pulse $d$ and get $d[n-1]$. Add $Ad[n]$. Delay twice and get $d[n-3] +Ad[n-2]$. Add $Bd[n]$. Add $Cd[n-1]$. So the FIR coefficients of the impulse response $h$ will be $[1,A,C,B]$, but in the reverse order (as correctly pointed out by colleagues), because $y[n] = \sum d[k]h[n-k]$. So your system is causal, and $h[0]=B$, $h[1]=C$, $h[2]=A$, $h[3]=1$.

However, at $n=20$, you are in a stabilized region of the signal (all the known samples are equal to $1$, so even if you make a mistake as I did first, the answer will be the same, $1+A+C+B$, as stated before by other colleagues.


FFT results are repeating? I am trying to detect notes out of a musical sound input. My program can be run here. The problem is the patterns in the frequency response is repeating:

I suspect that echoing in the room is causing the higher frequency components. If I try a sound generator such as this, repetition is much less, but it is there. Try for example 165 Hz (note E3).
My primary question is that, how come sometimes repetitions have higher spikes in FFT than the actual signal (see the above image). This happens once in a while when human voice is tested. I know the actual signal is the 1st one because I test it in another software (Singer's Studio on iPhone).

Interesting code! This was tough to follow from the pictures, but I did run the OP's code from the linked site, I can't play the sound loud enough to really test but appears to not be very sensitive for lower volume signals. There is no indication of the frequency axis on the graphic but I assume the left side of the screen is DC and the right side of the screen is half the sampling rate.
The sampling rate in the code was 4200 (I assume Hz) which would be way too low for audio and suggest the OP is seeing the aliasing of all higher frequencies above the Nyquist frequency of 2100 Hz. Since the audio system has a likely bandwidth of 15 KHz, any signals in multiples of 2100 Hz up through that frequency range would fold in and be a more obvious explanation of multiple frequencies due to aliasing.
I did change sampleRate to 44100 and it appeared (in my minimum ability to test) to behave more as expected.
Room echoes on their own should not be non-linear meaning they would not create new frequencies but only modify the phase and amplitude of the frequencies that already exist.


How do I filter audio with modified spectrogram? Say I computed Spectrogram of some audio input using scipy.signal.spectrogram.

Is there an easy way in python to use it as a filter for my audio input to produce an other audio which would basically be just silence or very near to that?

In practice I would modify that spectrogram before applying it back to the audio to filter out only some audio information.

Modify the spectrogram and then use scipy.signal.istft

scipy.signal.istft(Zxx, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, input_onesided=True, boundary=True, time_axis=-1, freq_axis=-2)
Perform the inverse Short Time Fourier transform (iSTFT).