How can one improve the robustness of adaptive beamformers to signal mismatches? Compared to standard (conventional) beamformers, adaptive beamformers such as minimum variance distortionless response (MVDR) have much better resolution and interference nulling capabilities. Briefly, the MVDR beamformer is as follows:

If $\mathbf{a}_\theta$ is the look-vector in the direction $\theta$ and if $\widehat{\mathbf{R}}$ is the covariance matrix of observations, then the adaptive weight vectors are calculated as

$$\mathbf{w}_\theta=\frac{\widehat{\mathbf{R}}^{-1} \mathbf{a}_\theta}{\mathbf{a}^H_\theta \widehat{\mathbf{R}}^{-1} \mathbf{a}_\theta}$$

and the output power of the MVDR beamformer is

$$B_{mvdr}(\theta)=\mathbf{w}_\theta^H \widehat{\mathbf{R}}\mathbf{w}_\theta$$

For the sake of completeness, the conventional beamformer output power is $B_{conv}(\theta)=\mathbf{a}_\theta^H \widehat{\mathbf{R}}\mathbf{a}_\theta$.

However, these adaptive beamformers are also prone to severe performance hits when there is a signal mismatch i.e., when the steering vector model for the signal of interest doesn't match what's in the data or when the actual array positions are perturbed from the array positions used in the steering/look vectors (which is what happens in practice). The performance degradation can be so large that it becomes inferior to the conventional beamformer.

How can one increase the robustness of adaptive beamformers (in this case, the MVDR, but applicable more generally) to signal mismatches?

Example:

Here's the result from a simulation with two plane wave signals, one at $+30^\circ$ and the other at $-22.5^\circ$ and 10 dB down from the first, impinging on a 10 element uniform line array (spaced at $\lambda/2$, where $\lambda$ is the wavelength) in a noisy environment (the exact value of the noise power is immaterial here and you can assume it to be much lower than the signal powers). The covariance matrix is formed from 50 observations (snapshots). The top two panels show the outputs of the MVDR and the conventional beamformers and you can clearly see that the MVDR has high resolution compared to the conventional. Both do a good job of estimating the signal power.

However, when the element positions are slightly perturbed (perturbation is unknown to the beamformer), causing a mismatch in the actual signal steering vector and the model, the signal power estimates of the MVDR beamformer goes kaput. It can hardly measure them correctly, although it still does a good job of localising the signals. Comparatively, the conventional beamformer gets the relative powers right, although its harder to localise the signal due to increased sidelobes. (Note: The plots have all been normalized by the maximum power so although you can't see it, the peak at $30^\circ$ in the top left plot is at 0 dB. So this doesn't truly show how the actual power estimate degrades due to mismatches, as I just wanted to paint a qualitative picture.)

http://dsp.stackexchange.com/questions/146/how-can-one-improve-the-robustness-of-adaptive-beamformers-to-signal-mismatches#comment-308 :

  The title for the lower left plot should read MVDR: Perturbed array.
  I'll correct the mistake sometime soon.
  – yoda Sep 3 at 5:28
  -->

Is it possible for me to improve the estimates of the signal powers using the MVDR beamformer in the perturbed case?

One approach that I'm aware of is to "load" the diagonal of $\widehat{\mathbf{R}}$ by some value $\sigma_d^2$ in order to stabilize the covariance matrix and make the adaptive vectors more robust to signal mismatches. This is often done in matched field processing and is referred to as diagonal loading or white noise constraint beamforming.

In other words, the robustified adaptive weight vectors are calculated as

$$\widehat{\mathbf{w}}_\theta=\frac{\left(\widehat{\mathbf{R}}+\sigma_d^2\mathbf{I}\right)^{-1} \mathbf{a}_\theta}{\mathbf{a}^H_\theta \left(\widehat{\mathbf{R}}+\sigma_d^2\mathbf{I}\right)^{-1} \mathbf{a}_\theta}$$

and the output beam power is now computed as

$$\widehat{B}_{mvdr}(\theta)=\widehat{\mathbf{w}}_\theta^H \left(\widehat{\mathbf{R}}+\sigma_d^2\mathbf{I}\right)\widehat{\mathbf{w}}_\theta$$

Here's the result comparing the no mismatch, mismatch due to perturbed array and diagonally loaded cases for the MVDR beamformer for the situation in the question. There is an improvement in the signal power estimation, at the expense of increased background noise power.


Advice needed for separating 3-channel signals I'm working on a personal project, where I'm trying to make a touch screen that uses force sensors as the basis of the technology. The rig I've made has have three transducers underneath a square platform, with the goal of identifying touch events on the platform (with both the force and location) via triangulation.

I've got the simple "single event" case working fine, but I'm having trouble when there are two overlapping touch events. This would be the basis of a multi-touch function.

Below is a reference graph to help me explain what I'm trying to achieve:

We have three input signals (A, B & C) shown in green, with their summed signal in black. These signals contain two touch events, which overlap.

What I need to do is generate two "clean" 3-channel signals, each with only touch event inside.

The ideal output signals are shown with:
- the first "touch" event signals (A-1, B-1 & C-1) coloured in blue
- the second "touch" event signals (A-2, B-2 & C-2) coloured in red

I'm pretty sure this is possible, as I can at least identify when there are two overlapping signals; by using the algorithm below:

identifying the activity period between settles
scan over the activity period
measure the ratio between each signals gradient and the gradient of the sum each signal
If the ratio is off by more than a set threshold then there's more than one event

Does anyone have an idea on how to accomplish this or something similar?

Any advice or steers would be appreciated (especially in layman terms) as I'm self taught in this area.

Thanks a lot!

I understand that you receive A, B, and C and you want to separate it into A1, A2; B1, B2; C1, C2.

From this, you can do the following: Model your expected signal A1 as some pulse shape (here, I assume its a tanh function, which is close to your example). Then, normally you can perform matched filtering on the received signal to find peaks at the position of your two pulses. However, since your pulses are non-finite (i.e. they remain constant non-zero), I decided to use their derivative as the signal of interest. Then, you can find something like this:

n = -10:0.01:40;

noise = 0.5;

ntest = -5:0.01:5;
testpulse = tanh(ntest);
testpulse_diff = diff(testpulse);

% Amplitude of first pulse
A1 = -4;
start1 = 0;
% Amplitude of second pulse
A2 = 11;
start2 = 3;
% The input signal (sum + some noise)
sumOfBoth = A1*tanh(n-start1) + A2*tanh(n-start2) + sqrt(noise)*randn(1, length(n));

subplot(2,2,1);
plot(n, sumOfBoth)
title('sum signal');

subplot(2,2,2);
plot(testpulse_diff);
title('Diff of required pulse');

subplot(2,2,3);
plot(diff(sumOfBoth));
title('1st derivative of input');

subplot(2,2,4);
metric = xcorr(diff(sumOfBoth), testpulse_diff);
plot(metric);
hold on;
[peaks,locs] = findpeaks(abs(metric),'MinPeakHeight',0.01,'MinPeakDistance',100);

plot(locs, metric(locs), 'or');
hold off;

title('Metric with the detected peaks that correspond to the two signals');

Which gives you the following output:

From the two detected peaks you can derive both the amplitude and the time position of the both peaks. Note that, even though the first derivative of your input signal looks like pure noise, it contains all the information that is needed. This information is extracted by means of the cross-correlation between the expected and the received peak.


Noise added to a Random Process if we have a discrete random process
\begin{equation}
x\left(n\right)\:=\:0.2x\left(n-1\right)+w\left(n\right)+w\left(n-1\right)
\end{equation}
where $ w\left(n\right)$ is a noise with a mean
$
m_w=0.5$ and variance $\sigma ^2_w\:=\:1$.
I  am trying to find the value of $E\left\{x^2\left(n\right)\right\}$.
Here's my work below:
\begin{equation}
E\left\{x^2\left(n\right)\right\} = E\bigg\{\big[0.2x\left(n-1\right)+w\left(n\right)+w\left(n-1\right)\big]^2\bigg\}
\end{equation}
for
$$E\big\{w\left(n-i\right)^2\big\}=\sigma _w^2+m_w^2=1+0.5\:=\:1.5.$$
After that, I get the following
\begin{align}
E\left\{x^2\left(n\right)\right\} =&\phantom{+} 0.04 E\left\{x^2\left(n-1\right)\right\}\\
&+0.4E\left\{w\left(n\right)x\left(n-1\right)\right\}\\
&+0.4E\left\{w\left(n-1\right)x\left(n-1\right)\right\}\\
&+E\left\{w^2\left(n\right)\right\}\\
&+2E\left\{w\left(n\right)w\left(n-1\right)\right\}\\
&+E\left\{w^2\left(n-1\right)\right\}
\end{align}
Am not sure what to do from here? Please help.

Hi: I don't know whether white noise can have a non-zero mean ( Dilipe can tell us that )  but, if it can't, then just call it whatever you call white noise with a non-zero mean.
But, assuming it's zero, this is how you can do it. ( If it's not zero, then the  solution is similar but the last step would be a little different. Use
$$\operatorname{Var}\left\{X\right\} = E\left\{X^2\right\} - E\left\{X\right\}^2$$
Note also I wrote it with "t's" and "epsilon's"  just because I'm more used to those but you can think of $t$ as $n$ and $\epsilon$ as $w$.
\begin{align}
x_t &= 0.2 x_{t-1} + \epsilon_{t} + \epsilon_{t-1}\\
\implies x_{t}(1 - 0.2L) &= \epsilon_{t} + \epsilon_{t-1} \tag{$\scriptsize{\text{using lag operator}}$}\\
\implies x_{t} &=  \frac{\epsilon_{t}}{1 - 0.2L} + \frac{\epsilon_{t-1}}{1 - 0.2L}\\
\implies x_{t} &= \sum_{i = 0}^{\infty} 0.2^{i} \epsilon_{t-i} + \sum_{i = 0}^{\infty} 0.2^{i} \epsilon_{t-1-i}
\end{align}
So, since the mean is zero, the $E\left\{x^2_{t}\right\}$ is the variance of the RHS.
So, it is
\begin{align}
\sum_{i = 0}^{\infty} \left(0.2^{i}\right)^2 \cdot \sigma^2_\epsilon + \sum_{i = 0}^{\infty} \left(0.2^{i}\right)^2 \cdot \sigma^2_\epsilon&= \sum_{i = 0}^{\infty} 0.04^{i} \cdot \sigma^2_\epsilon + \sum_{i = 0}^{\infty} 0.04^{i} \cdot \sigma^2_\epsilon\\
 &=  \frac{2 \cdot \sigma^2_{\epsilon}}{(1- 0.04)}
\end{align}


Periodic signal filtering If I convolve a periodic signal $x(t)$ with period = 1sec with an aperiodic signal $h(t)$ whose Fourier transform $H(f)$ is exactly equal to 1 at frequencies $f = 0,\pm 1,\pm 2, \ldots \textrm{Hz}$, but has some arbitrary and finite values at non-integers, will the resulting waveform $y(t) = x(t) \star h(t)$ be the same as $x(t)$?

I think it should, but I want to validate my reasoning.

My reasoning is that since $x(t)$ is periodic with period=1sec, its spectrum is discrete and is non-zero at integer values.

So if $H(f)$ equals 1 at integers, it should not distort $x(t)$ regardless of what its values are in between integers.

Note that I am referring to frequency in terms of Hertz and not $\omega$ (angular frequency).

Is this reasoning valid?

The convolution in the time domain is translated to simple multiplication in the frequency domain. This means that since the values of the transform of x(t) are at the integer locations, the result of the multiplication will be the transform of x(t), as you defined the integer values for h(t) transform to be 1 at integer bins.


When should the sum of all elements of a gaussian kernel be zero? I found an approximation of a 5x5 2D convolution kernel like this :

Here, the sum of the elements is zero and this one was used for Laplacian of Gaussian!

Another one here :

This one has all positive and the sum is not zero and is being used for smoothing.

Could anyone please tell me what is with the sum being zero in the previous case ?

References:

Link1 and  Link2

The sum of a gaussian kernel cannot be zero, because all the elements are going to be positive. The first kernel you have shown, is most likely an edge detection kernel, (which is a type of high pass filter), so the elements add up to zero because you want to completely null out any DC/constant component.

The second kernel you have shown however, is a low pass kernel, indeed used in smoothing. Here the elements are not going to add up to zero because you are not trying to remove the 'constant' part of the image, in fact, you are trying to remove any high frequency noise.

Hope that helps.


CFAR Detectors and Interference Constant False Alarm Rate (CFAR) is an automatic thresholding technique. It's situated after the FFT processing in other to detect target peak. This probability of detection depends on clutter and noise level. In automotive RADAR, there are many other issues, like multi-path and interference. My questions are:

Is CFAR or XX-CFAR stay efficient detectors in presence of multipath and interference? If No, could you suggest me other techniques/methods to efficiently detect true targets in presence of interference and noise?

The "best" detector is a highly subjective subject, and there is (in my opinion) not a definitive answer here. I work in radar processing and I've used everything that you've mentioned in one form or another.  All of these CFAR methods are tools, and all of them have good/bad applications.  For instance, asking what the "best" wrench is in the toolbox highly depends on what you're going to use that wrench for.  This is no different.  OS CFAR can help with non-stationary noise floors which are varying quite a bit from CUT (cell under test) to CUT, but sometimes CA-CFAR can work just fine and is much simpler to implement.

In summary, they're all tools, and therefore no single CFAR detector is inherently "better" than any other one; it all depends on your system, the noise floor, the interference environment, the clutter environment, etc.


Envelope of a signal via Convex hull / Edge detection Is there a method in order to compute some sort of edge detection / convex hull of a signal x[n] (in orange here). Here of course the orange signal is not convex nor concave, but only piecewise...

(I would like to know if this could be an alternative way of computing the envelope of a signal, alternative to the traditional method input -> abs() -> low pass filter)

You can use the Hilbert transform to get the analytic signal.  The absolute value of the analytic signal is the envelope.  It's available in SciPy and many other software suites, though it's usually poorly named.


Does it make sense to recompute Parks-Mccellan on a per sample basis? I'm thinking of using Parks-Mccellan for my equalization needs, but since I require that the filter be variable (gain, Q etc.) then does it make sense to recompute Parks-Mccellan on a per sample basis?

Reason for asking is that I'm still confused about whether different algorithms are suited for per sample computation or whether they need longer intervals to make their full effect.

I guess it would still need some sort of pre-buffer (since the filtering uses previous samples)?

My application is Parks-Mccellan style arbitrary magnitude response filtering that's dynamic (i.e. I want to be able to vary the magnitude response over some, even short sample intervals).

No, the Parks-McClellan algorithm, which is based on the Remez exchange algorithm, is a computationally quite expensive iterative method that is most commonly used for off-line filter designs. It is not at all well-suited for real-time applications, let alone for re-designing a filter each sample interval.

Depending on your application, you might be looking for an adaptive filter based on the LMS (least mean squares) algorithm.

Another approach, again dependent on the application, is to use an appropriate parametrization of the filter, where the parameters could be the gain, center frequency, Q factor, etc. This is easily possible with second order IIR filters.


Equalize Speaker using Frequency Response Curve I have a speaker for which I know the frequency response curve (-5dB at 600Hz, -3dB at 700Hz, ... +5dB at 10kHz). This response curve is not as flat as I would like, and so I want to write a script which takes an audio file as input, and then transforms it such that when this new audio is played out through the speaker, it sounds like what the original file would sound like if it were played through a speaker with an actually flat response curve. (Basically, I want to flatten the frequency response of this speaker in software.)

What would be the best way of going about this? My first instinct is to Fourier transform the audio signal, reweight each Fourier coefficient according to the inverse of the amplitude at the corresponding frequency in the frequency response curve, and then inverse Fourier transform the result. However I imagine that such a direct approach would result in problems such as ringing artifacts, and I know that there are ways of performing similar tasks through the use of specially designed filters.

Oops I misread your question and thought you have a microphone recording of the speaker output with given audio, and answered that. Well here is anyway, as some parts still apply:

Like others have suggested in comments, designing a filter is the way to go. I would recommend a windowed FIR filter design rather than Parks-McClellan. From MATLAB documentation:

  firpm exhibits discontinuities at the head and tail of its impulse
  response due to this equiripple nature.

Our ears analyze sound with kind of wavelets so they may detect those discontinuities before and after impulsive sounds. With windowed filter design however you can choose a window that smoothly starts from and ends at value 0, avoiding creating such discontinuities.

Doing division (inverse of the amplitude) is a bit dangerous, because if you divide something by a number close to zero you get a huge value. Discrete Fourier transform results often contain values close to zero. If you have a long audio file then you can split the audio into shorter pieces with say 50 % overlap, multiply each with a window function, and calculate the mean squared magnitude of their discrete Fourier transforms (discarding phase). The square root of this is the root mean square spectral power, which is unlikely to be zero anywhere. Do this for both the speaker and raw audio and bin-wise divide the root-mean-square spectral power of the latter by the first. Still, you may get very large values if the speaker is unable to reproduce very low or very high frequencies. You can set a (soft) limit to how much correction you do at most, and you can even softly ramp the correction to nothing at frequencies you know the speaker simply cannot reproduce. Assuming zero phase for all bins, you can directly inverse discrete Fourier transform the division result to get an unwindowed zero-phase impulse response. You would then multiply the unwindowed impulse response with a window function which can be much shorter than the transform result, giving you a practical, short filter impulse response.

One more trick is to convert the filter to its minimum-phase counterpart, which will completely eliminate pre-ringing (taking place before the main peak in the impulse response). There are specialized methods for doing this that preserve the magnitude frequency response. I would not worry about post-ringing because real-world acoustics will give you plenty of post-ringing and echos anyhow and we are used to masking them.


Polar form of eigenvalue I greatly appreciate the opportunity afforded by this forum to submit a query.

It has been suggested that for a sinusoidal input, represented as 2 phasors, into a Linear Shift-invariant system, the output...

...can also be written as follows:

where $$H(\Omega) = \int_{-\infty}^{\infty} h(\lambda) e^{-j\Omega\lambda} d\lambda$$

Accompanying brief comment suggests the above to be on account of "polar form representation of H".

Although I am cognizant of polar form representation of complex numbers i.e.

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

...I am completely unable to fathom the alternative above.

All advice would be greatly appreciated.

Best regards,

wirefree

I think all it's saying is that we can write $H(\Omega)$ as:

$$
H(\Omega) = \left| H(\Omega) \right | e^{j \angle H(\Omega) }
$$

It also seems to be assuming that $h(t)$ is real-valued so that $H(-\Omega) = H^*(\Omega)$.

The function $H(\Omega)$ is just a map from $\mathbf{R} \mapsto \mathbf{C}$ (from the reals to the complex numbers).  For any specific value of $\Omega$ this generates a single complex number $H(\Omega)$. And, as you say, any complex number resolves to its exponential notation.


How to find inverse of convolution integral? If $x^{-1}(t)$ and $ y^{-1}(t)$ denote the integrals of x(t) and y(t)  defined by

$x^{-1}(t)=\displaystyle\int_{-\infty}^{t} x(\lambda)d\lambda$

$y^{-1}(t)=\displaystyle\int_{-\infty}^{t}y(\lambda)d\lambda$

then,

$(x*y)^{-1}$ which denotes the inverse integral of (x*y )is equal to:

(a) $x^{-1}*y^{-1}$

(b) $x*y$

(c)$x^{-1}*y$

(d)$x*y^{-1}$

i tried it follows:

from both conditions above we can easily infer following

$\dfrac{d}{dt}x^{-1}(t)=x(t)\implies sX^{-1}(s)=X(s)$

$\dfrac{d}{dt}y^{-1}(t)=y(t)\implies sY^{-1}(s)=Y(s)$

then i multiplied both(keeping in mind multiplication in freq. domain corresponds to convolution in time domain)

$s^2[X^{-1}(s)Y^{-1}(s)]=X(s).Y(s)$

$s[X^{-1}(s)Y^{-1}(s)]=\dfrac{X(s).Y(s)}{s}$

$[sX^{-1}(s)Y^{-1}(s)]=\dfrac{X(s).Y(s)}{s}$

$[X(s)Y^{-1}(s)]=\dfrac{X(s).Y(s)}{s}] \implies x*y^{-1}=(x*y)^{-1}$

and also,

if i would've written as

$s[X^{-1}(s)Y^{-1}(s)]=\dfrac{X(s).Y(s)}{s}$

$[X^{-1}(s)sY^{-1}(s)]=\dfrac{X(s).Y(s)}{s}$

$[X^{-1}(s)Y(s)]=\dfrac{X(s).Y(s)}{s}\implies x^{-1}*y=(x*y)^{-1}$

so, options (c,d) both should be correct but in the book the answer is given as

only option (c)

so , where i am doing mistake, and what should be the correct answer, i'm stuck

with this book's faulty answers...so,please help

The integral of a signal $x(t)$

$$f(t)=\int_{-\infty}^tx(\tau)d\tau$$

can be written as the convolution of $x(t)$ with the unit step function $u(t)$:

$$f(t)=(x\star u)(t)$$

Consequently, the integral of a convolution

$$g(t)=\int_{-\infty}^t(x\star y)(\tau)d\tau$$

can be written as

$$g(t)=((x\star y)\star u)(t)$$

And due to associativity of convolution we have

$$g(t)=((x\star u)\star y)(t)=((y\star u)\star x)(t)$$

So options $(c)$ and $(d)$ in your question are both correct.


The difference about WPCN and TDMA or FDMA I am a little bit confused about the difference WPCN and TDMA or FDMA,i mean the level difference.

WPCN is a wireless powered communication network,as the name implies,it can transmit the information without wire,so is it a transmitted technology or a transmitted paradigm?TDMA and FDMA are the method that transmitt information by dividing frequency and time.

What i really want to ask is that do them just like that i know a method which can transmit signal without wire,and i called this method "WPCN",and in this method,WPCN, i use TDMA or FDMA to transmit information,am i right?

TDMA and FDMA are paradigms; they provide different methods for sharing communications resources between multiple users.  WPCN is a technology for using wireless communications resources, and can be used with TDMA or FDMA depending on the particular application, though neither is "best suited" to the technology.


Create Matched Filter Reference Signal when Scaling Unknown I am working on a problem which raised what should be a simple question, but I'm stuck on it. The objective is to detect a shifted and scaled version of a known signal in the presence of Additive White Gaussian Noise (AWGN).

Think of a radar system. In a radar system, you know the transmitted signal. The returned system is delayed in time and scaled. I understand that a time-delay, only shifts the location of the matched filter's peak output. No problems there. But if the scaling is unknown, then how does one construct the reference signal?

Using a reference signal with the an amplitude that doesn't match the amplitude of the signal in the noisy data will mean that a threshold tests based on that reference signal will be incorrect.

I'm guessing that in a radar system, the noise power is approximately known, and so the power of the received signal is estimated by taking the power of the received signal and subtracting the noise power. That value could then be used to scale the reference signal. Still, that would only be and approximation, and it would still effect the probability of detection and false alarm in some way. Any help is appreciated. Thanks!

I think the problem is not as bad as you suspect it is.

I wasn't around at the time, but from what I've read, early radar systems essentially connected the matched filter's output to an oscilloscope, and a trained operator would look at the phosphor and decide, from experience and intuition, when the signal raised above the noise ("the grass") indicating a reflection.

In practice, you never know in advance the amplitude of the pulse you're looking for. There is attenuation in the channel that can change over time, even on very short time scales. Sometimes you don't know how much power the transmitter is using.

There are at least two solutions to this problem. One is to use automatic gain control in the receiver, so that you know in advance the power or the peak amplitude of the pulse going into the matched filter. The second is to use a training sequence or a pilot signal, so that the receiver can calibrate the thresholds appropriately.

To add a bit of detail: to simplify the math, assume the signal $p(t)$ is time-symmetric and has energy 1. The transmitter emits signal $kp(t)$, with $k$ unkown at the receiver. The receiver correlates the received signal with a matched filter with impuse response $lp(t)$. When the transmitted signal is present at the receiver, the matched filter's output is $$r=kl+n,$$ which is a Gaussian random variable with mean $kl$ and variance $N_0/2$.

When nothing but noise is present at the receiver, the matched filter's output is just $$r=n,$$ which is a Gaussian random variable with mean 0 and variance $N_0/2$. If the receiver can estimate $l$, then the problem reduces to standard hypothesis testing (see for example these notes).


How can I plot the frequency response on a bode diagram with Fast Fourier Transform? Assume that we have an unknow dynamical system and we only want to estimate its parameters. The system can be discribed as:

Continous time:
$$G(s) = \frac{3s + 5} {5s^2 + 3s + 2}$$

Discrete time with sampling time $h = 0.5$:

$$G_d(z) = \frac{0.36735z -0.15316}{z^2 -1.6551z + 0.74082}$$

My goal is to find the amplitudes of $H(z)$ and its frequencies by using  Fast Fourier Transform(FFT). Later I'm doing system identification to find the model $H(z)$ of the input and output data.

I first begin to do the frequency input response that changes over time.

$$u(t) = A(t)sin(2 \pi \omega(t) t)$$

Where $A(t)$ is the amplitude over time and $\omega(t)$ is the angular frequency rad/s over time.

Now I want to find the amplitudes and frequencies:

$$H(z) = \frac{Y(z)}{U(z)} = \frac{\mathcal{F}(y(t))}{\mathcal{F}(u(t))}$$

Question:

How can I do that?

I first begin with to create the transfer function and it's response $y(t)$

t = linspace(0.0, 10, 3000);
w = linspace(0.0, 100, 3000);
a = linspace(0.0, 10, 3000);
u = a.*sin(2*pi*w.*t);
G = tf([3 5], [5 3 2]);
Gd = c2d(G, 0.5);
y = lsim(Gd, u, t);

And the input $u(t)$

And how I try to do FFT on both $u(t)$ and $y(t)$

fy = fft(y);
fu = fft(u);
H = fy./fu;
plot(w, H);

But that does not work for me! Why?

Update:

I have made this MATLAB code.

% Data
close all
t = linspace(0.0, 10, 3000);
w = logspace(-1, 2*pi, 3000);
a = linspace(0.0, 10, 3000);
u = a.*sin(2*pi*w.*t);
G = tf([3 5], [5 3 2]);
y = lsim(G, u, t);

Ts = t(2)-t(1); % Sampling time
Fs = 1/Ts; % Sampling rate
fy = abs(fft(y, Fs));
fu = abs(fft(u, Fs));

% Cut the amplitudes and frequencies
freq = (0:Fs-1)(1:end/2+1);
fy = fy(1:end/2)/length(fy)*2;
fu = fu(1:end/2)/length(fu)*2;

freq(1) = freq(2); % freq(1) = 0
semilogx(freq, 20*log10(fy./fu));

And it will plot this "bode" diagram. I think it looks like it's reversed bode diagram. Should not look like this:

And if I use this input response with amplitude 5 and frecquency 10.

u = 5.*sin(2*pi*10.*t);

I get this plot. It seems that I'm on right track anyway.

Frequency Response of Unknown System from Freq Chirp and FFT's

My understanding from further discussion with the OP that he wants to specifically use an approach of providing a swept sine wave stimulus and use the FFT of this input and system output response to derive the transfer function. This may be for a system identification problem where the swept tone will provide more power per frequency bin than a direct impulse response can provide, so would be more practical for experimental test purposes. This would be an alternative approach for comparison to traditional channel estimation approaches using noise-like stimulus (PN codes) and least squares estimation techniques such as those described here: Compensating Loudspeaker frequency response in an audio signal

That said this approach should work well with certain considerations that I will outline below:

Frequency Ramp Generation This post and specifically the derivation from @MattL is a useful reference on setting the start frequency and stop frequency within a FM chirp (frequency ramp) signal to create the desired instantaneous frequencies accurately.

Simulation of a Frequency ramp

Here he provided the solution copied below for the values of $f(t)$ in the chirp  function $\cos\big(2\pi f(t) t\big)$ such that the instantaneous frequency will start at $F_1$ at time $t_1$ and end at $F_2$ at time $t_2$

$$\begin{align}f(t)&=F_1-\frac{\Delta F}{\Delta t}t_1+\frac12\frac{\Delta F}{\Delta t}t,\quad t_1<t<t_2\tag{1}\end{align}\\$$ with $\Delta F=F_2-F_1$ and $\Delta t=t_2-t_1$.

Deriving this for phase versus sample time $n$ including a constant frequency at the start and end results in the following formulas specific to DFT parameters with the resulting sinusoidal frequency ramp given as $\cos(\phi(n))$. Expressing in units of phase versus time ensures phase continuity at the transitions:

$\phi(n)=\begin{cases}\omega_1n,&n \le n_1\\ \omega_1 n+\frac{\Delta \omega}{\Delta n}\big(\frac{n^2+n_1^2}{2}-n_1n\big),&n_1<n<n_2\\\omega_2n+(\omega_1-\omega_2)n_2+\frac{\Delta \omega}{\Delta n}\big(\frac{n_2^2+n_1^2}{2}-n_1n_2\big),&n\ge n_2\end{cases}\tag{2}\\$

With the resulting frequency ramp as $cos(\phi(n)$), in the time interval $[n = 0,N-1]$

Where:
$n_1$ : index for start of ramp with frequency $\omega_1$ radians/sample
$n_2$ : index for end of ramp with frequency $\omega_2$ radians/sample

The parameters for the frequency ramp are further detailed in the graphic below:

Windowing Windowing will be important to minimize distortion in the DFT results. However given we are sweeping the input frequency with time, tapering the signal at the boundaries will reduce the signal levels at these test frequencies. An excellent window choice for this application is the Tukey window as we can selectively taper just the outer edges, while the majority of the window is flat, offering significant performance in frequency even with a relatively small $\alpha$ which is the ratio of the taper portion of the window to the flat portion.  Additionally with a low $\alpha$ the resolution bandwidth is minimally impacted.

Sampling Rate Given the application will be to use a real tone, the sampling rate needs to be higher than twice the highest frequency over which we would like to measure the transfer function.

Number of Samples The number of samples is set based on the resolution bandwidth desired for the transfer function measurement. The number of samples will set the total time duration $T$ of the test signal, which will set the resolution bandwidth of the test according to $1/T$ (as mentioned above, the Tukey window if low $\alpha$ is used will not significantly impact the resolution bandwidth.)

Transfer Function With the above considerations, the transfer function is derived by the ratio of the output FFT to the input FFT. This would be a complex function with its associated magnitude and phase components.

Demonstration

This was done using the same transfer function as the OP used in his own answer  for direct comparison to his results (as well as those responses to this question Why does this transfer function estimation not work? System identification) instead of the original one given in the question, repeated here as:

$$G(s) = \frac{3}{s^2 + 0.5s +30}$$

In application the transfer function is the unknown, and the objective is to determine the frequency response of this unknown transfer function, and specifically in this case using a real sinusoidal frequency ramp stimulus and FFT's. So this transfer function was used to generate the actual output which was then used along with only the input to predict the frequency response.

Below is the resulting estimated frequency response versus ideal showing excellent agreement between the two:

Here is the MATLAB/Octave code for the optimized sinusoidal (real) ramp for FFT processing. This generates a cosine frequency chirp with Tukey amplitude taper at the start and end of the chip, along with constant frequency in the taper for optimized FFT performance for use with determining a transfer function.

function out = chirp(N, r = 0.05, w1=0, w2=pi);
  # Dan Boschen 4/18/2020

  n= 0:N-1;                       # sample index
  n1 = ceil(r*N/4);               # start ramp at 50% in window rise
  n2 = N-n1-1;                    # end ramp at 50% in window fall

  # phase versus time for linear ramp
  psweep = w1*n + (w2-w1)/(n2-n1)*((n.^2+n1^2)/2-n1*n);

  # constant frequency outside of Tukey window
  psweep(1:n1) = w1*n(1:n1);
  psweep(n2:end) = psweep(n2)+w2*n(n2:end)-w2*n(n2);

  win = tukeywin(N, r)';        # Tukey Window
  out = win.*cos(psweep);

end

The following script demonstrates proper use of the chirp:

N = 4096;             # number of samples
fs = 60;              # sampling rate
n = 0:N-1;            # sample index
t = n/fs;             # time index

# frequency chirp stimulus
u = chirp(N, r = 0.05, w1=0, w2=pi);

G = tf([3], [1 0.5 30]);    # OP's trasfer function
y = lsim(G, u, t)';         # Generate output as OP has done
fu = fft(u);
fy = fft(y);
resp = fy./fu;              # derived frequency response

# plot
figure;
faxis = n/N*fs;
half = floor(N/2);
semilogx(2*pi*faxis(1:half), 20*log10(abs(resp(1:half))));

The reason this works so well is due to the excellent flatness in the DFT of the chirp signal through careful planning in generating such a signal for use in the DFT.  Following the considerations listed above, the chirp stimulus has the following characteristics versus frequency and time. The taper duration is greatly exagerated in this diagram below; a 5% ratio was actually used for the demonstration above.  Importantly at no time during the chirp does the frequency repeat; given the requirement for a real signal, any repetition of a frequency later in time would result in deep nulls in the response. Going below $0$ or above $\pi$ would essentially create such a repetition, therefore the frequency was made to be constant at the edges while allowing the ramp to extend fully from $0$ to $\pi$ in normalized radian frequency (cylces/sample). Sidelobe ringing was minimized by extending the ramp 50% into the taper of the amplitude window.

An FFT of the chirp signal itself shows these performance features including complete frequency coverage and overall flatness approaching that of an ideal impulse (while allowing for much larger overall signal power thus we would expect much better performance versus an impulse response test in the presence of noise.)


Difficulties while understanding Convolution I have just started convolution sum. I am able to do Mathematical convolution sum. But I am unable to make a figure of it in my mind. So I tried to go through the graphical convolution.
Below is the picture

I understand decomposed part but I don't understand that how each decomposed part is generating x[ 0 ]h[ n ],x[ 1 ]h[ n-1 ]....?

I guess a little help makes me to understand the Convolution sum.
Thanks

The convolution sum is
\begin{equation}
y[n] = \sum_k x[k]h[n-k]
\end{equation}
If you take the term corresponding to $k=0$, the contribution is $x[0]h[n]$. Similarly the contribution of $k=1$ is $x[1]h[n-1]$. Here, $h[n]$ is delayed by $n=1$. Similarly you can see how each term of the sum contributes to overall sum. The sum is taken over values $k$ where $x[k] \neq 0$ because only those terms contribute to the sum.


How does integration of a signal change the spectrum? It's a very basic question, but I havn't found the question on google or dsp.stackexchange.com:
How does integraion of a signal change the spectrum?

It is the basic property of Fourier Transforms. Given that $g(t)$ and $G(f)$ are Fourier transform pairs i.e.
$g(t)  \rightleftharpoons G(f)$
then
$\int_{-\infty}^tg(\tau)d\tau \rightleftharpoons \frac{1}{j2\pi f}G(f)$

This assumes that $G(0)=0$. If $G(0)$ is not zeros then the integral of $g(t)$ has a Fourier transform that includes a Dirac delta function.


Cross-Correlation to find the similartity of songs I have been working on a problem where I am given two songs which are of the same category but there is a time difference between them. Eg. one song starts at x sec and other song starts at y sec. So my task is to find the x and y.
The length of the two songs can vary but the content of the two songs will remain the same.

Currently I am using cross-correlation using FFT for finding the offset between them but it is not giving me correct results.

For making the length of signals same I am zero padding the signal which is shorter in length.

Can someone guide me what I am doing wrong in this problem ?

I really appreciate your responses.
Thanks

--

Thanks guys for replying..
I am doing in Python and using below algorithm for finding the offset:
1. Using librosa.load function to get the audio signals of both the songs.
2. Then I zero pad the shorter signal from last and make sure it is power of 2.
3. Find the FFT of the two signals using scipy library function of python.(spipy.fftpack.fft)
4. Once I get the fft , I take the cojugate of one of the signals.
5. Elementwise multiplication of the two fft signals
6. Doing ifft of the numpy array of the output signals after multiplication.
7. Finding the peak index of that array.

P.S. The song length can be around 5-6 mins so I am sampling at the rate of 44.1Khz . So you can expect very large sample space.
Also , the bitrate of the two songs can be different.(I am not sure that will effect it or not)

Any idea if I am doing wrong?

thanks

First of all, apply a normalized cross correlation, which is more robust to certain variations. MATLAB's xcorr function has an option to do that. Next, apply it in a sliding window fashion, where you apply the correlation as a running filter. Then, simply select the maximum peak.


What is the frequency of the Gaussian white noise model in control systems? The design of control systems, particularly for SISO systems, is made convenient by tools such as Bode plots and the corresponding Nyquist stability criterion. These concepts allow engineers to design sufficient margins for the given system in a model based paradigm, when it's transfer function is known. Specifically in the Bode plot, the frequencies correspond to sinusoidal input signals and the plot provides information about the gain and phase margin available for stability. The idea is that by knowing the worst case of frequencies one can design the controller to reject disturbances.
However, in the case of stochastic noise, say the commonly used model of Gaussian white noise, it is not clear to me what the corresponding frequency representing this type of noise would be on the Bode plot. What is the frequency of the Gaussian white noise model in control systems? The first variation (time derivative) of this noise is infinite indicating very large frequency. Additionally, does this fact point to the evolution of stochastic control as a separate sub-discipline within robust control (control methods to reject adverse effects of disturbances and noise)? Finally, is there a frequency-domain representation of Gaussian white noise?

No such thing as a single frequency of the noise. That's exactly why it's called white; it has power in all frequency ranges, but not at a single frequency.

Finally, is there a frequency-domain representation of Gaussian white noise?

Yes, a constant power spectral density for all frequencies. That's like white light (which contains also a continuum of all frequencies), which is where it got it name from.
Dealing with continuous spectra, like that of white noise, but also those of actual signals is absolutely bread and butter of control engineering – we very rarely have an input to a system that's always constant or just a sinusoid, which means that its spectrum is never the actual one tone you mention in your question.
I'm not quite sure how to proceed here – I can't really paste three thirds of "Otto Föllinger – Regelungstechnik" (my German standard control engineering textbook) as an answer, and I'm pretty sure that would also not be super helpful to you – because probably, you'll touch upon how you can work with signals and systems in Laplace or Fourier domain, where this really becomes easier.
I'll thus venture a quick "outlook":
Your control system is linear (nod here, please, if you agree, and if you're not certain what linearity is, look it up), thus, if you input a sum of two different signals, the output is equal to the sum of outputs you would have gotten if you used them individually.
Thus, you can start approximating your white noise. Noise has a measurable power. So, pick one arbitrary frequency and put a tone of the same power here. A single line in frequency domain is not a good approximation of a noise that has no lines but is a constant "floor", but at least we've got the same power.
To make the approximation better, you now split the power (i.e. half it, which means multiplying the amplitude with $\frac1{\sqrt2}$) and make two tones that divide the frequency span you're looking at evenly. You calculate the output of your system for both inputs, add it, and get a possible output.
Small problem: You need to assume some phase of your tones. Well, noise can (and will) have independent random phase at different frequencies, so you calculate the output power for a large set of uniformly chosen random phases. You calculate the average power and claim it's an estimate for the output power frequency distribution for white noise.
Instead of halving the power and giving it to two tones, you could of course also had divided the power by 1,000,000 and given it to 1,000,000 tones. Same principle applies. Your estimate of the output power distribution will probably be much better than with 2 strong tones, since 1,000,000 weak tones with random phases are much more alike white noise than 2. However, your math senses are tingling intensely, and you recognize that you're basically doing a large $\sum\limits_f$ over 1,000,000 frequencies and really, you need to do it over uncountably many frequencies, and that sounds like a job for an $\int\cdots\,\mathrm df$.
Thus, you end up with an integral transform, such as the Laplace or Fourier transform. With these, you can describe signals that are deterministic, but have a continuous spectrum, which your Bode plot method seemingly cannot! Think about it: What is the "frequency" of "Input is 0 for eternity, until it turns 6 for a duration of 2, after which it's 0 again for eternity"? There is no single frequency in aperiodic things!
(and by the way, that's why I wrote it has no power at a single frequency: you divide a finite amount of power by an uncountably infinite amount of frequencies. Every frequency gets exactly 0 power. If you integrate over a range of frequencies, even an extremely small one, you get some power, so the result of your division is some kind of power density!)
The fact that the input is stochastic means you're in trouble if you assume fixed phase relationships, and that's why we usually have to resort to power spectral densities. However, from a simple "how does the output power look like if I feed white noise into a low-pass RC filter with a bode plot that looks like ...", you'll see that this is a very useful info to have!
Don't be afraid: The math you've learned for single-frequency inputs is absolutely necessary to build the continuous spectrum understanding, and it will be used all over the place in your continued studies of control engineering (or in fact, all disciplines of electrical engineering, if that's what you're studying). More often than not, the problem really boils down to "well, then just multiply the Bode plot point-wise with the input spectrum and sum / integrate things up"; and that's not really hard to do.


Why are these results different about 'frequency shifting'? I think 'spec_ori' & 'spec_shifted' spectrum amplitude and phase are same in below code. Because of 'spectrum shiting property' (What effect does a delay in the time domain have in the frequency domain?).
But both are different. why?

K = 4 * 10^10;
Fs = 10 * 10^6;
t_len = 40 * 10^(-6);
t=-t_len/2:1/Fs:t_len/2; % time axis in range
omega = -Fs/2:1/t_len:Fs/2;
chirp = exp(1i*pi*K*t.^2);
spec_ori = fft(chirp);
spec_shifted = fft((chirp .* exp(1i*2*pi*Fs*t)));
plot(angle(spec_ori)-angle(spec_shifted))

These angular errors are in the order of $10^{-10}$. These are numerical inaccuracies caused by the fact that exp() can't be arbitrarily exact.

"Exact" would mean that exp(1i*2*pi*Fs*t) would just yield a vector [1 1 1 … 1]; if you'd hit that, the input to the fft would be identical, because multiplication with $1$ doesn't do anything.

So, always consider the error relative to the phenomenon you're observing. You wondered about something that was about 200 dB lower in amplitude than the observed quantity.


FFT on continuous fifo data - Is y axis magnitude sum of all samples? I'm running an fft once a second on a buffer of data 60 seconds long. The data is sampled at 558Hz and is placed into the fifo buffer once per second. FFT is numpy.rfft. The data is scaled to psia prior to fft. From the research I have done so far (found several previous questions on this) the unit of the fft magnitude will also be in psia. But I am unsure of how I should interpret this magnitude. Is it a sum of the entire sample? Something else? My intuition says it is a sum of the entire sample. Where if I had say a peak at 25Hz of 30psia on the 60 second buffer, if I wanted to know the peak over a 1 second sample I would divide by 60. My sensor data is only varying by about 2-3 psi when plotted in time domain so that got me wondering how the fft y-data could be orders of magnitude larger than the pk-pk difference of input signal. Or maybe I'm thinking about this totally wrong, I am pretty ingorant on this : )

In general the  output of a Fourier Transform is "spectral density". For a continuous Fourier Transform that's implicit: if the signal in Volts, the unit of the FT are Volts/Hz. For the discrete case, the units are the same but the frequency interval is implicit: it's the sample rate divided by the FFT length.
If you want to extract "physical" meaning from the FFT, you need to watch your scaling. In your specific case it would probably be best to scale both the forward and backward FFT by $1/\sqrt{N}$ since that maintains Parseval's theorem, i.e.
$$\sum |x[n]|^2 = \sum |X[k]|^2$$
which means that total power in the time and frequency domain are the same.


IIR with asymptotic impulse response I have plot some of the IIR filter's impulse responses and some of them are exponentially increasing, increasing, decreasing etc but some IIR's impulse response approaches 0. If for example I have an IIR filter defined by its lccde:

$$ 6 y[n] + y[n-1] - 2y[n-2] = x[n] - 2x[n-1] $$

or

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

when I plot its impulse response in matlab:

[h,n] = impz([1 -2],[6 1 -2]);
stem(n,h);

its impulse response would approach 0, and this zero value will bestarting at n = 23. Does this mean that the IIR is an FIR since I can cut off the impulse response at n = 23?

As per the exact solution of a discrete LTI system:

$$y(t)=C A^k x_0+C\sum_{k=0}^{t-1}A^{t-k-1} B u(k)+D u(t)$$

and considering that $u(t)=\delta(t)$ and $x_0=0$ for the impulse response, we have:

$$h(t)=C A^{t-1} B+D \delta(t)$$

where $A$,$B$,$C$ and $D$ are the state space discrete matrices:

$$x(t+1)=Ax(t)+Bu(t)\\y(t)=Cx(t)+Du(t)$$

Which in this case are obtained as:

b=[1/6 -1/3 0];
a=[1 1/6 -1/3];
[A,B,C,D]=tf2ss(b,a);

Hence, aside numerical issues, both methods gave the same results:

[h,t] = impz([1 -2],[6 1 -2]);
stem(t,h);

t=(0:23)';
h(1)=D*1;
for k=2:length(t)
    h(k,1)=C*A^(k-2)*B;
end
stem(t,h);

So the $t$=23th term on the $h(t)$ signal is:
$$h(23)=C A^{t-1} B|_{t=23}=C A^{22} B=-3.39702260269318e^{-05} $$

Though small, the $t$th term never vanishes, hence, transfer functions with denominator different than the identity vector, will always be an IIR.


How to resample in MATLAB one signal that is supposed to be convolve with another one? i have two audio signals .wav and i want to convolve them, the trouble is that the frequency do not fit, they aren´t just the same frecuency-lenght in the wavread arguments. i am supposed to use the lesser frequency to resample the other one using interp1 and splines. I really appreciate your help.
Thanks

If the ratio of the sample rates is reasonably rational (i.e can be represented as the quotient of reasonably small integers) then you can use the MATLAB function resample().

Convolving to audio signals is rather unusual. Typically and audio signal is convolved with filter impulse response. The impulse response could be stored in an aduio file, but that would be a bit unusual. I would double check if you understood the task correctly.

Both interp1() and spline() are poor choices for sample rate conversion as they don't control the mirror images in any predictable way. I recommend reading up on "sample rate conversion".


Linearity and shift-invariance of 2-D system on lattice I know how to check a 1-D system for these conditions, but am confused about translating this to a 2-D system over a lattice. The system $H$ is define as:

$$w[\mathbf x] = H\left\{u[\mathbf x]\right\} = \sum_{\mathbf y  \in \mathbf Q} \varphi[\mathbf y] u[\mathbf x - \mathbf y],\quad\text{with}\quad\begin{cases}
&\mathbf x \in \Lambda \subset R^2,\\
&\mathbf Q = \left\{\mathbf q \in  \Lambda : \lVert\mathbf q\rVert \le \rho\right\},\\ &\quad \rho \in R,  \ 0 < \rho < \infty
\end{cases}$$

I have mainly done problems with very simple 1-D systems. I know that the summation part represents convolution and that the check for linearity is:

$$H\left\{af_{1}[\mathbf x] + bf_{2}[\mathbf x]\right\} = aH\left\{f_{1}[\mathbf x]\right\} + bH\left\{f_{2}[\mathbf x]\right\}$$

and the check for shift invariance is:

$$w[\mathbf x - \mathbf y] = H\left\{f[\mathbf x-\mathbf y]\right\}$$

However I can't figure out how to actually do these checks on this particular system.

Well, actually you only need to insert the system $H$ into the definition of linearity and check, if it holds:

$$\begin{align}H\{af_1[x]+bf_2[x]\}&=\sum_y\phi[y](af_1[x-y]+bf_2[x-y])\\
&= a\sum_y\phi[y]f_1[x-y]+b\sum_y\phi[y]f_2[x-y]\\
&=aH\{f_1[x]\}+bH\{f_2[x]\}\end{align}$$

So, the linearity holds. There is nothing special about the fact that the system is 2-dimensional or similar. You just follow the rules for the calculation, and you end up with the statement of linearity.

A similar thing can be carried out for the shift-invariance.


Which step response matches the system transfer function A system has the following open loop bode plot: -

Which one of the plots below describe the closed loop step response for the entire system?

My attempt
My initial thought was to look at the static gain from the open loop bode plot, which is $-4\text{dB} = 0.63$. I interpret this as this will be the final value of the step response when it is settled. The only graph that sort of matches this description is graph 3.
However, it turns out that this is not the correct answer. And therefore my reasoning about the static gain can not be correct either. So my question is, how can I identify the closed loop step response from the open loop bode plot?

The final value of the step response is the DC gain of the closed-loop transfer function, which is generally different from the open-loop DC gain.
Assuming unity gain feedback, the feed-forward transfer function $G(s)$ equals the open-loop transfer function, and the closed-loop transfer function is given by
$$C(s)=\frac{G(s)}{1+G(s)}\tag{1}$$
The final value of the closed-loop step response equals $C(0)$. With the given value of $G(0)$, which I would interpret as $-3$dB, the final value of the closed-loop step response should be around $\sqrt{2}-1\approx 0.4$. The step response in Figure $1$ seems like a good match.


Adding Poisson noise to a time series signal I'm working on a radiation detection problem. I have simulated some measurement count data and want to add an amount of Poisson background noise to achieve a certain SNR. Radiation count measurement (due to the background radiation) is given by the Poisson random variable.

After searching the web, I'm not seeing any clear indication of how to go about it doing this. Is there a similar method as in the case of AWGN? (https://www.mathworks.com/help/comm/ref/awgn.html).

If your data $\mathbf{x}$ has variance $\text{var}(\mathbf{x})$ and you want to add Poisson noise to it with signal to noise ratio of $\text{SNR}$ dB, then you can just generate Poisson random numbers with appropriate variance. The variance of the noise $\mathbf{n}$ you need is $\text{var}(\mathbf{n})=\frac{\text{var}(\mathbf{x})}{10^{\text{SNR}/10}}$. So now you just need to use your Poisson random number generator and then add the noise to get the noisy data: $\mathbf{x}+\mathbf{n}$


End of chirp in phase 0 I would like the chirp to end in phase zero.
Chirp time or end frequency may vary slightly.
Now I'm checking the maplot output.
I observe the same thing at the output of the sound card.
In this solution, the chirp does not always end in phase zero.
help find a solution.
Thanks for the help.
Andrew
from pylab import *
import numpy as np
from scipy.signal import chirp

f0 = 7000
f1 = 17000
samplerate = 192000
T = 0.0013
T = np.ceil(T*f1)/f1 # new T
t = np.arange(0, int(T*samplerate)) / samplerate
w = chirp(t, f0=f0, f1=f1, t1=T,phi=270, method='linear')

fig, ax = subplots(figsize=(6,1))
ax.set_title("Chirp ")
ax.plot(w)
show()

Derivations here - you can pick any tmin, tmax, and fmin, fmax for some number of samples or sampling rate N.
We adjust the code one line toward the end to rescale phi to end as an integer multiple of $2\pi$ to yield zero phase; this has the effect of nudging fmin and fmax, slightly or greatly depending on all other parameters - see here.
An alternative variant that forces only the end of the chirp to be zero phase will exactly preserve fmin and fmax, by subtracting. Forcing all, zero phase at tmin and tmax without changing fmin and fmax is impossible.
Taking your parameters, with cosine and sine:

Code
import numpy as np
import matplotlib.pyplot as plt

def _lchirp(N, tmin=0, tmax=1, fmin=0, fmax=None):
    fmax = fmax if fmax is not None else N / 2
    t = np.linspace(tmin, tmax, N, endpoint=True)

    a = (fmin - fmax) / (tmin - tmax)
    b = (fmin*tmax - fmax*tmin) / (tmax - tmin)

    phi = (a/2)*(t**2 - tmin**2) + b*(t - tmin)
    phi *= (2*np.pi)
    return phi

def lchirp(N, tmin=0, tmax=1, fmin=0, fmax=None, zero_phase_tmin=True, cos=True):
    phi = _lchirp(N, tmin, tmax, fmin, fmax)
    if zero_phase_tmin:
        phi *= ( (phi[-1] - phi[-1] % (2*np.pi)) / phi[-1] )
    else:
        phi -= (phi[-1] % (2*np.pi))
    fn = np.cos if cos else np.sin
    return fn(phi)

#%%######################################################################
f0 = 7000
f1 = 17000
samplerate = 192000
T = .0013

N = int(samplerate * T)
tmin = 0
tmax = T

t = np.linspace(tmin, tmax, N, endpoint=True)
for zero_phase_min in (True, False):
    for cos in (True, False):
        x = lchirp(N=int(samplerate * T), tmin=tmin, tmax=tmax, fmin=f0, fmax=f1,
                   zero_phase_tmin=zero_phase_min, cos=cos)
        plt.plot(t, x)
        plt.title("cos={}, zero_phase_tmin={}".format(cos, zero_phase_min),
                  weight='bold', fontsize=17, loc='left')
        plt.show()


DFT exercise in the book Understanding digital signal processing 3 Ed I am trying to solve exercises from the book Understanding digital signal processing 3 Ed - Richard Lyons.
I will repeat the question as it is in the book:

3.3 We want to calculate an N-point DFT from a one-second compact disc (CD) audio signal of duration $x(n)$, whose sample rate is $f_s = 44.1 \text{kHz}$, with a sample spacing of $1 \text{Hz}$.
(a) What is the number of time samples $x(n)$ required, N?
(b) How long is the sequence time $x(n)$ measured in seconds? Tip: This part (b) of the problem is more complicated than it looks at first. Think carefully.

My answer to question (a) is $N = 44100$, since the DFT frequency resolution is $\frac{f_s}{N}$. So, for $1 \text{Hz}$ resolution, I think it's quite obvious.
But, in question (b), I'm little insecure. As the hint suggests, looks like a trickier question, but I don't know exactly why.
My first guess is that the duration of the sequence is 1 second.
Would anyone know why it wouldn't be?
Thanks!

@Henrique Luna. Please forgive me. In Part (b) of the problem the words "sequence time" should "time sequence". Sorry for the confusion!
Years ago when I created that Part (b) question I was thinking about the answer to the question; "What is the time duration of an N-sample time-domain sequence?" Back then I believed the time duration under consideration to be the time difference between the sequence's first and last samples, rather than N times the time interval (1/fs) between the samples. Thinking this way leads to the conclusion that Part (b)'s answer is (N-1)*1/fs = 0.999977 seconds.
But subsequent discussion, years ago, with Dilip Sarwate convinced me that my concluded answer to 3.3 Part (b) may be somewhat valid when thinking about the time-domain but my concluded answer leads to confusion when thinking about signal durations in the context of spectrum analysis.
I now believe the 3.3 Part (b) question leads to more confusion than enlightenment. That Part (b) question was not my most well thought-out homework problem.
By the way Henrique, if you have not done so already, please have a look at the following web page: https://www.dsprelated.com/showarticle/1094.php


Additive noise meaning What is additive noise? What actually additive means?
I have tried to search in the internet but the only answer that I get is the noises gets added therefore it is additive, which make me think that noises are not destructive in nature am I right?

noises gets added therefore it is additive

correct!

which make me think that noises are not destructive

incorrect :(
simple thought experiment: You flip a fair coin $X$ (Head = -1 / Tail = 1) and tell me the result. The entropy here is 1 bit, i.e. the (expected) information ($I(X=\xi) = -\log_2 \left[P(X=\xi)\right]$) of each outcome is 1 bit.
Then there's additive noise $N$ that takes one of the values $\{-2,0,+2\}$ with equal probability.
When you receive a -1, you can't know whether the coin was Head and there was 0 noise, or the coin was Tail and there was -2 noise. Both are equally likely!¹
So, your additive noise is absolutely able to destroy information and hence is very destructive to your signal.
If you're more coming from a wireless communications background: Your $X\in\{-1,+1\}$ can be interpreted as BPSK. Now you see how even benign Gaussian noise destroys your reception when its sign is the opposite of your transmit symbol!

¹ we can even formalize that. Since $X$ (2 options) and $N$ (3 options) are independent, and each of them equidistributed, there's six possible combinations, each of them equally likely
 X | N  | Y = X+N
------------------
-1 | -2 | -3
-1 |  0 | -1
-1 | +2 | +1
+1 | -2 | -1
+1 |  0 | +1
+1 | +2 | +3

Thus, we have four possible outcomes for the sum of signal and additive noise, -3, -1, +1 and +3.

If we see +3 or -3, we get 1 bit of the 1 bit info in the coin toss. (It must have been +1, otherwise we couldn't get +3, or -1 for -3, respectively.) That happens in 2 out of 6 times, so with probability 1/3.
If we see -1, we don't know whether it was +1-2 or -1+2, so we have zero bit of the 1 bit of the coin toss. Same for +1. That happens 4 out of 6 times, so with probability of 2/3.

Thus, the expected information to get out of this channel is 1/3·1 + 2/3·0 bit = 1/3 bit, where put in full 1 bit! That's a very destructive additive noise channel.


Fourier transform 4 times = original function (from Bracewell book) I was glancing through "The Fourier Transform & Its Applications" by Ronald N. Bracewell, which is a good intro book on Fourier Transforms. In it, he says that if you take the Fourier transform of a function 4 times, you get back the original function, i.e.
$$\mathcal F\Bigg\{ \mathcal F\bigg\{  \mathcal F\big\{  \mathcal F\left\{  g(x) \right\} \big\} \bigg\} \Bigg\} = g(x)\,. $$
Could someone kindly show me how this is possible? I'm assuming the above statement is for complex $x$, and this has something to do with $i^0=1$, $i^1=i$, $i^2=-1$, $i^3 = -i$, $i^4=1$?
Thank you for your enlightenment.

I'll use the non-unitary Fourier transform (but this is not important, it's just a preference):
$$X(\omega)=\int_{-\infty}^{\infty}x(t)e^{-i\omega t}dt\tag{1}$$
$$x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(\omega)e^{i\omega t}d\omega\tag{2}$$
where (1) is the Fourier transform, and (2) is the inverse Fourier transform.
Now if you formally take the Fourier transform of $X(\omega)$ you get
$$\mathcal{F}\{X(\omega)\}=\mathcal{F}^2\{x(t)\}=\int_{-\infty}^{\infty}X(\omega)e^{-i\omega t}d\omega\tag{3}$$
Comparing (3) with (2) we have
$$\mathcal{F}^2\{x(t)\}=2\pi x(-t)\tag{4}$$
So the Fourier transform equals an inverse Fourier transform with a sign change of the independent variable (apart from a scale factor due to the use of the non-unitary Fourier transform).
Since the Fourier transform of $x(-t)$ equals $X(-\omega)$, the Fourier transform of (4) is
$$\mathcal{F}^3\{x(t)\}=2\pi X(-\omega)\tag{5}$$
And, by an argument similar to the one used in (3) and (4), the Fourier transform of $X(-\omega)$ equals $2\pi x(t)$. So we obtain for the Fourier transform of (5)
$$\mathcal{F}^4\{x(t)\}=2\pi\mathcal{F}\{X(-\omega)\}=(2\pi)^2x(t)\tag{6}$$
which is the desired result. Note that the factor $(2\pi)^2$ in (6) is a consequence of using the non-unitary Fourier transform. If you use the unitary Fourier transform (where both the transform and its inverse get a factor $1/\sqrt{2\pi}$) this factor would disappear.
In sum, apart from irrelevant constant factors, you get
$$\bbox[#f8f1ea, 0.6em, border: 0.15em solid #fd8105]{x(t)\overset{\mathcal{F}}{\Longrightarrow} X(\omega)\overset{\mathcal{F}}{\Longrightarrow} x(-t)\overset{\mathcal{F}}{\Longrightarrow} X(-\omega)\overset{\mathcal{F}}{\Longrightarrow} x(t)}$$


Recursive version of DFT as presented in Cooley-Tukey paper The seminal paper of Cooley and Tukey  provides an iterative method for computing the DFT for a sequence of length $N$. Specifically, they mention a method which utilizes the fact that $N$ can be written as $N = r_1 r_2$. The two intermediate results during computations themselves are DFT computations, therefore, their computation in turn can be recursively estimated.

Is there a reference/method of recursive formulation for computing DFT using the $N=r_1 r_2$ factorization ?  The versions of FFT I come across use factorizations where $r_1 = \frac{N}{2}$ and $r_2=2$. I would like to know if there are recursive formulations for general factorization of $N$.

The Cooley Tukey method allows general factorizations $N=N_1N_2$. The easiest way to see this is to use index mapping. Consider the length-$N$ DFT of a sequence $x[n]$:

$$X[k]=\sum_{n=0}^{N-1}x[n]W_N^{nk}\tag{1}$$

with $W_N=e^{-j2\pi/N}$. Let's assume that $N$ can be factored as $N=N_1N_2$. Now use the following index mapping

$$n=n_1N_2+n_2,\quad n_1=0,1,\ldots, N_1-1,\quad n_2=0,1,\ldots,N_2-1\\
k=k_1+N_1k_2,\quad k_1=0,1,\ldots, N_1-1,\quad k_2=0,1,\ldots,N_2-1\tag{2}$$

Plugging (2) into (1) gives

$$X[k_1+N_1k_2]=\sum_{n_2=0}^{N_2-1}\sum_{n_1=0}^{N_1-1}x[n_1N_2+n_2]W_N^{(n_1N_2+n_2)(k_1+k_2N_1)}\tag{3}$$

The complex exponentials in (3) can be simplified as follows:

$$W_N^{(n_1N_2+n_2)(k_1+k_2N_1)}=W_N^{N_2n_1k_1}W_N^{N_1N_2n_1k_2}W_N^{n_2k_1}W_N^{N_1n_2k_2}=\\=
W_{N_1}^{n_1k_1}\cdot 1\cdot W_{N}^{n_2k_1}W_{N_2}^{n_2k_2}\tag{4}$$

because $W_N^{N_1}=W_{N_2}$, $W_N^{N_2}=W_{N_1}$, and $W_N^{N_1N_2}=W_N^N=1$. Plugging (4) back into (3) and rearranging terms gives

$$X[k_1+N_1k_2]=\sum_{n_2=0}^{N_2-1}\left[W_N^{n_2k_1}\left(\sum_{n_1=0}^{N_1-1}x[n_1N_2+n_2]W_{N_1}^{n_1k_1}\right)\right]W_{N_2}^{n_2k_2}\tag{5}$$

Equation (5) can be interpreted as follows: the result can be computed by computing $N_2$ length-$N_1$ DFTs of subsampled versions of the input sequence $x[n]$ (this is the inner sum in (5)), then multiplying the result by twiddle factors $W_N^{n_2k_1}$, and finally computing $N_1$ length-$N_2$ DFTs (the outer sum in (5)). This procedure is most easily visualized using matrices:

read the data row-wise into a $N_1\times N_2$ matrix
compute DFTs over all $N_2$ columns (this can be done in-place)
multiply the result element-wise with twiddle factors
compute DFTs over all $N_1$ rows (in-place)
read out the data column-wise to get the final length-$N$ result

This little octave/Matlab script illustrates the procedure:

N = 1536;          % non-power of 2 FFT length used in LTE
N1 = 512; N2 = 3;  % 512 x 3 = 1536
x = randn(N,1);    % some input signal

% generate N1 x N2 matrix of input data, reading them row wise
X = reshape(x,[N2,N1]).';

% DFT over each column
X = fft(X);

% Multiplication with twiddle factors
X = X.*exp(-1i*2*pi/N*(0:N1-1)'*(0:N2-1));

% DFT over each row
X = (fft(X.')).';

% read out column wise
X = X(:);


simulator for dsp kit tms320c6713?without using actual hardware kit tms320c6713? As you know, most universities across the globe are closed due to pandemic  and people are unable to acess research facilites ,so in such a scenario one wants to improve his dsp hardware programming skills, how can he make up for dsp kit tms320c6713? I mean, is there any software proteus/matlab or Labview or anyother software that enables to simulate tms320c6713 without using actual hardware tms320c6713

Usually, even before someone builds a processor, they write a specification of the instruction set, and then they build a simulator for that.
The same applies to TI's C67xx series: https://processors.wiki.ti.com/index.php/List_of_Simulator
The wiki article strangely says that TI is moving away from simulation in Code Composer Studio 6, which would honestly be a terrible, terrible, terrible idea and put them back in the 1980s with respect to software engineering practices (no chance to do unit tests on any machine).
The only reason for doing so I could think of is lack of manpower to maintain that simulator, and/or a plan to deprecate the C67xx series altogether.


Coarse summary of various methods for detecting when a signal 'starts' My question is basically in regards to how to detect when (in time) a signal 'starts', and it might also be called 'transient detection', or, as I have read in some papers regarding music, the 'onset' of a signal.

The motivation here is ultimately a Time-Delay-Of-Arrival (TDOA) estimate between two signals, but suppose you do not have the luxury of doing a cross-correlation because the signals are either way too long, and/or there is so much multi-path corruption that the cross-correlation method doesnt work. In this case my belief is to look for where a signal 'starts' and go from there.

I am looking for more breadth than depth here, as in, a survey of possible/common methods used in doing such a thing.

Thanks!

If both signals are audio signals from the same source (just different paths to get there) you can do

Calculate the envelope (various methods for that, depends on the signals and the accuracy you want)
Down sample
Do a running cross correlation (at the down sampled rate)


Inverse polyphonic chord recognition Polyphonic music transcription does not currently appear to be a solved problem.

How about the inverse of a small portion of the problem.  Are there any kind of spectral characteristics (from an STFT) that can be used to eliminate some musical chords from the probability space?  (e.g. this snippet of sound most likely does not contain any C# chord, or any kind of diminished minor chord, or this is a single note not a chord, etc.)

Assume the audio snippet is more-or-less stationary (transient attack removed, etc.), and that overtones for most or all individual notes are very likely present.  (And this question is not about inverted chords.)

The problem of polyphonic music transcription has received a huge amount of attention in the research community in recent years and I would go as far as to say that for single instrument polyphonies (piano, guitar etc...) the results are very good. Here are a few papers/authors who have looked at this problem deeply. Derry Fitzgerald has done a lot in the area, a lot of his NMF work on source separation produces accurate transcriptions. Anssi Klapuri has looked even more specifically at the problem of chord detection within multi instrument polyphonies. And finally, although not published, Mikel Gainza has developed very accurate chord transcription algorithms for commercial music soon to be released in a guitar based software product Riffstation.
The publications in the linkes here should give you a good idea of how the Polyphonic music transcription landscape currently lies.


Why do we need 4 points for homography but 7/8 points for fundamental matrix calculation? Sorry if this is a trivial question but I could not find a suitable answer elsewhere.

I recently learned about homography and fundamental matrices in a CV class. I understand how these matrices are calculated and how the the formulae are derived. I also understand that both homography and fundamental matrices have 8 degrees of freedom.

But I'm unable to understand why we need 4 point correspondences for finding the homography while 8 points for calculating the fundamental matrix.

Shouldn't the number of correspondences needed be same intuitively?

It is because in the case of fundamental matrix, each correspondence point relates to only one constraint(i.e it maps a point to a line in other image). Hence 8 correspondence points are required.

But in the case of homography, each correspondence solves two constraints. Hence only 4 correspondence points are sufficient.


Filtering PCM data in java I'm working on an audio based application for android platform , I'm trying to use a low pass filter to a fixed or variable max frequency "depends on the cost of development".
first I've been playing around with FFT and I believe I can get results out of it "by muting the frequencies above the required frequency and then do inverse FFT".
but the processing time seems high so I've searched again and found that a real time LPF can be implemented and applied to PCM data "IIR , Butterworth".
any suggestions on the best approach keeping in mind low processing speed "CPU" and small available memory "around 20MB for dalvik virtual machine"

Audio EQ in by FFT?!? Please don't do that. CPU usage will be the least of your worries. http://blog.bjornroche.com/2012/08/why-eq-is-done-in-time-domain.html There was also a related discussion recently on the linkedin dsp group.

This post covers how to do a bell filter in the time domain, but it's extremely similar to create a low pass (actually a low pass is easier), and will be MUCH easier than trying to do an FFT (which won't do what you want anyway): http://blog.bjornroche.com/2012/08/basic-audio-eqs.html

You won't have to worry at all about memory (it only needs 2 memory slots) and CPU usage is very low. Plus latency is minimal. The final implementation is like 10 lines of code.

The original "cookbook" reference is here: http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

Now, that said, this filter may not do for your purpose, but you haven't really made it clear what that is, so please feel free to clarify.


Does a filter add oscillations to a signal? this is a more general version of a question which i asked previously.

From what I understand, the purpose of a filter is to change the amplitude of a specific frequency band.

As I went through the analytic calculation of a filtered (2nd order) signal I observed that the output signal contains an oscillating component (a sinusoid, with the frequency being the imaginary part of the complex conjugate pole pair).

A complex pole (in the left half plane) does correspond to a damped oscillation. I really cannot understand how the addition of a specific frequency to the signal conforms with the job of a filter. Thanks for your insights!

EDIT:

The filter I am considering is a second order low pass Bessel Filter. The input signal is an exponential decay. I tried to calculate the output signal and observed there was a exponentially decaying cosine involved with its frequency being the imaginary part of the two poles. It really surprised me since I did not expect an oscillation in a filtered signal, which it did not express before.

I do understand, that this sinusoid was present in the signal already before filtering, it just struck me that it appears in the filtered signal as a single term.

A linear filter can increase the magnitude of frequency spectra already present in a signal.  Note that the Fourier Transform of any (reasonable) waveform or signal (including ones that do not look at all like they are oscillating) will decompose it into sinusoids (which do oscillate).  Amplify any one frequency, and an oscillation may start to stand out and make it look like it just appeared; but a lower amplitude form of that spectra was already hidden in the input signal to start with.

It wasn't added.  It was amplified.  Which is possible with an LTI filter.


Problem with Generalized Cross Correlation (PHAT filter) I have to calculate time delay between two signals.I am using the GCC-PHAT method. First, I cut the signals into smaller parts (~0.5s) and then try to calculate the delay. When I display those two signals on the graph, the delay is very obvious.

The problem is that for some parts of that long signal, my program calculates the correct delay, but for some parts the delay is totally incorrect. In both cases, when I display those two signals on the graph, the delay seems very obvious.

Relevant code:

duration=20000
a=hanning(duration+1)
[signal1,Fs] = audioread(mic1);
[signal2,Fs] = audioread(mic2);

start=1000;
stop=start+duration;
x1=[a.*signal1(start:stop)];
x2=[a.*signal2(start:stop)];

A=fft(x1);
B=fft(x2);

PHAT=A.*conj(B)./(abs((A).*conj(B)));
result=abs(ifft(PHAT));

Examples:
1) Correct result:

2)Incorrect result:

I recorded those signals in a closed room, so there are probably also reverberations.
Am I doing something wrong? Or is this method just not very effective in my circumstances?

I also tried to prefilter signals, but the results were the same.
There is also no improvement, if I change the length of the signals.

Using the gcc-phat method, you can use the matlab function "gccphat(s2,s1,fs)" (R2016a) with s1 as the reference signal, s2 the other one and fs the sampling frequency.

Basically it will do a similar correlation calculation :

  Ncorr = 2*N-1;  % N size of the signal
  NFFT = 2^nextpow2(Ncorr);
  R12 = bsxfun(@times,fft(s2,NFFT),conj(fft(s1,NFFT)));
  r12_temp = fftshift(ifft(exp(1i*angle(R12))),1);
  r12 = r12_temp(NFFT/2+1-(N-1)/2:NFFT/2+1+(N-1)/2,:);

The only difference with your code is that the matlab function use the nearest power of two sequence length for the fft (as Peter.K suggest with "Nfft = length(x1) + length(x2) - 1") and doing an exponentiel stuff (and this part make the difference with a regular correlation, what you were doing)

But to return the time delay, the function is doing an other trick :

  lags = (-(Ncorr-1)/2:(Ncorr-1)/2).';
  lags = lags/fs;
  [~,idx] = max(abs(r12));
  tau = ts/(2*fs)+lags(idx); % ts as the total sample of the signal

I'm doing similar experimentation and this method is working great for long and short non-moving signal. the reverberation of the room may cause some issue, but if you can see an obvious delay when displaying the signal it should work fine.

Right now, i'm having difficulty to return a good time delay with a moving source (the function retourn 0 or +/- the length of the signal) And i'm still working on it... So if anybody got a lead, it will be much apreciate ;-)

Hope this will help you.

Peace

PS : Also, when doing an fft/ifft you really should use fftshift/ifftshift because the matlab fft/ifft is mixing some stuff during the process.


Zero padding and 2D Fourier transforms: how does zero-padding affect phase? It's pretty clear that zero-padding an image before performing Fourier transform simply enlarges the magnitude image (stretching it to the new, padded size).

What I can't understand is how it affects the phase... At least to me it does not look like simple enlargement of the phase image. Could somebody clarify?

Symmetric zero-padding (in the center of an image around the N/2,N/2 sample) does not affect the FFT phase result.  Or after an 2d fftshift before the 2DFFT, symmetric zero-padding around the edges (circularly around the 0,0 sample) does not add phase shift.

An FFT phase measure the even to odd ratio around 0,0.  This won't change with any padding that maintains the identical even and odd decomposition of the input.  An fftshift moves the phase reference to the middle of a 2D image, which may make more sense with reference to an odd and even function decomposition of the original input.


Real-time single speaker recognition I'm developing an application for PC that records microphone data in real-time, and plays it back after applying some effects.

I'd like to apply these effects only when a single known speaker is talking. The environment might include other people's voice and indoor noise (for example, a classroom or a restaurant).

Since that single speaker is also holding the microphone, I tried to simply filter input data based on a (configurable) volume level. This method is working but not very well, since sometimes the environmental noise might go above the threshold, or the speaker's voice might go below it.

Is there a better way to recognize a single speaker? I've read about speaker identification using feature extraction and HMM/GMM, but I'm not sure whether it is an overkill for my case.

I apologize if this it too vague. I'm new to DSP so I'm not really sure what additional information is needed to answer this.
Thanks!

Edit:

Speakers usually speak one after another, but background chatter is possible.
In addition, although the speaker is known, I do need some kind of an initial setup procedure (to adjust the application to that speaker).

I was going to suggest the machine learning route but a quick search turned up OSS so you might not have to re-invent the wheel: Identify speakers with sndpeek

Do the speakers talk over each other (which complicates things) or is it mostly one after another?


Extended Kalman Filter in mechanics, electronics and hydraulics? Extended Kalman Filter is most used in GPS and navigation systems. But how much do I gain to switch to Extended Kalman Filter to the linjear Kalman Filter if I do LQGI controller for the industry?

All physics are nonlinear like hydraulics and mechanics.
What do you prefer?

It's the application which dictates whether you should use a linearized, extended, unscented Kalman filter or you can go easy with the simpler linear Kalman filter.

In other words, the extended Kalman filter is not preferred over the linear one because it provides better performance than the latter. Just the opposite; its estimation error is larger than what could have been provided by an equivalent linear Kalman filter if the problem could be solved by the linear.

When the physical nature of the problem does not permit a linear mathematical model, then you should either abondon the Kalman filter at all, or you should find a way around (such as the linearization of the nonliner model) to continue using the Kalman filter and still expect its benefits unless you have some alternative algorithm to go with.

Most fundamentally the accuracy of the extended Kalman filter depends on the the sampling interval (step size) used to update the linearized matrices and other parameters. Hence you should use dense enough sampling so that the nonlinearity does not change too much to add excess linear apprximation error. However smaller sampling interval also means more computations per second.

Also subtle issues such as noise de-coupling, independence assumptions, dimension reduction have thier own solution strategies which can also affect the accuracy of the filter estimates per given sampling rate.

When your problems are inherently nonlinear, then you cannot use the linear Kalman filter, therefore you can either use the extended Kalman filter or think about other algorithms that could be useful.


Face Detection + Cropping Project: Face Detection

Description: I want to detect and crop a face in an image. The image is captured through webcam and only one face per image.I used OpenCV face detector, but I was not satisfied with the cropping. So, I started using STASM (http://www.milbo.users.sonic.net/stasm/) face landmark detector to crop the image.STASM uses OpenCV face detector to find face in an image and STASM locates landmarks in faces. In bad light conditions, the cropped image from STASM is not good as it is not exactly detecting the face alone.

1) I want to know any better algorithm for face detection. My main aim is to crop the face from an image.

2) currently I am using STASM for cropping. In bad light conditions or when in an image, if the whole or complete face (forehead to chin) is not captured, STASM cropping is not reliable (The output will be only eye or lips). And in my application, if there is no proper output from the stasm or if the face is not cropped prpoerly then i should reject the images. How to do that? So I am planning to validate the face in an image by finding the Eyes. If I am right in my approach, how to detect the eyes from the cropped image?

FYI... STASM uses OpenCV's face detector (albeit very slightly modified for better image-boundary detections). If you really want a better face detector then I recommend hunting down detectors from the FDDB dataset challenge (most probably proprietary) or create your own custom detector to fit your needs.

If you still plan on going the route of verification via eyes, then OpenCV has a cascade for that, see http://docs.opencv.org/doc/tutorials/objdetect/cascade_classifier/cascade_classifier.html


Practical question about DFT in Matlab Please consider this piece of code:

clc;clear all;
fs=10^3;
t=linspace(0,1,fs);
f=2;
x=sin(2*pi*f*t);

W=dftmtx(fs);

X=fft(x);

X(3)
X(end-1)

%check : imaginary parts are 180 degrees dephased but right frequency
%plot(imag(W(end-1,:)));hold on; plot(imag(W(3,:)))

Xrec=(1/length(X))*(X(end-1)*W(3,:)+X(3)*W(end-1,:));

plot(x,'+g');hold on;plot(Xrec,'r','LineWidth',1.25);

The idea is to show (for my personal understanding) that we can reconstruct the pure sine wave x (of frequency 2 cycles) with the only 2 Fourier coefficients (corresponding to the 2 spikes we see in the spectrum) i.e. a negative frequency = -2 and a positive frequency = 2, for this example.

It seems to work well with the line:

Xrec=(1/length(X))*(X(end-1)*W(3,:)+X(3)*W(end-1,:));

Where I want to show that the signal x can be reconstructed by the weighted contribution of W(3,:) and W(end-1,:) (which I found to correspond to +freq and -freq (+2/-2)). This is the point of view of linear algebra where the vector x is a linear combination of the basis vectors W(3,:) and W(end-1,:) with the complex weights X(3) and X(end-1).

Is this explanation correct?

Moreover, there is the perspective of "matched filter", i.e. correlation. From this point of view we should have only contribution of the 2 frequencies aforementioned here 2/-2. But I still get some non-zero correlation when doing e.g. :

W(555,:)*x' =  -0.0063 - 0.0011i

... i.e. W(k != 3 or 999) (which are the indexes corresponding to freq=2 or freq=-2)

its small but not zero, so I would then to think that these have some contribution to the signal yet the should not ?

I augmented and "fixed" your code a little bit.

Your explanations are basically correct. Your main mistake is the following: your signal needs to be periodic in time domain, i.e., a whole number of periods must fit into your time window. If you use linspace this will not work: your first data point is 0, your last data point is 1, however cos(2*pi*f*0) = 1 = cos(2*pi*f*1). Hence your last and your first point are equal, they appear twice.

Try it with a very coarse spacing: cos(2*pi*linspace(0,1,5)) = [1, 0, -1, 0, 1]. Periodic copies give [1,0,-1,0,1,1,0,-1,0,1]. What you want is [1,0,-1,0,1,0,-1,0,1] though and for this you'd need to remove the last point.

A better way of doing it is this:

clear;
fs = 10^3;     % sampling frequency (Hz)
t0 = 1/fs;     % sample spacing in time (s)
NSamp = 1000;  % number of samples
f0 = fs/NSamp; % spacing in frequency
t_axis = (0:NSamp-1)*t0; % time axis
f_axis = (0:NSamp-1)*f0; % frequency axis
f_axis_symm = (-NSamp/2:NSamp/2-1)*f0; % frequency axis (symmetric)

f = 2;
x = sin(2*pi*f*t_axis);
X=fft(x);
C_two = X(2 + 1); % C[2]: need to add one since Matlab counts from one
C_minus_two = X(-2 + NSamp + 1); % C[-2]: need to add N (periodicity) and one (s.a.)

figure(1);
clf;
% stem(f_axis,abs(X)); % "asymmetric version [0,fs)
stem(f_axis_symm,fftshift(abs(X))); % symmetric" version [-fs/2,fs/2)
xlabel('f [Hz]');

W = conj(dftmtx(fs))/NSamp; % IDFT matrix: see help DFTMTX
Xrec = W(:,3)*C_two + W(:,NSamp-1)*C_minus_two; % reconstruct

figure(2);
clf;
plot(t_axis,x,'s','MarkerIndices',1:10:NSamp);
hold on;
plot(t_axis,Xrec,'b--');
xlabel('t [s]');

The code also shows the spectrum, plotted in a symmetric way so that you can see exactly two lines popping up at -2 Hz and 2 Hz. It also shows how to reconstruct the signal. Note that you need an IDFT matrix for reconstruction. See help dftmtx for a comment regarding this.


Radio Signal Stationarity A radio signal recording of a wireless communication system (e.g: Wi-Fi traffic) is beaconized, channelized and subject to noise.

When working with such an RF signal, numerically transformed to a stochastic timeseries, to which extent can I consider the signal to be stationary?
What properties are relevant from a communications point of view?

When working with such an RF signal, numerically transformed to a stochastic timeseries, to which extent can I consider the signal to be stationary?

That depends on your signal model. We can't tell you that – but for example, in time-slotted system, obviously the received signal can't be stationary – its variance (which is basically its power) depends on the time you're looking at.

We tend to model noise as (weak-sense) stationary; same for data coming out of the source coding; we often even technically enforce that using whitening scramblers.

However, that really depends: some systems use whitening. Others don't.

Some techniques assume all disturbances to be stationary. Others are built around the specific notion that interference is not stationary.

So, model your signal.

  Any pointer papers about signal stationarity in telecommunications are appreciated.

Not an issue for papers: this is subject for standard textbook discussion; Goldsmith and Viswanath/Tse would be my go-to books here.

(Generally, if you feel like something is so important that it is a basis for all that you do, it's almost certain people have looked at it a long time ago and have already written a good book on it. The trick is finding the good books; in this case, things are so fundamental that basically all wireless communication books concern themselves with it.)

  What properties are relevant from a communications point of view?

All and none. Defined in the signal model.


Recommendation for courses / studies on digital signal processing I hold a master's degree in mechanical engineering. However at my job I am more and more diving into topics of signal processing and data science. I find it great to discover about new topics and to learn on the go as I work on projects. During my study I had some basic courses on which I can build on. But sometimes I wonder if it wouldn't be easier and also more fun to for instance do a 2nd master or something similar on the topics I am trying to learn. I think my employer would support me if I want to and as I liked studying in the past I think that could be a great experience.
So with that the question arises what type of course/study would be a good fit? Ideally it would cover relevant topics* for my work/interest, be able to partake in part-time, be online (or taught block-wise in central europe) and include collaborative elements (I don't want to just learn and work by myself without interactions).
*some topics I recently dealt with that I would like to know more about:

Singular Value Decomposition
Kalman Filters (especially in the context of inertia measurement units and nonlinear problems)
Principal Component Analysis
Digital Filter Design
Control theory (implementing a control algorithm and tuning it for a problem)
Running a control algortihm on a microprocessor
Machine Learning for Audio Processing
Wavelet Transforms
...

I've kind of grouped your subjects into larger overall subjects.
Note that there's a lot of overlap here, with the possible exception of actually making it work in a microprocessor (except -- in my opinion the best person to implement something is someone who understands it.  So -- overlap).
Specifically, you could claim that it's all applied math.  Or all signal processing.  You could even stretch this list a bit and wedge it into control theory, but it wouldn't be an easy fit.

Applied Math

Singular Value Decomposition
Principal Component Analysis
Machine Learning for Audio Processing

Signal Processing

Digital Filter Design
Kalman Filters (especially in the context of inertia measurement units and nonlinear problems)
Wavelet Transforms

Control Theory

Kalman Filters (especially in the context of inertia measurement units and nonlinear problems)
Control theory (implementing a control algorithm and tuning it for a problem)

Embedded Software (or systems) Engineering

Running a control algorithm on a microprocessor

A weakness of the US tertiary educational system is that all of these subjects will be used in engineering departments, but they're all really applied math.  So different universities will group them differently, and it can be hard to find a program that will cover them all.  I'm not sure if that applies to where you live.  I happened to end up -- by luck more than design -- at an institution that let me get enough of a grounding in the related subjects that I could learn the rest on my own (Worcester Polytechnic Institute -- it's a great place).
I would suggest that you approach this by finding a university program that includes the classes you need to get a grounding in this, and get your second master's.  Don't expect to learn all of this within a Master's program -- I can lay claim to about 90% of this list, but about half of it are things that I learned on the job, over the last three decades.
Classes to look for are below.  Note that you're not going to have time for them all, so you'll need to pick and choose.  If you want to grow into your full list, take as much math as you can and expect to do self-study:

Stochastic signals and systems (basically, this is the study of finding the optimal system to process a signal with random components).

If you haven't taken it, a 4th-year course on statistics from the math department.

Real analysis (AKA "Advanced Calculus", but I think that's an obsolete name).  It's an invaluable help when you're sitting at your desk at work wondering if the nifty new thing you're applying is actually mathematically sound.
If it's offered, estimation and detection theory.  This is basically a follow-on to stochastic signals and systems.  It's a deep dive into the math underlying the detection of signals, signals in noise, etc.  The one that I took had you derive the basic Kalman filter from first principles in a homework problem -- which gives you an idea of the depth of the subject.
Any class on optimal state estimation (i.e., Kalman filtering and all its variations).
Signal processing, assuming you haven't taken a class in it.
Digital signal processing, ditto.
Control theory classes will just be called "control theory", but you may not be able to fit much in.  If there's a class offered in state-space systems and you've already got classical (transfer function) control under your belt -- take it.
If you want to specialize in control theory, go for nonlinear control or whatever is offered.  But if you do that, you probably won't have time for the Detection & Estimation class, so think hard.

I'm not sure what to suggest on the implementation front -- that's basically yet another specialty.  I just picked that up on my own, on the way.  But, I started out as a kit with a hobby in electronics and a microprocessor board, so I was already doing basic embedded programming when I was 13.  Then when I got an EE degree, any class that involved digital circuits was either "oh, yes, I've done this", or "oooh!  what a neat formal way to do what I've struggled through intuitively".


Real Time Depth Estimation in MATLAB I want to estimate depth from live video stream with two calibrated webcams in MATLAB but has no results. I have Webcam Support Package and Image Acquisition Toolbox but I don't know how to use live video stream via stereo calibrated webcams against recorded video file. In this post depth is estimated with two recorded video files.

How can I estimate depth in live video stream?

You would have to connect both cameras to your computer, and either create two webcam objects or two videoinput objects to read images from them.


"Undersampled" Polyphase Filter Bank I'm working on an application where I want to grab high-resolution pieces of the spectrum, but I don't need the whole spectrum-- just a few bands. I've looked into using the Goertzel algorithm (too slow; with the number of points we want the regular FFT is faster), pruned FFTs (doesn't look like much of an improvement in performance), and polyphase filter banks (PFBs) (the problem being that it requires an FFT of the same size that I would have otherwise already computed, defeating the point).

The books that I've looked at so far ("Digital Signal Processing, 4th ed" by Proakis/Manolakis and "Signal Analysis: Wavelets, Filter Banks, Time-Frequency Transforms and Applications" by Mertins) both focus on either "critically sampled", where the decimation rate $D$ equals the number of subbands $N$ or "oversampled" where $N > D$.

Is there a way to produce a computationally efficient PFB such that $D > N$? In my application the loss of information is unimportant because I only need a few subbands of small bandwidth.

To analyze a slice of the spectrum:

Use an IQ mixer to translate the centre of the slice, $F_c$, to 0 Hz. This is done by multiplying the signal with $e^{-j2\pi F_c/F_s}$, where $F_s$ is the sample frequency.
Use a decimation low-pass FIR filter to remove the signal components outside the slice. The pass band of the filter is half the slice width. Output from this process is the down-converted slice sampled at a low sampling rate (just meeting the Nyquist criteria).
Apply your favourite window function and FFT to the time signal.

This process is repeated for each of the wanted slices.


Finding an Image's Orientation With DFT Frequency Amplitudes I'm looking to get some general idea of an image's orientation so that I can then rotate it to the nearest 90 degree angle.

My idea on how to do this is to take the DFT on the image, then do a "radar sweep" to find which angle has the greatest sum of frequency amplitudes along the line, and call the winning angle the orientation.

I scan from 0 to 180 degrees instead of from 0 to 360, and make the radar sweep double sided since negative and positive frequencies seem to belong together.

For example, if i take the DFT of this image:

The frequency amplitudes are:

And if I do a radar sweep, I'd find the winning line / angle here:

Using that information, I can know where the nearest 90 degrees is for that blue line, figure out how much i'd need to rotate the dft image to make that blue line vertical, but instead rotate the source image by that amount, to make that image orient towards the nearest 90 degrees.

Something like this, but actually level (I'm working out some bugs at the moment, in my implementation)

Is this method sound?  Are there any complications that I'm going to encounter when moving to more complex images such as complex patterns or photographs?

Edit: In the above I didn't explain the reasoning for what I am doing.

If you take the dft of simple images - like a single line, or a repeating pattern of black and white stripes - the resulting frequency amplitudes image has data that runs perpendicular to the direction of that line or those stripes.

Line this:

DFT'ing to this:

May (possibly naive!) understanding is that I should then be able to tell which direction (from the center) has the highest sum, and that will tell me the angle perpendicular to the largest represented orientation of the image.

From there, i can find how much i'd have to rotate that line to get to the nearest 90 degrees, and apply that rotation to the source image, to make it line up to it's nearest 90 degrees.

Is my thinking too simplistic, or is this a solid approach?

Have you heard of the Hough Transform? It is a possible alternative for tackling your problem.