Recover IQ samples from signal I have a data acquisition device that provides me with an array of integers that represents a modulated signal (let's say QAM16).

I would like to recover the I and Q samples from the raw time-domain signal, so that I can demodulate the data. Is this possible? If so, how would I do it.

If you just want to demodulate, then it is not necessary to recover I and Q components. What you have is most likely the QAM symbols. For example, for $16$-QAM they are $16$ different integer values that are assigned to different constellation points. You need to know the type of assignment (i.e. coding) that is used in the QAM constellation at the modulator.

At the $16$-QAM modulator, every $\log_2 16=4$ data bits are mapped to a certain constellation point (which ultimately is identified by an integer). For instance, $0101$ can represent the symbol $5$, and so on. Gray coding is usually used for this purpose, but it is not always the case. So you need to find out the exact constellation mapping and then, the demodulation process is nothing more than the inverse of that mapping. So when there is a $5$ your data bits are $0101$.


How to label sound events in noisy .wav files I'm trying to find a way of automatically identifying birdsong in phone-recorded sound files (1min long). Currently the algorithm I'm using doesn't label all of the events I want it to. (it is designed for non noisy audio files).

I don't have much experience in DSP. I'm hoping that upon looking at my wave-forms someone may have an idea of how to pick out these 'birdsong events'!
I've shown one of the waveforms. The rest can be seen here: imgur.com

Any method would be appreciated!

I've thought about applying a noise reduction algorithm but I suspect this could just interfere with things.
Here is a link to the algorithm I'm using in case you're interested: Github/kylemcdonald

The labelled events are the cross hatched ones.

Thanks in advance for taking the time to look at this!

Automatic detection of birdsong is a well studied problem in machine-learning applied to ecoacoustics/bioacoustics. State of the art methods use Convolutional Neural Networks on spectrograms. In binary classification of bird/no-bird they reach AUC ROC of around 0.90 in noisy and mismatched conditions. The DCASE2018 Bird Audio Detection challenge is a good starting point. Open-source code is available for many of the systems, including the challenge winner. Several of the systems could be adapted to Audio Event Recognition to give time-information of predictions.

Note that to verify your approach is working well, you should always also have a manually labeled test-set.


Is the fundamental basis function of the DFT not an intregral number of cycles? When the DFT is defined and described I typically see the basis functions described as an integral number of sine and cosine waves.

But the basis functions are $\sin(2\pi k (n/N))$ or $\cos(2\pi k(n/N))$ where $k$ is the harmonic, $N$ is the number of samples, $n$ is the time domain sample number.

So the right hand boundary of this argument is  $2\pi k(N-1)/N$, which seems to be always a bit short of a full cycle.

Could someone please confirm or explain my error?

The picture below shows a sine basis function, apparently ending short of a full cycle.

First, a small correction:
The basis functions of all Fourier transforms are of the the type $e^{2j\pi k\frac nN}$, which can also be interpreted as $\cos(2\pi k \frac nN) + i\cdot\sin(2\pi k \frac nN)$, but never only real-values sines or cosines; that'd lead to large problems on the algebraic structures around. But that really doesn't matter here:

You're absolutely right, the "completing" sample of the oscillation is always missing – but that's absolutely logical, if you think about it – when you periodically continue the signal you're transforming, you'd get twice the same value (the "start" value) after each other, and that mustn't happen.

Your picture is an excellent illustration of this:

The oscillation is observed for exactly two periods, and each period is 16 samples long; that means that for the $n$th sample, the $(n+16)$th sample must have the same value. For example, the $n=16$th has the same value as the $(16+16)=32$nd sample – which, by principle of 32-periodicity is the same as the $0$th sample.


Optimization of harmonics calculation I need to compute the sine and cosine of an argument along with n "harmonics"
\begin{matrix} \sin(x) & \cos(x) \\
  \sin(2x) & \cos(2x) \\
\cdots \\
\sin(nx) & \cos(nx)
\end{matrix}
This takes a lot of computing time and since this is a real-time system, I need to optimize this calculation.
Here's my solution so far :
I defined a rotation matrix $$
R = \begin{pmatrix}
    \cos(x) \ -\sin(x) \\
    \sin(x) \ \cos(x) \\
    \end{pmatrix}
$$
I compute $[\cos(nx); \sin(nx)] = R  [\cos((n-1)x); \sin((n-1)x)]$; I simply iterate to compute the missing $n-1$ harmonics.
This method is 3 times fast than computing all harmonics separately.
I'm curious, are there better solutions? I know a look-up table could be even faster but I would prefer not to go this way.


Detect repetitive units within signals I have several signals that consist of repetitive units. In the figure you'll clearly see the variability of the signals, that increases top down. The first signal is super repetitive and units are indicated with green lines.
In the third, you'll see peaks and in the middle a little insertion that I know consists of rather diverging units, which still are units however (red lines).
The remaining three signals display the variation more.

Which signal processing / machine learning tool should I use in order to detect these units? With thresholding it works to find the significant peaks, but once signals get funky it's really difficult to accurately detect unit positions.

Edit:

I made some progress with spectral analysis and plotted a filtered signal over the initial signals. Indicated by the blue arrow is a region variable repeats. Here, the amplitude of the original and filtered signal do not match. Same for the red arrow. Similar effects for signal further down.

Your repetitive signal seem to be in a lower rate than the noise. I would run an FFT (why short) on all 6 sequences and filter out (zero) the frequencies that are clearly attributed to noise and then run an inverse FFT.
If you know the repetition rate you might select the one frequency that has the correct information and use the phase information to locate the desired pick. I will try to include exactly 20 cycles if possible (or 16 or 32).


Audio compression with STFT Let STFT[frame, freqbin] be the short-time Fourier transform of some audio data.
(here frame is the index in time, and  freqbin is the index in frequency).

Very often, STFT[k, f0] will be close to STFT[k+1, f0]...

... then it could be interesting to store, for each frequency bin f0, only the differences :

STFT[0, f0], STFT[1, f0] - STFT[0, f0], STFT[2, f0] - STFT[1, f0], STFT[3, f0] - STFT[2, f0], etc.

(such processing is known as delta encoding) ; these differences will be small on average, and this could lead to a reduction of the average bit usage for storage of the audio data.

Can this idea be turned into an efficient audio data compression ?

(I don't speak here about psychoacoustic effects that could lead to use different levels of quantizations for each frequency / or frequency masking ; I know that such ideas are used in MP3 compression).

Edit: as mentionned by @pichenettes, the graph above showing continuity between adjacent frames is the "magnitude plot". The "phase / angle plot" is a bit less continuous :

I see two big drawbacks to your method:

The image you posted is not the STFT, but its magnitude. The difference between the magnitude of adjacent STFT frames is small; so for your scheme to work you would have to store the magnitude and phase separately. But then, there would be no compression at all on the phase data that would still represent 50% of the data and would be required for reconstructing the signal. Try plotting the phase (angle) of the STFT to see how ugly it is.
Assuming the deltas are encoded using variable-length integers or Golomb codes (I assume that's the idea behind it...), your method would compress large and temporally stable values in the spectrum (the red "lines" on the spectrogram image). It's obvious that these represent only a small fraction of the values. There would be no gain for the background "noise" which amounts for a very large portion of the spectrogram magnitude data.

A much more efficient way of getting rid of these stable/predictable peaks in the spectrogram is to whiten the signal and code the residual+the filter coefficients. This is (among other things) the principle of the FLAC lossless codec.


Parameter extraction from DFT results, possibly with leakage effects Let's start with a simple temporal law, that I want to analyze in the frequency domain :

$x(t)=\frac{1}{2}(A+Bcos(2\pi \nu_0 t+\phi_0))$, with :

$A=4$
$B=10$
$\nu_0=10^8$
$\phi_0=10^8$

I made a simple DFT (in Ruby) with $N=256$ samples at $f_s=6.10^8$ and got the following results, that are self explanatory, hopefully :

{:freq=> 0, :ampl=>1.9661708933907966, :phase=>0.0}
...non significant amplitudes...
{:freq=>41, :ampl=>0.8468197908638829, :phase=>0.502386551741981},
{:freq=>42, :ampl=>2.0871309982776522, :phase=>0.5146361931223131},
{:freq=>43, :ampl=>4.115492878192129,  :phase=>0.5268855235929356},
{:freq=>44, :ampl=>1.0144367909506748, :phase=>0.5391345296503323},
...non significant amplitudes...

I clearly retreive A : $1.96 = \frac{A}{2}$ gives $A = 3.9 \approx 4$. Good.

As for $\nu_0$, this seems correct also :
$f_{analysis}(43)=\frac{43f_s}{256}=100781250.0\approx 10^8=\nu_0$.

For $B$ and $\phi_0$, things are more difficult, because (I guess) of DFT frequency leakage. Direct reading of my result gives $$2\times 4.11=8.22 \not\approx B=10$$

As experts, how would you get acceptable $B$ and $\phi_0$, from my results ?

Interpolation is one method used to estimate parameters of non-periodic-in-FFT-aperture signals.  Parabolic or, better yet, Sinc interpolation will help with estimating both the frequency and the amplitude peak B.  For phase interpolation, I find it easier to use something like a pre or post fftshift to move the phase 0 reference to the middle of your data, then interpolating an estimated phase from the fft, then offsetting the phase back relative to wherever using the frequency estimate.


Compressive Sensing: Reconstruct Gap in Antenna Array I have an antenna array with $N$ elements spaced half a wavelength apart.
I have a second, identical antenna array that is the distance $D$ apart from the first one.

Could I use compressive sensing to reconstruct the full antenna array?

If so, how big can $D$ be and what are my limitations?

In general: no. What you can and cannot reconstruct will depend a lot on your model.

If you assume that you have:

point sources in the far field of the antenna array
these point sources emit narrow-band, mutually uncorrelated signals
you can observe them over a long window in which they are stationary
then you can use co-array processing to virtually fill the missing elements between your two arrays.

Let's see how it works for a single source and let us emit noise for now for clarity (it's not hard to add it later). The observation vectors in your two subarrays can be expressed as $\mathbf x_i(t) = \mathbf a_i \cdot s(t) $, $i=1,2$, where $\mathbf a_i$ are the array steering vectors and $s(t)$ is the emitted signal. Let $\mathbf x(t) = [\mathbf x_1(t)^T, \mathbf x_2(t)^T]^T$ be the vector of all observations. Then, its covariance matrix for a single source is given by $$\mathbf R = \mathbb{E}\{\mathbf x(t) \mathbf x(t)^H\} = \mathbf a \cdot \mathbf a^H \cdot P_s, $$ where $P_s$ is the variance of $s(t)$ (remember it was assumed stationary) and $\mathbf a = [\mathbf a_1^T, \mathbf a_2^T]^T$. Now, the trick is to vectorize $\mathbf R$. Using the Kronecker-vectorization identity, this gives $$\mathbf r = {\rm vec}\{\mathbf R \}= (\mathbf a^* \otimes \mathbf a) \cdot P_s,$$ where $^*$ denotes complex conjugation and $\otimes$ is the Kronecker product. Now, this looks like a single observation with an array that has an array steering vector given by $\mathbf a^* \otimes \mathbf a$, which we refer to as the co-array of $\mathbf a$.

For a ULA with half-wavelength spacing, your array steering vectors take the form $\mathbf a_1 = [{\rm e}^{0}, {\rm e}^{\jmath \mu}, {\rm e}^{2\jmath \mu}, \ldots, {\rm e}^{(N-1) \jmath \mu}]^T$ and $\mathbf a_2 = [{\rm e}^{D \jmath \mu}, {\rm e}^{(D+1)\jmath \mu},  \ldots, {\rm e}^{(D+N-1) \jmath \mu}]^T$, where $\mu = \pi \sin(\theta)$ is the electric angle. Therefore, $\mathbf a^* \otimes \mathbf a$ contains terms ${\rm e}^{p \jmath \mu}$, where $p=-(N-1), \ldots, (N-1)$ (from $\mathbf a_1^* \otimes \mathbf a_1$ and $\mathbf a_2^* \otimes \mathbf a_2$) and $p=-D-N+1, ..., D-N+1$ and $p=N-D+1, ..., N+D-1$ (from $\mathbf a_1^* \otimes \mathbf a_2$) [give or take a $\pm$, I guess you get the idea].

Given that, for $D<N$, you've created a virtual ULA that is completely filled. I'll let you do the math how many elements it has in total.

Note that you will in general need to apply spatial smoothing if you want to estimate more than one source.

That said, compressive sensing is the wrong term here. It is co-array processing, which is somehow linked to sparsity-based processing (which became popular in the compressive sensing field, leading to many people confusing the two).


Do FIR filter have any application in control theory? I know about IIR filter as discrete pendent to transfer functions in the Laplace domain. So it is actually quite simple to convert the function of the control system and finally receive the discrete output vector. (Bilinear Transform...)

But, how about FIR filter? Can we use them  to approximate functions of control systems.

FIR filters are fairly common in some areas of control theory. As they usually incur a lot of added phase/time-delay, they are not really usable in the feedback path of regular control systems, but they are useful when the added phase/time-delay is not affecting the system in an adverse way, or when the particular phase response and time-delay is desired.

Examples:

Feed-forward control. FIR filters are useful for producing filters that approximate arbitrary frequency responses, hence they can be used to shape a reference signal. A typical example is to use an FIR filter with the inverse frequency response of the plant -- trying to counteract the dynamics of the plant in order to get a desired output. Phase/time-delay is not interfering with the stability or performance since the computation can be done offline. FIR filters can often produce higher performance than IIR filters, especially where there are non-minimum phase zeros.
Learning systems (iterative learning control [ILC] and repetitive control [RC]). In these systems, again, FIR filters are often used to generate an inverse response filter of the plant. This is done to improve the performance of the learning, as this enables a higher gain to be used for a larger bandwidth. The time-delay is here desired (RC), or it does not matter, as the computations are done offline (ILC). Low-pass FIR filters are also used here for their linear-phase properties, as the linear-phase property makes the filter not interfere with the phase of other filters, and this improves performance.


How do I find the phase response of a filter from the filter equation? I am new to this forum. I have been teaching myself DSP for the last two years, so my experience and knowledge has a long way to go.

I am trying to find the phase response of a specific filter that I am using on Csound for any given fundamental frequency of the input (and I guess I'll need to be able to do that for any cuttoff frequency and ideally any sampling rate as well).

The filter is a 1st order lowpass IIR, and the filter's equation is as follows:

$$y[n] = c_1  x[n] + c_2 y[n-1]$$

$h_p$ - half power point
$b = 2 - \cos(2 \pi h_p/s_r)$
$c_2 = b - \sqrt{b^2 - 1}$
$c_1 = 1 - c_2$

Three questions:

Can I find the phase response for any frequency based on this equation alone?
Can I derive a phase response equation from this equation alone?
Is there a standardised way of finding the phase response for any lowpass filter, or does it depend on what kind of filter it is?

I have found some info about finding phase response both in these forums and elsewhere online, but without knowing the answer to question 3, it's hard to know if they are relevant or not.

I have also read most of the scientists and engineers guide to DSP (which I notice is also the boards DSP guide, awesome!), and I have been going back over it the last few days to try and get some clues. The best solution I can find is to fft the impulse response, and I know how to find the impulse response for a cuttoff frequency, but it'll take a long time to calculate the impulse response for all possible cuttoff frequencies.

EDIT:

I have been reexamining ths online resouce from analog.com:

http://www.analog.com/library/analogDialogue/archives/43-09/active_filters.html

The website tells me this about the phase response of a 1 pole lowpass filter:

Is this true of any 1st order low pass filter, or are they refering to a specific filter type? If this is true of my filter, I think I might be able to do this (though I won't say for sure yet!)

You need to determine the frequency response of the filter. The phase of this complex function is the filter's phase response. If you have the filter's difference equation, you need to transform it to the $\mathcal{Z}$-transform domain, where a delay of $k$ (with integer $k$) corresponds to a multiplication by $z^{-k}$:

$$Y(z)=c_1X(z)+c_2Y(z)z^{-1}\tag{1}$$

From $(1)$ you can determine the transfer function

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

The frequency response is the transfer function evaluated at $z=e^{j\omega}$, where $\omega$ is the normalized frequency in radians:

$$\omega=\frac{2\pi f}{f_s}\tag{3}$$

with $f_s$ the sampling frequency, and $f$ the actual frequency in Hz.

From $(2)$, the filter's frequency response is

$$H(e^{j\omega})=\frac{c_1}{1-c_2e^{-j\omega}}=|H(e^{j\omega})|e^{j\phi(\omega)}\tag{3}$$

where the phase response $\phi(\omega)$ is given by

$$\phi(\omega)=\text{atan2}\left(\frac{-c_2\sin(\omega)}{1-c_2\cos(\omega)}\right)\tag{4}$$

For more complicated filters it's usually easier to determine the frequency response (and the phase response) numerically by using an FFT. The (sampled) frequency response is obtained by element-wise division of the FFT of the numerator coefficients of the transfer function by the FFT of the denominator coefficients. The (sampled) phase is then obtained as the argument of this complex-valued vector. In Matlab/Octave, the frequency response can be obtained directly from the filter coefficients using the command freqz.


Please help me understand this paper about Discrete Wavelet Transform! I am trying to understand this paper says.

I have found this paragraph to be vague:

  (a) Preprocessing Stage. The first step in this stage was to apply the
  Discrete Wavelet Transform (DWT). This tool is based on the
  decomposition of a signal in subbands by means of the use of a pair of
  digital filters (low pass and high pass filters). The outputs of the
  low pass filter are named approximation coefficients (), while the
  outputs of the high pass filter are named detail coefficients (),
  where  represents the level of the subband. This process of
  decomposition through filtering is repeated  times. In each iteration
  the signal is subsampled by a factor of 2. In practice, the DWT is
  implemented with the Mallat pyramid algorithm [28]. Some studies have
  shown that the use of a four order Daubechies wavelet is one of the
  most effective when processing ECG signals [29]. Afterwards, the
  energy percentage of each level was calculated (see (1)-(3) and Figure
  3), then the four ones with more energy were selected to reconstruct
  the signal (Figure 4). In this way it was assured that the levels with
  more information of the ECG signal were selected because noise or some
  interferences such as those of the electrical network or the artifacts
  are usually found at low energy levels (high frequencies, generally
  between  and ) (see Figure 3). Therefore, using energy levels to
  discriminate the noise of QRS complex was a good option. Finally we
  proceeded to remove the offset of the signal by leaving out the
  approximation coefficients (). In the present study a db4 mother
  wavelet with 7 levels of decomposition was selected. The number of
  levels was selected because the data were processed in buffers of 1024
  and thus the number of iterations allowed was 7.

Specially this part:

  Afterwards, the energy percentage of each level was calculated (see
  (1)-(3) and Figure 3), then the four ones with more energy were
  selected to reconstruct the signal (Figure 4).

the energy percentage of each level was calculated - the energy of what, of the low frequency coefficients or the high frequency coefficients or both? I guess the high coefficients, or detail, right?

then the four ones with more energy were selected to reconstruct the signal

how?

On the text the researcher says the signal was decomposed in seven levels, so, he have obtained seven low frequency coefficients and seven high frequency coefficients.

He says he have used detail level 4, 5, 6, and 7 but he mentions these numbers in the following order: 5, 4, 6, and 7. Why?

What is he talking about? Reconstructing the signal just using 4 levels of IDWT? If he does that the reconstructed signal will have 3 levels missing, or in other words, it will be 1/8 of the original frequency. I the signal is being sampled at 30 Hz, the final reconstruction will have 3.75Hz.

Is this what the paper is saying? and why have the researcher mention the levels out of order?

Basically, the method decomposes a signal into one approximation $A_7$ and seven details $D_7$,..., $D_1$. Since the wavelet is orthogonal, the energy of the signal is equal to the sum of the energy of each subband. To me, equation (1) is not correct, one should not sum approximations over levels. Approximation coefficients $A_7$ are discarded, they kind of remove the low frequency or baseline parts. They keep the four most energetic detail subbands, which they order by decreasing energy. Then they set the others to zero, and reconstruct the signal by inverse wavelet transform. So no level is missing, they are just erased. This turns out to be a very crude form of denoising, as they finally only keep the central frequencies, in the interval $[1/2^7\; 1/2^3]$. This could be a misuse of wavelets, as a simple bandpass filter could have been better.


Calculate maximum of filter kernel I'm sure there must be an easy way to do this, but given the Fourier transform of an isotropic filter kernel, $\hat{f}(\mathbf{u}) = \mathcal{F}f(\mathbf{z})$, can one calculate the value of the kernel at $\mathbf{z} = 0$?

Since $$f(\mathbf{z})=\int_{\mathbf{R}^n}\hat{f}(\mathbf{u})e^{2\pi i\mathbf{z}\cdot\mathbf{u}}\;d\mathbf{u}$$

$$f(\mathbf{0})=\int_{\mathbf{R}^n}\hat{f}(\mathbf{u})\;d\mathbf{u}$$

So you simply integrate (or sum in the discrete case) over $\hat{f}(\mathbf{u})$.


Phase shift between two signals So I got some oscilloscope captures for a project I'm doing and I'd like to find a phase shift between them because I don't trust the scope calculation.

So I extracted data from a .csv file and loaded it into Matlab and now I'd like to find a phase shift using FFT (I think). Is there a way to do this even though the second signal, output, is not a perfect sinusoid?

There are many different ways to extract phase shift, the simplest one I think is using normalized cross-correlation using instruction 'xcorr' and then finding the index where the maximum correlation is placed. Other technique as you mentioned is using FFT and spectrum multiplication and ifft(equals to correlation in time) and/or calculation of  arctan(complex/real) of bins. Other methods like Hilbert transform is comprehensively detailed in below links:

https://stackoverflow.com/questions/27545171/identifying-phase-shift-between-signals

Best method to extract phase shift between 2 sinosoids, from data provided


How to compare binary signals with different lengths I have bunch of 1d arrays of 0 and 1 with specific length as data. also i have a key signal(another 0,1 array) with shorter length than my data arrays.
Now how can i have a good measure of similarity or dissimilarity between those long arrays with the shorter key array?

you need to elaborate on what you mean by similarity because there are many possibilities.

If you key sequence is embedded in any of your test sequences, the cross correlation is as was suggested by one of the commentators is a good approach, which forms a point by point product between the test and key sequence and sums the results of the product and provides a number, which if the key has $n$ ones the sum of products equal to $n$ indicates an exact match.  The key sequence is shifted relative to the test and the sum of products is calculated for all possible shifts.  This also works for situations where there are error bits in the tested sequence, when the sum product is near $n$.

One might presume that a sum product equal to zero is an indication of dissimilarity, and could be depending on your application but it is usually preferable that measures of similarity should be invariant under information preserving transformations of the data. The replacement of each one by a zero is an information preserving transformation, so a cross correlation value near zero in this case should be considered to be very similar.

A reverse ordering of test sequences is another information preserving transformations.

It is usual to consider dissimilarity as a situation where  a one in the reference would correspond to either a one or zero in the other. In this case dissimilarity to equivalent to statistical independence.

It is sometimes useful to transform $(0,1)$ to $(-1,1)$ when calculating cross correlation.

You can also compare sequences by the ratio of ones to zeros. Run length statistics are also useful.

Cross correlation is compatible with statistical hypothesis testing and confidence bounds. These calculations are covered in any good calculus level statistics text book (like Bhattacharyya and Johnson)


Isn't the rect pulse filter implicit at the transmitter? I've heard about shaping the symbol with a rect filter at the transmitter, but why? Isn't this just implicit from the fact that it's going to be transmitted as a symbol truncated to the symbol period regardless? Or is that simply what they're referring to by shaping it with a rect filter i.e. isolating the symbol to that window. Am I right in thinking it's just the default. Also does the receiver extend the signal as if it were infinite when performing the fourier transform and therefore the rect pulse shaping at the receiver is simply the receiver not doing that?

In most waveforms for data communication we do NOT want to transmit a rect for each pulse given the excessive bandwidth required (since a rect is a Sinc function in frequency which is relatively very wide-band). This is the motivation for pulse shaping (such as the commonly implemented root-raised cosine pulse shape): to limit the transmitted spectrum and improve the overall spectral efficiency of the waveform.


How to efficiently compute a sliding FFT/DFT with a variable window? I aim at performing a sort of sliding FFT/DFT on an audio input with a variable time-window, that depends on frequency.

I am aware of sliding FFT/DFT algorithms. But instead of having a window of constant duration 100 ms or whatever, I would like to make the window shorter for higher frequencies and longer for lower frequencies. So: the idea is that higher frequencies would appear and vanish quickly in the FFT output than lower frequencies. Is such a thing doable as efficiently as traditional sliding FFT/DFT? Are there algorithms for this already?

When using the sliding DFT instead of an FFT, you can actually have a different window length per bin. This is because each bin is computed separately anyway. The only constraint is that each window length must be a multiple of its bin wavelength (or half of its bin wavelength if you modify the transform a little). For example, you can do a constant-Q sliding DFT by using the same factor for each bin, and spacing the bins appropriately. More complex variants are possible but require a lot of tweaking, depending on what you want to do with the result. The different windows should be aligned at their center.

I once found a paper about this (after coming up with essentially the same thing myself); it's called "Sliding with a constant Q" and can be googled. For resynthesis, alternatively to the method described in that paper, the trivial resynthesis algorithm for the sliding DFT (essentially just summing everything) also works for variable window lengths.


Is there any literature discussing PDF after quantization? Consider a random variable with a known distribution. I want to know if we quantized the variable the how the probability density function(PDF) would change? Specifically if the original PDF is Rayleigh distributed what would be the distribution would look like after the quantization?

any discrete random variable has a p.d.f. that is a summation of dirac delta functions.

$$ p_\mathrm{y}(\alpha) = \sum\limits_i P_i \ \delta(\alpha - y_i) $$

where $\sum\limits_i P_i = 1 $.

if $y[n]$ is the quantization of $x[n]$:

$$ y[n] = \Delta \bigg\lfloor \frac{x[n]}{\Delta} \bigg\rfloor $$

where $ \lfloor \cdot \rfloor$ is the floor() function and $\Delta$ is the quantization step size.

then $y_i = i \Delta $ for some set of integers $i$.  the p.d.f. weighting constants become

$$ P_i = \int\limits_{i \Delta}^{(i+1)\Delta} p_\mathrm{x}(\alpha) \, d \alpha $$


Slope of 2nd order low pass vs band pass Butterworth filters I have noticed that some manufacturers of audio processing products label a 2nd order Butterworth low pass as having a slope of 12dB per octave, but then also label a 2nd order Butterworth band pass as having a slope of 12dB per octave. In actuality each side of a 2nd order Butterworth band pass has a slope of only 6dB per octave. Is there a convention that because it is symmetrical that the slope of each together constitute the "slope" or is it just that they are wrongly assuming that order is equivalent to slope for all the Butterworth filter types?

By convention, the pole count refers to the number of poles in the design polynomial (Butterworth, Chebyshev, etc) and this defines the amount of roll off in the stop band.

If we maintain this convention, then all filters roll off at the rate of 6 dB / octave / pole in the stop band, whether it is a low pass, high pass, band pass, or notch.

Also, a nth order Elliptic has the same ultimate roll off rate in the stop band as a nth order Butterworth. Of course, the Elliptic's transition band will be much narrower.


How one apply correctly FFT in image denoising I'm writing program (Qt widgets/c++) for removing noise from images. As denoising method, i selected non local means method. This method has incredible quality of restored images (that's why it's the only denoising method in OpenCV), but has huge computation cost, so i made a lot of modified variants of this method (some with multithreading, some algorithmic). But, i'm having problem with the one, involving FFT

I followed all the steps of this article (only one page, 1430) and all works perfectly, except for FFT part, there just 2 lines about it in the paper and i can't understand, HOW should one use fft

This problem has bothered me for months, any help or insight would be greatly appriciated.

Shortened version of question:
  How can i get summed squared difference of two arrays on the image (the one at top and one in the middle, values are colors) quickly? ( O(n^2) is huge cost, there are lots of this kind of operations, paper above states, that it can be done via FFT with O(n*log n) (says that this 2 arrays forming circular convolution somehow))

The trick inside the paper is the following:

What you want to compute is $\sum_{i \in W} |I(x+i)-I(y+i)|^2$, where $I$ is an image, $x$ and $y$ two noisy pixels and $i$ is a 2D offset used to define a patch.
Expanding the expression yields: $\sum_i I^2(x+i) + \sum_i I^2(y+i) - 2 \sum_i I(x+i)I(y+i) = A + B - 2C$.
$A$ and $B$ are computed using a squared integral image, i.e., an integral image from the squared original image.
$C$ is the convolution between the two patches centered on $x$ and $y$. Thus, it can be computed in the Fourier domain, where it becomes a multiplication. You get the value of $C$ by computing the Fourier transform of the patch around $x$, the patch around $y$, pointwise-multiplying these results and taking the inverse Fourier transform of the multplication result.

The Fourier transform is obviously a 2D transform since you are working with 2D data. What you obtain for a given patch is a 2D array of complex values.

Additional notes

In my opinion this article is not the best NL-means speedup strategy.
Experiments I did way back in 2007/2008 show that pre-selection of patches are better (both in terms of speed and quality of the results).
I have started blogging about these here, but unfortunately I am looking for time to finish the posts.

The original NL-means papers mention blockwise implementations that can be interesting. There are fundamentally 2 ways in implementing NL-means:

writing a denoising loop for every pixel in the image
writing a denoising loop for each patch, then back-project the patches to form an image.

The first impolementation is the original approach, because in 2005 memory and multicore CPUs were expensive. I chose on the other hand number 2 on recent hardware in the past 2 years. It depends on your typical image size and if you want to be able to compute domain transforms like DFT/DCT (as in the proposed paper and in BM3D).


How does Overlap-Add work for IIR filter? So let's say I'm trying to implement something like an LPC vocoder. I analyse a speech signal by breaking it up in small chunks and determining their LPC coefficients, which are by design, the coefficients of an all-pole (IIR) filter.

So far, I understand that for an FIR filter, the impulse response is a finite M samples, so when performing an overlap-add (OLA) block convolution, each signal of length N is extended into a sequence of N+M-1 length, and for a naive implementation, we can overlap and add these M-1 samples in the tail, with the next chunk.

If I want to process a buffer of audio with the IIR filter coming from the LPC analysis, the IR should be infinite, so the only way I can think of using a block convolution is by convolving a delta with the IIR to get a long FIR and truncating it (which I think is terrible).

I also see pages like this that say to use the Block Recursion method, which seems to establish the IIR filter coefficients as columns of a matrix (states?), and that if it is causal the upper-right triangle of this matrix is filled with zeros. I don't understand how this works or why such a matrix is needed?

According to wiki, it seems that overlap add technique usually works for FIR filters. I hasn't read the link you posted and the matrix form for IIR before. But when it comes to buffer for IIR, I have an alternative (hope simpler) method which might works in your case.

  If I want to process a buffer of audio with the IIR filter coming from the LPC analysis, the IR should be infinite, so the only way I can think of using a block convolution is by convolving

The concept is filtering in sections and utilizing the memory in filters. Let's first look at the testing code and usage with following steps,

Dividing the x in to x1 and x2,
Filtering x1 and x2 seprately (with memories in filterszf)
Combining the ouput together.

The results confirm that filtering in section can be the same as x directly.

%% Filter in section
x = randn(10000,1); % psuido signal

x1 = x(1:5000);
x2 = x(5001:end);

b = [2,3];
a = [1,0.2]; % IIR filter coeffs

[y1,zf] = myFilter(b,a,x1); % The zf is the memory in IIR filter!
y2 = myFilter(b,a,x2,zf); % assign the zf when filter the next section.

y = myFilter(b,a,x); % filtering x directly.

isequal(y,[y1;y2]) % Verification

Here is the implementation of myFilter Function. In fact, myFilter imitates(implements) the built-in function fitler in matlab. The implementation architecture is Direct form II, which has advantages over least usage of delay operator and disadvantages over not free for numerical overflow for fixed points.

function [Y, z] = myFilter(b, a, X, z)
    n = max(length(a),length(b));
    a = [a zeros(1,n-length(a))];
    b = [b zeros(1,n-length(b))];
    z(n) = 0;
    b = b / a(1);
    a = a / a(1);
    Y    = zeros(size(X));
    for m = 1:length(Y)
       Y(m) = b(1) * X(m) + z(1);
       for i = 2:n
          z(i - 1) = b(i) * X(m) + z(i) - a(i) * Y(m);
       end
    end
    z = z(1:n - 1);
end


Object detection in binary image This task comes from tracking object on a steady background. So far I was able to remove the background and obtain binary masks like this:

I need to get bounding rectangle of the toy without the cord attached to it (on some frames it's quite thick). Sounds easy enough but I can't think of a good way to do it.

I tried cross-correlation with template of the toy I cropped. It wasn't really nice, plus I needed rectangle to be able to change it's size. And the cord added some mess.

meanShift() also didn't show good results because of the cord and occasional noice.

For now I stopped on findContours() with some tweaks, but it still produces some weird boundaries on what seems like okay image.

I keep thinking there must be some more simple and effective solution.

First I would recommend filling in the contour of the toy - in case it looks like the one in the second image. You could do this by analyzing the hierarchy output from  findContours: make  white all regions having a parent or by using an iterative morphological operations (not directly implemented in OpenCV).

Once the toy is nice and fat (not just the edge), you can remove the cord. I would use here some erosion with a flat kernel (horizontal line).

See the code below (this is in Python, but you'll have a clear overview):

import cv2

lemBGR = cv2.imread("lem.png")
lem = cv2.cvtColor(lemBGR,cv2.COLOR_BGR2GRAY)

# Dilate the image in order to close any external contour of the leming
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(5,5))
lem = cv2.dilate(lem,kernel)

# Identify holes in the leming contour
# This could be done by iterative morphological operations,
# but this is not directly implemented in OpenCV
contour,hier = cv2.findContours(lem.copy(),cv2.RETR_CCOMP,cv2.CHAIN_APPROX_SIMPLE)
# And fill them
for c,h in zip(contour, hier[0]):
    if h[3]!=-1:
        cv2.drawContours(lem,[c],0,255,-1)

# Now bring the leming back to its original size
lem = cv2.erode(lem,kernel)

# Remove the cord by wiping-out all vertical lines
kernel = cv2.getStructuringElement(cv2.MORPH_RECT,(15,1))
#lem = cv2.erode(lem,kernel) # first wipe-out
#lem = cv2.dilate(lem,kernel) # then bring back to original size
# erode and then dilate is the same as opening
lem = cv2.morphologyEx(lem,cv2.MORPH_OPEN,kernel)

# Find the contour of the leming
contour,_ = cv2.findContours(lem.copy(),cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)

# And draw it on the original image
for c in contour:
    # enter your filtering here
    x,y,w,h = cv2.boundingRect(c)
    cv2.rectangle(lemBGR,(x,y),(x+w,y+h),(0,255,0),2)

# Display the result
cv2.imshow("lem",lemBGR)
cv2.waitKey()

cv2.imwrite("lem-res.png",lemBGR)

And below are my results.


What is the maximum value that can result from a 2D DCT? Assuming that I have an 8x8 matrix of 8 bit unsigned values and carry out a 2D DCT on them, what is the biggest value that can result from it? This will help be decide how many bits to use. I am aware that the output coefficients are signed.

As suggested in the comments by Marcus Müller, you have to start from the equation of the DCT:
$$
S_{uv} = \frac{1}{4} C_u C_v \sum_{y=0}^7 \sum_{x=0}^7 s_{xy} \cos\frac{(2x+1)u\pi}{16} \cos\frac{(2y+1)v\pi}{16}
$$
The argument of the summation is of the magnitude of $s_{xy}$ times two values which are at most 1 (absolute value). So still 8 bits. You are summing $8 \times 8 = 64$ of those values, so $2^6 \times 2^8 = 2^{14}$. But you also divide this by $\frac{1}{4}$ so the maximum value requires 12 bits. You also have to consider that this value is obtained only when u anv v are 0, and in that case $C_u=C_v=\frac{1}{\sqrt{2}}$, so $C_uC_v=\frac{1}{2}$, that is another bit is unneeded. Some other checks are required when only $u=0$ or $v=0$, but overall you will need 11 bits at most.
This is exactly seen in the JPEG category tables for encoding the DCT coefficients.
For the skeptics, as you can see in the table below, when the input is in [0,255] the maximum values are always below 2048, thus 11 bits are enough. If your input is in [-255,255] you will need 12 bits.

0
1
2
3
4
5
6
7

0
2040
924,25
942,3571
924,25
1020
924,25
942,3571
924,25

1
924,25
837,4883
853,8957
837,4883
924,25
837,4883
853,8957
837,4883

2
942,3571
853,8957
870,6245
853,8957
942,3571
853,8957
870,6245
853,8957

3
924,25
837,4883
853,8957
837,4883
924,25
837,4883
853,8957
837,4883

4
1020
924,25
942,3571
924,25
1020
924,25
942,3571
924,25

5
924,25
837,4883
853,8957
837,4883
924,25
837,4883
853,8957
837,4883

6
942,3571
853,8957
870,6245
853,8957
942,3571
853,8957
870,6245
853,8957

7
924,25
837,4883
853,8957
837,4883
924,25
837,4883
853,8957
837,4883


How multi layer map to different beams in LTE Beamforming I understand for SU-MIMO Beamforming single layer maps to all antennas with phase shift applied to create beam. Also I understand in this case all beams will be carrying same signal (say s1).
But in case of MU-MIMO how can s1 and s2 be transmitted from multiple antennas at the same time so that one of the resulting beam carries s1 and other carries s2?

The answer is simply orthogonalization of radio resources. For any given OFDM symbol, exactly one user will be allocated certain sub carriers. This allocation is dynamic and informed to the user in what is known as the DCI, in the downlink control channel. Each user will decode data only on the subcarriers on which it has a grant allocated by the base station.

So each user is assigned a unique set of subcarriers for each OFDM symbol, hence when multiple users are multiplexed you do not get interference from other users. Because you are allocated a "unique orthogonal resource". Now if the base station wants to send multiple streams to the same user, then this is done by precoding matrices. 3GPP defines unique precoding matrices, the user estimates the channel and tells the base station about a probable precoding matrix and it is upto the base station to make the final call on the precoding matrix and inform it to the user.

Note: to your notion of beamforming on this question,

Beam forming, here is not in the sense of steering a beam in the direction of the user. But, precoding to transmit the stream in a given direction in the "vector space sense", where the vector space is formed by the N dimensional receiver space formed by the receive antennas. It is not a sense of physical direction as in conventional beamforming, it is more of beamforming in the receiver signal/vector space.


Echo Cancellation over Voice Telephone Network I'm a bit confused about the way echo canellation over a voice telephone network works like in the following figure

My issue is with the fact that the point of the adaptive filter and LMS algorithim is to drive the error e(t) to zero correct?
Then surely, Speaker B would hear nothing as the algorithim converges.

A usual arrangement of echo cancelers in a full-duplex system is shown in the following figure (from Digital Communication by Lee and Messerschmitt):

The echo cancelers try to remove any residual transmitted signal in the receiver of the same end. Since the input to the echo canceler is the transmitted signal, it can only remove signal components in the received signal that are correlated with the transmitted signal. It cannot remove any other signal components, such as the signal that is transmitted from the far end (which is the desired received signal). Consequently, the error signal can never become zero, and the desired signal (transmitted from the other end) will not be influenced by the adaptive filter.


Motion detection one threshold over dataset We're working on a project where we have to detect the moments where motions happen in a medical scene. I have a data set of 20 videos that are closely resemble (Same camera & configuration used, same shooting/environment conditions etc...). Here is what I did for each video:

Compute the optical flow (Farneback) between 2 consecutive images (t-1, t)
For each couple of images (t-1, t), compute the Norm1 of the displacement vector of all the pixels and then get the mean of these values. Then, I represent each couple of images by the mean of the Norm1 of optical flow results.
By doing so, I could draw a signal that represents the variation of these values as you can see below (Example video #1)

We zoom on the first 100 seconds of this signal to see the following:

To be able to detect the motions in video #1, we can just use as threshold the value 0.2 so we can differentiate between motion and not. But to be able to handle similar cases, we have to set a threshold at application/project level not on a video level: the threshold should be applicable over all the videos and for any new video of the same type. Let's take a look over the same kinda signal for video #2:

Obviously, the threshold in this case is 0.25 which is different from 0.2, the one we used for video #1 . To have a robust solution, we shall have one threshold per application so my question is about how can we use these thresholds (choosed manually) to have a global threshold for this application?

P.S: We already tried to use the mean of all these thresholds but it didn't work.. It failed to detect some important motions. We're asking such a question here because we just want to know if there is a fundamental principle to follow in such cases.

you should filter the data before you use the threshold. Run a low pass filter on the 1D signal (convolution with a Gaussian) and than normalize the data to the range between 0-1 after that i think you will be able to find a global threshold.


IFFT Synthesis from an empty spectrum I'm trying to generate a sawtooth wave using the IFFT. I can do this successfully by first performing an FFT and then putting that analysed spectrum into the IFFT, however I would really like to be able to write a UI that enables me to draw the bin magnitudes and generate a full spectrum ready for IFFT from that. I have been partially successful, however the waveforms don't look/sound quite right even though the spectrum input is per the formula for synthesising a sawtooth wave.

My question is how do I generate a spectrum array of complex values for an IFFT when all I have are the magnitude values?

I know that bin 0 and 1 are zero valued. I've been trying to render the wave by only calculating the magnitudes and phases for the significant bins (bins with magnitude > 0).
Am I able to put white noise in the remaining magnitude bins and then calculate the accompanying phases at the centre frequency or alternatively calculate these non significant bins  as random phase?  my start phase is -90, bin size 1024, complex fft.   Thanks in advance for your contribution.

I believe that if you will keep the phases and not convert it into magnitude (power spectrum) the IFFT will work. When you use the center frequency, the phase is 0, in other words, the real part is also the magnitude. As you move from the center frequency the error in magnitude growth as the phase growth. As result your error in recovery by means of IFFT growth.


Bilinear Transform (Tustin's Method) applied to the Derivative I hope that I have not misunderstood something terribly wrong, but the continuous derivative $D=d/dt$ can be considered a transfer function in Laplace space $D(s) = s$, right?

So when I try to discretize it using the bilinear transform (Tustin's method) I trivially get

$D(z) = \frac{2}{T} \frac{1-z^{-1}}{1+z^{-1}}$

When I apply this to a series containing one discrete impulse, the response oscillates at the Nyquist frequency. Even worse, the spectrum around $\omega=0$ is quadratic and not $\sim i\omega$ like it would be expected from the derivative. (Edit: the latter was just due to roundoff error, because the low-freq amplitude got swamped by the Nyquist-peak)

Although I know of course how well the bilinear transform works for discretizing all kinds of filters, it is hard for me to understand, why it is considered superior if it seems to fail so miserably for something as simple as the derivative, which can otherwise be easily represented by first order finite differences:

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

or even second order (symmetric) finite differences

$D(z)=\frac{1-z^{-2}}{2Tz^{-1}}$

I am sure it all has a very simple explanation, but I can't see it.

PS: What is all the more confusing: when I apply the bilinear D(z) to a step function, the result is (correctly) a single peak. Consequently the inverse of the bilinear D(z) applied to an impulse yields the step function, like it has to be. What is going on there?

It's a matter of perspective if the bilinear transform "fails so miserably" when trying to approximate a derivative. First of all, it's not true that the approximation is quadratic for small $\omega$, it's linear as it should be:

$$D(e^{j\omega})=\frac{2}{T}\frac{1-e^{-j\omega}}{1+e^{-j\omega}}=\frac{2j}{T}\tan\left(\frac{\omega}{2}\right)\approx \frac{j\omega}{T}\quad\textrm{for}\quad |\omega|\ll 1$$

So for small frequencies the discrete-time system behaves as a (discrete-time) differentiator. For larger frequencies we expect the effect of frequency warping, which is exactly what happens here: the infinite continuous-time frequency range is mapped to the interval $\omega\in[0,\pi]$.

In sum, the bilinear transform behaves as expected, but it is also clear that there are better ways to approximate a derivative in discrete-time than applying the bilinear transform to the continuous-time transfer function $H(s)=s$.


Downsampling an image by an integer factor When downsampling an image by an integer factor $n$, the obvious method is to set the pixels of the output image to the average of the corresponding $n \times n$ blocks in the input image.

I remember vaguely having read somewhere that this method is not optimal (sorry I don't remember any details)

Is it true that there is a better method (and if so, where does the above method fail, although it seems "obviously" correct)? I do not know a lot about signal processing, this question just interests me.

Downsampling an image reduces the number of samples that can represent the signal. In terms of frequency domain, when a signal is downsampled, the high-frequency portion of the signal will be aliased with the low-frequency portion. When applied to image processing, the desired outcome is to preserve only the low-frequency portion. In order to do this, the original image needs to be preprocessed (alias-filtered) to remove the high-frequency portion so that aliasing will not occur.

The optimal digital filter to remove the high-frequency portion (with the sharpest cutoff) is sinc function. The reason is that the Sinc function's frequency domain representation is a nearly constant 1 over the entire low-frequency region, and nearly constant 0 over the entire high-frequency region.

$$\text{sinc}(x)=\frac{\sin(\pi x)}{\pi x}$$

The impulse response of the sinc filter is infinite. Lanczos filter is a modified sinc filter which attenuates the sinc coefficients and truncates them once the values drop to insignificance.

However, being optimal in frequency domain does not imply being optimal in human eyes. There are upsampling and downsampling methods that do not obey linear transformations but produce better results than linear ones.

With regard to the statement about $n \times n$, it is important to keep in mind that during image sampling, the choice of coordinates correspondence between the high-resolution signal and the low-resolution signal is not arbitrary, nor is it sufficient to align them to the same origin (0) on the real or discrete number line.

The minimum requirement in the coordinates correspondence is that

Upsampling an image containing arbitrary random values by an integer factor, then downsampling by the same integer factor, should result in the same image with minimal change numerically.
Upsampling/downsampling an image consisting of just one uniform value, followed by the opposite operation, should result in an image consisting of the same value uniformly, with minimal numerical deviations.
Repeatedly applying pairs of upsampling/downsampling should minimize the shift in image content as much as possible.


Fourier transform of exponent: Delta pulse or hyperbola? Why do some tables say that Laplace (or Fourier?) inverse of exponential is a time-shifted delta pulse

\begin{align}
\delta (t) &\overset{\mathcal F}{\Longleftrightarrow} 1\\
\delta (t-t_0) &\overset{\mathcal F}{\Longleftrightarrow} e^{-j2\pi f t_0}\\
1 &\overset{\mathcal F}{\Longleftrightarrow} \delta (f)\\
e^{j2\pi f_0 t} &\overset{\mathcal F}{\Longleftrightarrow}\delta (f-f_0)
\end{align}

And table:
\begin{array}{c|c}
f(t)=\mathcal L^{-1}\left\{\mathcal F(s)\right\}&\mathcal F(s)=\mathcal L\left\{f(t)\right\}\\\hline
e^{at}&\displaystyle\frac{1}{s-a}
\end{array}

whereas others say that it is a (pole) hyperbola rather than a delta-pulse?

$$
e^{at}\overset{\mathcal F}{\Longleftrightarrow}\frac{1}{s-a}
$$

As pointed out by Peter K., it is important to distinguish between Laplace and Fourier transforms. The first few transform pairs in your question are Fourier transform pairs, whereas the last pair is a correspondence of the unilateral Laplace transform:

$$F(s)=\int_{0}^{\infty}f(t)e^{-st}dt$$

In the last transform pair in your question the $\mathcal{F}$ symbol is wrong because it would imply that it is a general Fourier transform pair. To be precise, the (Laplace transform) correspondence is actually

$$e^{at}u(t)\Longleftrightarrow \frac{1}{s-a}\tag{1}$$

because with the unilateral Laplace transform we only consider causal time functions, which satisfy $f(t)=0$ for $t<0$. Note that this is not necessarily the case with Fourier transforms.

Let's now consider the Fourier transform of the (causal) time function in (1):

$$F(j\omega)=\int_{-\infty}^{\infty}e^{at}u(t)e^{-j\omega t}dt=\int_{0}^{\infty}e^{-(j\omega - a)t}dt=\frac{1}{j\omega-a},\quad \textrm{for Re}\{a\}<0\tag{2}$$

Comparing (1) and (2) we see that, for $Re\{a\}<0$, the two transforms are indeed the same for $s=j\omega$. For $Re\{a\}>0$ the Fourier transform does not exist because the region of convergence of the Laplace transform $F(s)$ does not contain the imaginary axis.

Let us finally consider the time function $f(t)=e^{j\omega_0 t}$. If we were to multiply it with the step function $u(t)$, its unilateral Laplace transform exists according to equation (1) with $a=j\omega_0$. However, if we consider the function for $-\infty<t<\infty$, then its (bilateral) Laplace transform does not exist, whereas its Fourier transform does exist. It is given by a delta impulse in the frequency domain as shown in your table.

How can we now make sense of the relation between the Fourier and the (unilateral) Laplace transform? We must consider three cases:

The region of convergence of the Laplace transform $F(s)$ contains the $j\omega$ axis: then the Fourier transform is simply the Laplace transform evaluated at $s=j\omega$.
The region of convergence of the Laplace transform $F(s)$ does not contain the $j\omega$ axis: then the Fourier transform does not exist.
The region of convergence is $Re\{s\}>0$ but there are singularities on the $j\omega$ axis: both transforms exist but they have different forms. The Fourier transform has additional delta impulses. Consider the function $f(t)=e^{j\omega_0 t}u(t)$. From (1), its Laplace transform is given by

$$F(s)=\frac{1}{s-j\omega_0}$$
However, due to the singularity on the $j\omega$ axis, its Fourier transform is

$$F(j\omega)=\pi\delta(\omega-\omega_0)+\frac{1}{j\omega-j\omega_0}$$

It might seem that the Laplace transform is more general than the Fourier transform (when looking at the second point above), but this is actually not the case. In system theory, there are many important functions which are not causal, e.g. the impulse responses of ideal band-limiting (brick-wall) filters. For these functions the Laplace transform does not exist, but their Fourier transform exists. The same is of course true for sinusoidal functions defined for $-\infty<t<\infty$.


Uniform quantizer for gaussian input signal I have a gaussian signal $x$ with zero mean and unit variance. I am quantizing it by using a uniform quantizer of step size $q$ calculated by

$$q = \frac{x_\mathrm{max} - x_\mathrm{min}}{2^{M_\mathrm{ENOB}}},$$

where $M_\mathrm{ENOB}$ is the effective number of bits. Is this the correct way of quantizing a Gaussian signal?

The problem with quantizing Gaussian distributed signals (like the real/imaginary part of an OFDM signal) is that they can take any value in theory. It is thus necessary to clip such signals at threshold $C$ prior to quantization.

Low $C$ increases the distortion noise in this process, while large $C$ will lead to strong quantization noise. It can be shown, that there is an optimum value for $C$ that depends on the variance $\sigma^2$ of the signal and the quantizer resolution $M$. For details see [1,2]. After $C$ has been determined, the step size $q$ is calculated by
$$
q=\frac{2C}{2^M}
$$

By the way, I don't think it makes sense to use the ENOB (of a DAC) for calculating the step size. The ENOB is a kind of signal-to-noise ratio combining quantization noise and other effects like frequency selectivity and spurious frequencies of the DAC device. The number of input bits is still the nominal number of bits $M$.

Example: let $\sigma^2=1$ and $M=8$. Then with eq. (7) from [2] we can calculate the SNR $\gamma$ after clipping and quantization as a function of the clipping ratio $\zeta=C^2/\sigma^2$. The SNR should be maximized. I don't know a closed-form expression for the optimum $\zeta_\mathrm{opt}$ yielding $\gamma_\mathrm{max}$ but we can determine it numerically. For $M=8$ and with Fig. 2 from [2] we have $\zeta_\mathrm{opt} \approx 11.87\,\text{dB} \approx 15.38$. Furthermore, $C_\mathrm{opt}=\sqrt{\zeta_\mathrm{opt}\sigma^2}$ and with $\sigma^2=1$ we finally get $C_\mathrm{opt}\approx3.92$.

Remark: as you tagged this question digital-communictions you might be dealing with complex signals. In that case be careful with the definition of $\sigma^2$. In my answer it is the variance of a real-valued, Gaussian distributed signal. If the signal in question is complex-valued with independent and identically Gaussian distributed real and imaginary part, it has variance $2\sigma^2$, where $\sigma^2$ is the variance of the real or imaginary part.

References

[1] D. J. G. Mestdagh, P. Spruyt, and B. Biran, “Analysis of clipping effect in DMT-based ADSL systems,” in Communications, 1994. ICC’94, SUPERCOMM/ICC'94, Conference Record,'Serving Humanity Through Communications.'IEEE International Conference on, 1994, pp. 293–300.

[2] M. Bernhard, D. Rörich, T. Handte, and J. Speidel, “Analytical and numerical studies of quantization effects in coherent optical OFDM transmission with 100 Gbit/s and beyond,” in ITG Symposium on Photonic Networks, 2012, pp. 34–40. Online


How Is Laplacian of Gaussian (LoG) Implemented as Four 1D Convolutions? I am trying to understand the four 1D convolution operations involved in implementation of Laplacian of Gaussian(LoG). I have read this answer and I am also reading this pdf (See slide# 62 and 63). My current understanding is:

1) Pre-compute LoG and separate to 1D filters in x and y: gxx(x) and gyy(y).

2) Take Gaussian (g) and separate to: g(x) and g(y).

3) First apply g(y) and gyy(y) to the image. That makes two 1D convolutions.

4) Then apply g(x) to the result from gyy(y) and gxx(x) to the result from g(y). That adds two more 1D convolutions, making the total four.

Questions:

1) Is the above understanding correct?

2) How is this same as LoG?

3) Also, the diagram in slide #62 contradicts slide #63, in the pdf. Which one is correct?

Edit:

It's actually slide#61 and slide#62.

Edit2:

Please see the very last comment, posted Dec 19, 2018, on the first answer.

I think Chris Luengo's answer is perfect.
The trick is that you can calculate the 2nd derivative of the image (Using Finite Differences -> Convolution) and then blur it with Gaussian Filter.
Since both are seperable kernels you can do that by 4 1D convolutions.

If the input image is given by I

Calculate Ixx (The 2nd derivative on x direction) using convolution.
Filter Ixx with 1D Gaussian Kernel along the x direction.

Do the above for the y direction as well.
The LoG image is the sum of both.

To answer your questions:

See my explanation above. I think it matches your understanding.
This is the very definition (The Discrete) of the LoG.
I don't see the contradiction between #62 and #63.


Fixing for windowing function in FFT I've been learning about the FFT and about the usefulness of windowing functions. A problem I'm running into is the affect of windowing functions on the amplitude of signals in the FFT. For example, lets say I have a signal that is equal to:

y = 5.0 * sin(50.0 * 2.0 * pi * x) + 1.0 * sin(80.0 * 2.0 * pi * x)

And then I want to perform an FFT and receive the following plot:

Then if I perform a kaiser windowing function of the 14th degree, I receive the following FFT:

What I'm wanting to do is the following: Given a signal, I want to extract a wave out of it of a frequency I approximately know. To do this, I'd like to apply a windowing function to make it clearer what the frequency is I'm trying to extract but not lose information about the amplitude of wave I'm trying to extract. Is this possible? Like a way to reverse the amplitude decrease from the windowing function.

Perhaps the information at the following web page will be of some use to you:

http://www.dsprelated.com/showarticle/155.php


change the frequency component of a spectrum in an audio file I am working on an audio dynamic range compressor using Matlab.
After calculating the FFT of the audio file, I use several bandpass filters to amplify specific bands separately.
What I want to do then is remove components of a band and put them in another band.
for example, remove the components of the band 8khz-10khz and put them in the band 6khz-8khz.
would you please give me the general picture.

Judging by your question, it sounds like you are developing hearing aid solutions. Given that you are calculating an FFT and are now in the frequency domain, the easiest route forward is to use phase vocoder techniques to shift spectral peaks. Here is the original paper on this technique > http://505606.pbworks.com/f/LaroD99-pvoc.pdf


Determining Toy Helicopter Orientation From Image I am working on developing an autonomous control system for toy IR remote-controlled helicopters (like the Syma s107, for example) using a webcam and an Arduino with some IR LEDs. I am using C++ with OpenCV on a laptop running at about 2Ghz and I need to keep the system as close to real time as possible.

I am unsure of a way to effectively determine the rotation of the helicopter in the 2D image (more precisely its yaw with respect to the camera image plane). This information is important to being able to stabilize the helicopter in the screen because otherwise it floats out of view.

I have considered, and at least partially tried to determine the location of the body and the tail and use some tricks of geometry and perspective to get an estimation of shape but accurately finding the tail and body with low enough noise proved difficult. I have also considered using colored markers of some sort on certain spots of the helicopter but I'm not sure how to get a more accurate estimation than approximately left and approximately right.

For yaw it's sufficent to have three distinct data points fixed to the helicopter in a plane parallel to the rotor rotation: different colored leds for example. You could put them on three thin radial arms fixed to the landing gear. Maybe colored reflective materials would work, too, but you'd have to compensate for environmental light changes.

Then getting the yaw angle is simple.  Suppose the lights are red, green, and blue 120 degrees apart and ccw when viewed from the top.

Then in the image you will see these lights in various x-coordinate orders. As the helicopter rotates ccw fromt the position shown in the diagram:

G R B   // green and blue in front; red in back
R G B   // green in front; red and blue in back
R B G   // red and gree in front; blue in back
B R G   // etc...
B G R
G B R

(Of course at the borders of these regions two lights will merge to a single point.  These can be handled as special cases.)

So with 3 lights you have split all posible yaws into 6 x 60 degree segments.

You can further get down to 30 degree resolution by just comparing the left and right "gap" widths. If pointing directly away from the camera is azimuth theta zero degrees and the red light is on the nose (as in the picture), then the algoritm is:

Xr = -sin(theta)
Xg = -sin(theta + 60)
Xb = -sin(theta + 120)
if (Xg < Xr && Xr <= Xb) {
  // Green and blue are in front of red
  dxLeft = Xr - Xg
  dxRight = Xb - Xr
  yaw = (dxLeft < dxRight) ? 15+-15 : -15+-15;
}
else if (Xr < Xg && Xg < Xb) {
  // Green is in front of red and blue
  dxLeft = Xg - Xr
  dxRight = Xb - Xg
  yaw = (dxLeft < dxRight) ? 30+-15 : 60+-15;
}
else ... four more cases

The X+-15 above mean X degrees plus or minus 15: a 30 degree segment.  To get the exact offset, you can construct a preprocessed table that converts the ratio

R = dxLeft < dxRight ? xLeft / xRight : xRight / xLeft

into an an exact azimuth offset between -15 and 15.  This is pretty simple trigonometry to work out, and the nice thing is that it works independent of how far away the helicopter is, so long as it never turns upside down!

Another side benefit: the image distance between the outer two lights gives you inverse distance to the camera.


Can disparity range from negative to positive values for parallel cameras? In my stereo-rectified camera pair (not physical, but two positions of the same camera), when I project some 3D points in both images to find a rough disparity range, I often get a range that goes from negative to positive (e.g. -11, 24). I have no problem with ranges like (-22, -10) or (10, 15).

Is that kind of range even possible? When I do stereo matching how should I treat negative values? And what about the depth calculated from that disparity that can have negative values (and possibly areas with lowest cost at 0 disparity)?

It was the way I was setting principal points during rectification, to center the images back for the warping. If I have same principal points in both cameras I only get positive disparities, as expected. If cameras are parallel, it is not possible to have both negative and positive disparities. That only happens with converging stereo cameras (negative disparities for objects in front of point of convergence and positive for objects in the back).


Kalman Filter from 2D LookUp Table I have two sets of inputs, A (10 values) and B (20 values), and for each point (A,B) I have a measurement to make a (10x20) table of measurements C. Is there a way to use a Kalman Filter to improve these measurements if I have no physical rule to go from the inputs A and B to the value of C?

Is there a way to use a Kalman Filter to improve these measurements if I have no physical rule...

No.
Kalman filters add value by using a model of a process to refine an estimate of the states of that process.
If you don't have a model, and you don't have enough data to develop a model, then you can't get there from here.


Cut off frequency for low pass filter Hi I am trying to use a low pass filter in Matlab. I am trying to use the simple butterworth filter.

I have defined dataset say "data". Then I created the butterworth filter as follows:

[B,A]=butter(2,0.01)

Butterworth design require the user to give N and Wn. I used N as 2 and Wn as 0.01. Here, I randomly chose the value of Wn to 0.01.

My data is a hourly data and I want to apply 48 hour low pass filter. What would be the cut off frequency in my case ?

After that I used the following code to design a filter:

filter(B,A,data)

The filtered data looks good but I want to know what cutoff frequency should I use.

Thanks.

The butter command expects a normalized frequency from 0 to 1, where 1 is the Nyquist frequency.  In other words, Wn is in half-cycles per sample.

So you just divide your desired frequency (1 cycle / 48 hours) by half the sample rate ((1 sample per hour)/2), giving Wn = 0.04167 half-cycles per sample.


Amplitude distribution of Rician faded complex signal I have question regarding Rician fading simulation. I am trying to develop Rician fading channel simulator in LabVIEW, but I ran into some problems. I have complex sinusoid as my line-of-signs signal. Then I generated white gaussian noise samples for I and Q components which represent one scatered wave. I added this scattered wave to my LoS sinusoid. According to this I should get a faded signal with folowing amplitude distribution at different signal amplitudes:

When I calculate amplitude and its distribution with histogram I get following result:

I am using the same amplitude values and sigma value. There is significant difference in these results and my questions are:

What could be the problem here?
If I make multiple scatered waves, will it make difference?
Do I need to take into account Rician factor and somehow scale power of generated scatered waves?

This was a scaling problem. I did not have the same vectors for X-axis, thus creating uneven areas below PDF curves.


Forecasting noisy data with polynomials What's the best way to incorporate polynomials or sinusoidal equations into a linear model to predict a future amplitude in this sort of data? You can see that if those polynomials continued to go on, they wouldnt be anywhere close to fitting the final location. Im trying to use the greeen shaded area to predict the non shaded area.

Those light grey lines are polynomials like this: $$y_1 = a_1(x)^2 + b_1(x) + c_1\\$$ that are fitted to different sections of the data. Do I just use the a b and c variables from each polynomial, or do extrapolate out and use some other value?

Starting with your last question: to extrapolate you do just use your fitted coefficients and extend them to a greater $x$ range.

However, I think this approach is inherently flawed for data such as yours over any significant time period, especially when using a polynomial model. The problem is that any polynomial has a limited number of turning points. Once you pass these the data will tend to head off to $\pm \infty$ which is unlikely to model anything well.

Frequency domain analysis may be a possibility, but I expect you would run into similar problems. Another approach could be to use some sort of moving average or auto-regressive model.

Banks and similar financial institutions invest millions into modelling this sort of data so it definitely isn't easy.


Recovering signal from magnitude and phase of FFT I'm trying to code the basics of an analysis-synthesis system on MATLAB, but I'm getting incorrect results. From the Wikipedia page:

$$
X_k=\lvert X_k\rvert e^{i\angle X_k}
$$

Here's a simple MATLAB code:

y = speech(10:10+127);

yfft = fft(y,128);

spec = abs(yfft);
phase = unwrap(angle(yfft));

Y = spec .* exp(sqrt(-1)*phase);

yinv = ifft(Y,128);

yet y and yinv do not match. What's the issue here?

I tried running the same code, and it works fine, except for 2 minor points:

It's preferable to use exp(1i*phase) instead of exp(sqrt(-1)*phase)
The final result will deviate from the original signal due to rounding noise (typically around 1e-16). Also, some of this noise will appear in the imaginary part of the result (so, even if y was real, yinv will be complex, due to the small imaginary component). This will be particularly apparent if you try: plot(yinv), because Matlab will see that yinv is complex, and plot it on the real/imag axis


Can you use Bartlett method in reverse? I'm wanting to do an inverse Fourier transform.
Can I use the Welch method to generate this inverse, by replacing the FT used within with an IFT?

Not really. Welch's method is specifically designed to estimate spectral density and discards phase information.
IFT and FT are almost identical and have the exact same complexity. If you want to do an IFT, just do an IFT unless you have a really good reason not to.


Measuring microphone Signal to Noise ratio (SNR) with A-weighting Foreword: I have already read Applying A-weighting and user @endolith's  A-weighting code

I would like to measure the SNR of a microphone with the same method applied by the main manufacturers, that is displayed in their datasheets (example here with Neumann's U87 technical spec):

Here is how I plan to do it:

Play a continuous 1 kHz sine signal with a studio monitor loudspeaker (ideally in an anechoic room)
Adjust the playback volume such that the measured volume, close to the microphone, is 94 dB SPL, by using a sound level SPL meter device

Question here: should it be 94 dB SPL RMS or peak? (For a sine there's a difference of ~ 3 dB between RMS and peak)
Record the audio signal with a digital audio interface
Stop the playback of the sine and also record some silence (and thus some microphone noise: that is precisely what we want to measure!)
Normalize the recorded WAV file to 0 dB FS (so this will take the peak of the signal in consideration)

Is this step correct? It seems that a normalization is required here; if no normalization is done, the result SNR would depend on the preamp input gain!
Apply the A-weighting or CCIR-weighting (also called ITU-R 468 weighting), for example with endolith's code here
Compute the RMS of the A-weighted signal, during the silence/noise-only part.
This will be the final "Signal-to-noise ratio, A-weighted (re. 94 dB SPL) of the microphone".

Is this correct?

Here is a WAV file I recorded with the method described here. I'm curious to find how to find the A-weighted SNR of this microphone.

IEC 801-21-20 Defines sound pressure as the RMS of of the instantaneous pressure, and as far as I know all meters use RMS.  The problem you have in your implementation is that you are normalizing the peak value of the active region but measuring the RMS of the inactive part. As you described, for sine waves the RMS is 3dB less than the peak.  You could take the Difference between the RMS in the active to the inactive regions.  Or you could assume that the active RMS is -3dB and use that instead.  Otherwise you should be good to go.


Fourier-Analysis of Stationary Random Signals Let's say we have discrete-time stationary random signals with Gaussian PDF of mean value 0 and variance 1, whose individual signal values are uncorrelated.

For such a signal, how can we determine ACF and the PSD?

assuming ergodicity...

ACF:

$$\begin{align}
R_x[k] &= \lim_{N \to \infty} \tfrac{1}{2N+1} \sum\limits_{n=-N}^{N} x[n] \, x[n+k] \\
  &= \mathbb{E}\Big\{ x[n] \, x[n+k]  \Big\} \\
  &= \sigma_x^2 \ \delta[k] \\
  &= 1\  \delta[k] \\
\end{align}$$

PSD:

$$ S_x(e^{j\omega}) = \sum\limits_{n=-\infty}^{\infty} R_x[n] \, e^{-j \omega n} $$


Estimating sampling frequency offset I am searching a paper or article on this subject.
I have heard of a function that estimates the sampling frequency offset between nominal and actual (sample rate of samples) sampling frequency and corrects for the found value.
Particularly it is done by splitting the signal into blocks of block sized symbols and estimating the sampling time offset for each of the blocks.
Then, the sampling time offset is corrected for by cyclic time shift of the individual blocks.
I have searched all over but didn't find any.

As mentioned in a comment, what you describe looks like a timing recovery algorithm for digital communications.
I can recommend two sources, an article and a chapter in a book, which give a nice introduction and overview of the topic in general:

L. R. Litwin, «Matched filtering and timing recovery in digital receivers», RF Design, pp. 32–48, Sep. 2001.

M. Rice, «Chapter 8. Symbol Timing Synchronization», in Digital Communications: a Discrete-Time Approach, Upper Saddle River, N.J: Pearson/Prentice Hall, 2009.


How does taking the absolute value of a complex signal reflect in the frequency domain? I have a frequency-domain representation $X(e^{i\omega})$ of the complex discrete one-dimensional signal $x[n]$: $X(e^{i\omega})=\mathcal{F}\{x[n]\}$. Is there a frequency-domain transformation of $X(e^{i\omega})$ into $\hat{X}(e^{i\omega})=\mathcal{F}(|x[n]|)$?

Obviously, $\hat{X}(e^{i\omega})=\mathcal{F}\left\{\left|\mathcal{F}^{-1}\{X(e^{i\omega})\}\right|\right\}$ will do the trick, but I am interested if there is a dual to the absolute value function $|x|$ in the frequency domain.

Absolute value is not a linear operation, so no, there is not a consistent dual in the frequency plane.

Having said that, though, it is somewhat similar to calculating the signal power, which involves multiplying $x[n]$ with it's complex conjugate.  Multiplication in the time domain is equivalent to convolution in the frequency domain.  The result in the time domain is real, so the result in the frequency domain is even (i.e. symmetric around 0 Hz).

Absolute value will probably produce a similar effect to calculating the power since it is the square root of the power.  It is also real so it will be even in the frequency plane.


A-weighted Digital Filter for 1/3rd Octave Band Analysis I want to do a 1/3rd Octave band analysis on a noise signal. But the sound pressure levels should be in dB(A), i.e., and it should be A-weighted.
For this, I am using the A-weight function from this module.
From this link, what I understood is that the function A-weight takes the input as a signal(in Pressure(Pa)), and sampling frequency and gives the output as another signal(in Pressure(Pa)). If This output signal converted into pressure level should provide the A-weighted pressure levels.
Now, I use the A filtered sound pressure signal as my starting point and use the Butterworth filter to convert it into 1/3rd Octave Bands. Reference for this - https://stackoverflow.com/questions/56791652/1-3-octave-from-audio-file-with-python
The 1/3rd filter is as -

The output plot I get for the original signal(not weighted) is -

This is the correct output.
The output plot after applying the A-weights is:

This plot is not correct. The weights for each band should be:

Reference of this table
The method I am using is not giving correct results in the lower bands, and I cannot figure why.
I don't want to use the weights directly as I want a generally weighted signal for purposes other than just 1/3rd Octave Band Analysis too.
I am a beginner in DSP so that any help will be much appreciated. Feel free to suggest even alternative approaches for the A-weighting filter that can be executed in Python. Thank you


Detecting changes between voice and music I have a bunch of audio files that have both speech and music parts (think radio). For instance, a given file might consist of the following sections:

music --- speech --- music --- speech --- music --- speech --- music

My goal is to detect when the file changes between music and speech, and cut away the music parts leaving only the speech. However, I don't know what method to use to distinguish between speech and music. I was considering detecting the tempo, but some of the music is very slow or has a varying tempo.

(While we're at it, so that I don't have to write my own, are there any existing software libraries out there that can distinguish between speech and music?)

This is a well-studied problem, dating back from the mid 90s (DARPA/NIST broadcast transcription challenges). Search for "speech/music segmentation" or "audio segmentation" and you'll find thousands of research papers.

There are two broad approaches to solve this problem:

Supervised classification

Train a speech/music classifier, using a standard machine learning approach. You can use MFCCs as input features, along with other basic feature like zero-crossing rate, amplitude modulation at 4Hz, etc. Recently it became common to throw in as many features as possible, and using feature detection techniques to identify the most discriminant ones.

Any classification algorithm will do - support vector machines, gaussian mixture models, decision trees. Once the classification is done, you'll have misclassified frames (for example a tiny acapella segment in a song will be classified as speech; or a FX or jingle between speech will stand-out). This requires post-processing, the most common approach is to apply mode filtering (voting) on the sequence of classifier outputs. The classification/temporal smoothing are sometimes rolled into one through the use of hidden markov models for both classification and temporal smoothing.

Ref: Content-based audio classiﬁcation and segmentation by using support vector machines, Lu et al.

Unsupervised segment change detection

Consider a 10s window sliding over the signal. Compute audio features on the first half, on the second half, and use a statistical test to decide which hypothesis is the most likely: the two sets of audio features are drawn from the same distribution, or are drawn from two different distributions. The output of the test will tell you how likely it is that the middle of the window corresponds to the boundary between a speech and a music segment. Select the points with the highest scores as the segment boundaries.

The same audio features as for the supervised approach (MFCC, ZCR, amplitude modulation at 4 Hz...) can be used.

"Textbook" criterion for the statistical test: bayesian information criterion (BIC).

Ref: Unsupervised Audio Stream Segmentation And Clustering Via The Bayesian Information Criterion, Zhou & Hansen (for an introduction to BIC).

Combined Supervised and Unsupervised approaches for automatic segmentation of radiophonic audio streams, Richard, Ramona & Essid (for more exotic change detection tests).


Relation between the DTFT and CTFT in sampling- sample period isn't as the impulse train period I will quote an answer of Matt L. from this post
(I didn't comment there because I can't)

  If you have a continuous-time signal $x(t)$, then the two signals you're talking about are

  $$\begin{align}
x_c(t) &=x(t)\cdot\sum_{n=-\infty}^{\infty}\delta(t-nT) \\
 &=\sum_{n=-\infty}^{\infty}x(t)\delta(t-nT) \\
 &=\sum_{n=-\infty}^{\infty}x(nT)\delta(t-nT) \\ \tag{1}
\end{align}$$

  and you define

  $$x_d[n] \triangleq x(nT)\tag{2}$$

  The first signal given by $(1)$ is technically a continuous-time signal, even though it is only non-zero at discrete times $t=nT$. The reason why it is considered a continuous-time signal is because it can and must be transformed using the continuous-time Fourier transform (CTFT). So $(1)$ is the continuous-time representation of a sampled signal. Eq. $(2)$ is the discrete-time representation of the same signal. Here the sampled signal is represented as a sequence of numbers. You can't apply the CTFT to $(2)$, but you must use the discrete-time Fourier transform (DTFT).

  The nice thing is now that the CTFT of $x_c(t)$ given by $(1)$ and the DTFT of $x_d[n]$ given by $(2)$ are identical. So if the CTFT

  $$\begin{align}
X_c(j\Omega) &= \int_{-\infty}^{\infty}x_c(t)e^{-j\Omega t}dt \\
   &= \int_{-\infty}^{\infty}\sum_{n=-\infty}^{\infty}x(nT)\delta(t-nT) e^{-j\Omega t}dt \\
   &= \sum_{n=-\infty}^{\infty}x(nT) \int_{-\infty}^{\infty}\delta(t-nT) e^{-j\Omega t}dt \\
   &= \sum_{n=-\infty}^{\infty}x(nT)  e^{-j\Omega nT} \\
   &= \sum_{n=-\infty}^{\infty}x_d[n]  e^{-j n(\Omega T)} \\
\end{align}$$

  and the DTFT:

  $$X_d(e^{j\omega})=\sum_{n=-\infty}^{\infty}x_d[n]e^{-jn\omega}$$

  we have:

  $$X_d(e^{j\omega})=X_c\left(\tfrac{j\omega}{T}\right)\tag{3}$$

  In sum, the signals $(1)$ and $(2)$ are just two different representations of the same signal, and their spectra (one defined by the CTFT, the other defined by the DTFT) are identical.

My question arises due to a question I got when the period of the impulse train is not the same as the period of the sample:

$$x_c(t)=\sum_{n=-\infty}^{\infty}x(nT)\delta\big(t-n(\tfrac{2}{3}T)\big)\tag{1}$$

and

$$x_d[n]=x(nT)\tag{2}$$

does we still get the connection:

$$X_d(e^{j\omega})=X_c\left(\tfrac{j\omega}{T}\right)\tag{3}$$

if yes, does the $T$ in the connection imply to the impulse train period ($\frac{2}{3}T$) or to the sample time ($T$) and why?

When I set in  $X_d(e^{j\omega})$ the value $\omega=\frac{2}{3}T\omega$ I got an answer which should be true. but I don't know if it's by luck or really true and I don't know why.

We are given two signals:

$$x_c(t)=\sum_{n=-\infty}^{\infty}x(nT_1)\delta(t-nT_2)\tag{1}$$

and

$$x_d[n]=x(nT_1)\tag{2}$$

where $x(t)$ is some known signal whose Fourier transform $X(j\Omega)$ exists, and $T_1>0$ and $T_2>0$, but generally $T_1\neq T_2$.

We want to express the DTFT $X_d(e^{j\omega})$ of $x_d[n]$ in terms of the CTFT $X_c(j\Omega)$ of $x_c(t)$. The DTFT $X_d(e^{j\omega})$ is an aliased version of $X(j\Omega)$:

$$X_d(e^{j\omega})=\frac{1}{T_1}\sum_{k=-\infty}^{\infty}X\left[j\left(\frac{\omega-2\pi k}{T_1}\right)\right]\tag{3}$$

In order to compute the CTFT $X_c(j\Omega)$ of $x_c(t)$ we define a compressed ($T_1>T_2$) or expanded ($T_1<T_2$) version of $x(t)$:

$$\tilde{x}(t)=x\left(\frac{T_1}{T_2}t\right)\tag{4}$$

The CTFT of $\tilde{x}(t)$ is

$$\tilde{X}(j\Omega)=\frac{T_2}{T_1}X\left(j\frac{T_2}{T_1}\Omega\right)\tag{5}$$

With $(4)$ we have $x(nT_1)=\tilde{x}(nT_2)$, and $x_c(t)$ can be written as

$$x_c(t)=\sum_{n=-\infty}^{\infty}\tilde{x}(nT_2)\delta(t-nT_2)\tag{6}$$

The CTFT of $(6)$ is

$$X_c(j\Omega)=\frac{1}{T_2}\sum_{k=-\infty}^{\infty}\tilde{X}\left[j\left(\Omega-\frac{2\pi k}{T_2}\right)\right]\tag{7}$$

and with $(5)$ this becomes

$$X_c(j\Omega)=\frac{1}{T_1}\sum_{k=-\infty}^{\infty}X\left[j\left(\frac{T_2}{T_1}\Omega-\frac{2\pi k}{T_1}\right)\right]\tag{8}$$

Comparing $(3)$ and $(8)$ we finally obtain

$$X_d(e^{j\omega})=X_c\left(j\frac{\omega}{T_2}\right)\tag{9}$$


Unclear point about quantization Quantization of Continuous-Amplitude Signals

This point is not clear to me:

"The smallest quantization levels ($\pm\Delta$) correspond to the least significant bit of the binary code word."

Could someone explain it please?

Quantization is a process of representing real values of a signal using integer numbers. That process results in loss of accuracy, which is determined by the selection of the smallest value of quantization level.
If the smallest quantization level is $\Delta$, then every signal level will be rounded to the multiple of that $\Delta$. For example, if a signal level is $x=3,3 V$ and $\Delta = 1V$, then the corresponding quantized representation of the signal level will be $Q(x)=3\Delta$ and the resulting $\hat x = 3_{10} = 11_2$. If $\Delta = 0.5V$, then the signal level will be $Q(x)=7\Delta$ and the resulting $\hat x = 7_{10} = 111_2$.
So, to rephrase the original point: the weight of the least significant bit of the binary code word equals $\Delta$.


Can SIFT run in realtime? What is the best possible runtime for extracting SIFT keypoints and Descriptors?
I know it depends on #keypoints extracted. So, say for image of size 640x640, the code that I have been using requires around 3 seconds for description of ~2000 keypoints.

What is the best runtime achieved and how (on similar grounds)? How can I make it run in realtime?
Is hadoop with 100 clusters a solution?

SIFT implemented on a GPU can run in realtime!

http://cs.unc.edu/~ccwu/siftgpu/


Why is cepstrum used in MFCC instead of autocorrelation sequence? I know that the inverse Fourier transform of the power spectral density $P(k)$ gives the autocorrelation sequence $A(n)$ (By Wiener-Khinchin theorem):

$${\cal F}^{-1}\left\{P(k)\right\}=A(n)$$

And if we replace $P(k)$ by its log value, we get the cepstrum sequence $C(n)$:

$${\cal F}^{-1}\left\{\log_{10}P(k)\right\} = C(n)$$

My questions are:

What's the reasoning behind using log value?
Why does that make a difference?
Why isn't the autocorrelation sequence used in MFCC, LPCC, PLPCC, etc.. instead of the cepstrum?

Usually in speech signal processing cepstrum is used to represent low and high frequency components , which are multiplied with each other(in time domain its a slowly varying signal convolved with high frequency signal).

  What's the reasoning behind using log value?

logarithm is taken to remove the multiplicative effect and convert it to addition of two components. $$\log(a\cdot b) = \log(a) + \log(b)$$

  Why does that make a difference?

Now the two components are in addition form.

  Why isn't the autocorrelation sequence used in MFCC, LPCC, PLPCC, etc.. instead of the cepstrum?

Same reason, because we need to remove multiplicative effect of slowly varying component and fast varying component.

In speech signal processing slow varying component represents different phonemes and thus act as feature vectors for many applications(ex ASR). hence we need to remove the effect of fast varying component to isolate slow varying component.


Questions about MP3 bitrate conversions I was trying to analyse the lossy compression of MP3 audio files through their Spectrograms and the results were something like this (where the upper is 80Kbps and lower is 128Kbps).

It is obvious that the compression works by deleting higher frequencies from the audio. This is visible from the Spectrogram itself so I wanted to know

How exactly does the compression take place behind the simple ffmpeg commands?
Are spectrograms the best way to analyse the change in bitrates of an audio file?
The STFT used, gives us the phase and the amplitude from which we take the magnitude and plot the spectrogram. So how different can I expect the phases of the higher and lower bitrates to be?
Is it possible to somehow recover the lost frequencies because ffmpeg based bitrate upscaling doesn't produce accurate results on the Spectrogram

I am from a CSE background so I will be grateful if someone can answer these questions or atleast refer me to some article or paper regarding the same.

How exactly does the compression take place behind the simple ffmpeg commands?

MP3 is a perceptual coder. It creates a time/frequency representation and quantizes in the frequency domain primarily by calculating the spectral and temporal  masking thresholds of human auditory perception. It maximizes the masking of the quantization noise to minimize the audibility of the quantization noise.

Are spectrograms the best way to analyse the change in bitrates of an audio file?

No. Depends a bit on what you want to analyze, but a spectrogram doesn't correlate well with the audibility of compression artifacts. Most of problematic artifacts are temporal anyway, the loss of a few high frequencies is a comparatively smaller perceptual problem.

The STFT used, gives us the phase and the amplitude from which we take the magnitude and plot the spectrogram. So how different can I expect the phases of the higher and lower bitrates to be?

Phase will be fairly close for all spectral components that are above the masking threshold. Not so much for low energy components.

Is it possible to somehow recover the lost frequencies because ffmpeg based bitrate upscaling doesn't produce accurate results on the Spectrogram

No. If that would be useful in any way, the original inventors of MP3 (or AAC, DOlby, DTS, WMA, etc.) , would have included it in the codec. Many people have tried but it turns out the harm outweighs the benefit.


Pitch detection algorithms "give up" after very little I'm sorry if pitch detection has become a stale topic here, but I'm really interested if anyone has any advice on this topic. I'm looking for the best pitch detection algorithm for stringed instruments, but I think this question will turn into "What else needs to be done to most accurately and reliably detect a pitch of a stringed instrument?".

So basically I applied all algorithms from https://github.com/JorenSix/TarsosDSP/tree/master/src/core/be/tarsos/dsp/pitch and they all differ, but they all have one thing in common, on any device I try, the algorithm reports -1 Hz after around a second of a string on an instrument being played.

Basically the magnitude needs to be really high in order for the algorithm to detect anything. Now, obviously it's nothing on my part since I literally just copied the exact code and got those results. So there are obviously methods of still detecting a pitch in low magnitudes ( other applications, not sure what they use, can detect the pitch in the same situation for around 4-5 seconds which is a lot longer then these algorithms ), but I'm not sure what they are.
I'm not even talking about the fundamental being weak here, since the algorithm detects nothing.

Am I missing something, or is it the more likely option that something else needs to be done other than just applying a PDA to a sound buffer and if so, what is/are those things?

Thanks!

You might just have the magnitude or gain on your test input signals turned down too low.  Or the threshold set too high.  Or poor microphone placement.  Or perhaps you need an AGC controlled gain block before the pitch estimation algorithms to stabilize the decay portion of any ADSR-"like" amplitude envelopes or sound evolutions.

The above assumes you have a decent S/N ratio (e.g. little background noise and no harmony or accompaniment).

Also, you use a singular "best algorithm", when multiple algorithms (plural) may more appropriate for different portions of a note's envelope, and/or different sized stringed instruments.  In addition, your list of algorithms does not include any of the newer machine learning (ML or Convolutional DNN) inference methods of pitch detection/estimation or octave resolution.  Perhaps even RNNs to use melodic history to improve prediction statistics.


Detecting tennis courts in aerial images? Possible Duplicate:
  How to find tennis courts in aerial imagery

I am trying to detect tennis courts that are present in aerial images ..

  I have the following image:

  I tried color image segmentation on it but no success it gives a bad
  result :

What techniques should i apply for the same please suggest i think Hough transform or watershed could be used however watershed causes over segmentation and the results i get are not satisfactory

Tried using adaptive thresholding:

Some more things that i tried Connected components :

Output:

Since you are trying to discover structure in the images it's better to work in grayscale.

This is a very nice case where the court appears to be a nice rectangle, but in general, courts might come in different sizes and orientations. Also, white lines on green background is not a general rule, consider for example Roland Garos.

Having said this, you can try thresholding your image (you might have to use adaptive thresholding for this) to get something like this:

You can then go through the image creating chains of pixels that are 4 or 8 connected

Since you are looking for lines, you can get each of the chains and fit a line to it.

At this point, the "interesting pixels" in your image have been transformed to a set of line models and some of them will be intersecting. You can now apply simple rules to discover which sets of lines are intersecting (here is a good example ).

To discover the set of intersecting lines that make up a tennis court you would search for (this is an example) "A set of 9 lines where 4 are intersecting at right angles, two more parallel to the long lines intersecting them at 1/8 of the width of the field, two more parallel to the short sides and lying about 1/5 from the beginning of the court and a third running in the middle of the court and intersecting the two shorter ones"...In other words, you will have to express the tennis court as a set of intersecting lines with relative positioning to each other (so that your expression is free of the absolute dimensions).

For a very good application of the above technique you can have a look at this paper that outlines the recognition of shotcodes in an image.

For the theoretical background and more information about these techniques you can refer to the excellent "Digital Image Processing With Matlab" and specifically, chapters 10-13 which are the most relevant to your problem

I hope this helps.


Steady State of a Difference Equation using MATLAB I'm given this difference equation, and I'm being asked to plot its response to an input $x[n] = 4\sin^3\left(0.25\pi n+\frac \pi 3\right)u(n)$:

$$
y(n) = 0.9051y(n-1)-0.598y(n-2)+0.29y(n-3)-0.1958y(n-4)+0.207x(n)+0.413x(n-2)+0.207x(n-4)
$$

Using MATLAB's filter function, I implemented the following code:

num = [1 -0.9051 0.598 -0.29 0.1958];
den = [0.207 0 0.413 0 0.207];
x = 0:75; b1 = sin(0.25*pi*x+pi/3);
for i = 1:length(b1) sincub(i) = b1(i)*b1(i)*b1(i); end
inp = 4*sincub; z = filtic(num, den, [0 0 0 0]);
y = filter(num, den, inp, z);
stem(x, y, 'linewidth', 2);grid on
xlabel('Time-index, n'); ylabel('Amplitude'); title('Response to x(n)');

The plot of the input function is attached and so is the response.

The problem I'm facing is that I'm being asked to find the steady-state value of the output when the output appears rather periodic in nature (if x = 1:100 or 1:200 this will be clearer). Where is the transient response in such an output and what can be the steady state? Or is there something wrong in my code?

Any help would be much appreciated.

Given the trig identities you can rewrite your input as
$$x(t) = 3\cdot \sin(0.25\cdot \pi \cdot t+ \pi/3) - \sin(0.75\cdot \pi \cdot t+ \pi) $$

So it's just two sine wave going in, so it will be two sine waves coming out. You could potentially calculate the value of the transfer function $\omega = 0.25 \pi$ and $\omega = 0.75 \pi$ and manually calculate the output this way.