Derivative of delta function I've come across the below. Can you please help me understand the proof for this?

$$\int_{-\infty}^{\infty}x(t)\delta'(t-2.5)dt=-\frac{dx(t)}{dt}{\Big |}_{t=2.5}$$

You can see this by using integration by parts:

$$\begin{align}\int_{-\infty}^{\infty}x(t)\delta'(t-T)dt&=x(t)\delta(t-T){\Big|}_{-\infty}^{\infty}-\int_{-\infty}^{\infty}x'(t)\delta(t-T)dt\\&=x(T)\delta(t-T){\Big|}_{-\infty}^{\infty}-x'(T)\\&=-x'(T)\end{align}$$

where it is assumed that $x(t)$ and $x'(t)$ are continuous at $t=T$.


Estimating number for iterations for gaussian smoothing I have some data sets on which I applied Gaussian smoothing using [1 4 6 4 1] kernel. In my program I iterated this kernel 50 times on the data sets. But only a few of them could retain their shape, most of them just flattened and lost all the patterns

I want to know is there any way I can estimate the number of iterations of smoothing for each data set ?

Should I be using Normalized Squared Error for the data set and its smoothed version , or maybe use any type of Correlation between the two ?

Short Answer

Just compute a gaussian blur kernel for a given sigma, then perform single iteration of 2D convolution or two 1D convolutions (for rows, then for columns or vice versa).

Don't forget to normalize the kernel!

If the kernel comes out too large (it is typically cropped to $\lceil 3\sigma\rceil$), compute the convolution in frequency domain (i.e. IFT(FT(image) *. FT(kernel)) where *. is a point-wise multiplication).

If you don't care about precision, you can apply several iterations of a box filter (the more iterations, the better approximation to Gaussian blur, but the more blurring occurs). You can find a nice definition in SVG Specs. Box filter is separable (can be done as two 1D filterings) and can be implemented with running sum (filter size is not a problem) - this is probably the fastest approximate smoothing method possible.

Long Answer

SInce you have the kernel in integers, it is only approximation to gaussian. To imagine what happens after applying by source_kernel 50 times is equivalent to applying this target_kernel once:

target_kernel = IFT( FT(source_kernel) ^. 50 )

where ^.50 is point-wise power function.

To learn more about the iterative approach, note that Gaussian blurring is equivalent to applying heat equation on an image (imagine how light "hot" pixels diffuse into darker, cooler pixels).

The heat equation is partial differential equation (PDE) in case of images.

You can apply certain amount of smoothing by specifying time amount for simulating heat dissipation and then applying the heat equation in small time steps. I don't know how time amount relates to gaussian $\sigma$, but you probably find more info in any introductory material to diffusion image filtering using PDEs.

Each time step requires computing finite differences at every pixel and applying a small stencil - it is somewhat similar to convolution.

This iterative approach is called implicit filtering and can give you very precise amount of smoothing (the precision is limited by the size of your time steps).

Explicit formulation also exists but require solving very large system of equations.

This kind of filering is rather used for more comlplicated problems like anisotropic diffusion, inverse problems, total variation denoising etc.


deciphering the frequency response from the input data Newbie to Signal processing, sorry if the question is basic.
I have a data set consisting of 601 values. The values being time Vs voltage . The voltage values have transients in them. My aim is to filter out these transients using an FIR filter. But i'm still far from doing that.
This is my voltage Vs Time plot.

I plotted the fft of voltage and this is what i got.

I read up that fft basically divides my sampling frequency,Fs, into N equal parts (N= 601, Total number of values). According to my (voltage Vs Time) plot, my Fs=2500Hz (by calculating Ts as the difference between consecutive values of time from data set).Then,2500/601 gives 4.159 as the first point in my fft. Is this even correct?

Going by the above,the strength of the signal is greatest at the 11th value (ie,10*4.159 =41.59 Hz) and the 592nd value (591*4.159=2458 Hz).Now, What significance does this information hold if i need to design a filter to eliminate the transients? What do information do i take away exactly, by performing the fft of the voltage signal?

Thanks

Since the signal is most certainly real-valued, you can ignore the FFT values after $f_s/2$ (or sample 300 in your plot) since they will mirror the first ones. With a sampling frequency of 2500Hz, the highest frequency in your signal can be 1250Hz. If you have higher frequencies and don't filter them prior to sampling, you will get aliasing. Look up the sampling theorem for more information, it's quite important when working with sampled signals.

For your actual question I'm not sure what you want to achieve. With an FIR filter you can attenuate certain frequencies. If your transients are strictly not within your desired frequency range, that might work. If they are in the same frequencies as your actual signal, an FIR filter won't help. Is that supposed to be a transient in the first plot where the signal goes from being clean to jiggly?

If so, you could probably just use a lowpass filter which cuts off at 50Hz or similar.


Understanding of wireless communication system from MATLAB point of view I am dealing with a communication model described as

$$\overrightarrow{y(n)}=\pmb{H(n)}\overrightarrow{x(n)}+\overrightarrow{w(n)}$$

where $\pmb{H(n)}$ is the $M\times M$ fading coefficient matrix at time $n$. What is difference between time instant  $n$ and length $M$?

Below is a sample MATLAB/pseudocode to generate a 4x4 MIMO system, for 1000 Monte Carlo trials:

M=4; MC=1000;
H= randn(M,M,MC);
x=randn(M,1,MC);
w=randn(M,1,MC);
Y=zeros(M,1,MC);
for j = 1:MC
   Y(:,:,j)= H(:,:,j)*x(:,:,j) + w(:,:,j);
end

This would generate y according to the equation in the question. You can view the MC values as different time instances. Meaning different channel, transmit vector and noise matrix for each instance.


STFT: why overlapping the window? For STFT, we impose window of certain size onto the original signal, then we perform fft on each window. The uncertanty about frequency and time is determined by the width of the window, however, I can't understand what is the point of having overlap windows...

If we have a signal, for instance, why can't we just divide the signal into 6 trunks (non-overlapping window), and then we perform fft on each of those trunks?

Maybe let me make it more clearly in my application. I am going to mostly dealing with 60Hz power line wave, and occationally, we want to monitor the 180Hz transient effect at the power line. Since the signal will be mostly periodic, should I use window then?

We always want to apply some kind of a window function in order to minimize the effect of leakage. This makes rectangular window (lack of any windowing) case never used, this is why:
Any tapering function used is almost always decreasing to zero at boundaries.

This is why we are losing some data. In order to retrieve that somehow you will usually do 50% of overlap when processing. This will retrieve whatever was in between.

Another thing is that if you apply the Inverse STFT, you should use complementary window, that is summing to 1, i.e. Hanning with 50%.

Finalising - yes, you should pretty much always use windowing in your applications.

For more comprehensive informations please refer to great white-paper:

  Heinzel G. - Spectrum and spectral density estimation by the DFT,
  including a comprehensive list of window functions and some new
  flat-top windows


Linear response function for a system with derivative: $U=L \frac{d I}{dt}$, expressing $U=f(I)$ I have a super basic questions. I am a not really into signal processing (more about physics), but I would like to understand an aspect of linear response function (I think the question fits for this forum).

From my basic understanding, any linear, time invariant, causal system can relate the output $S(t)$ to the input $E(t)$ under the following relations:

$$\begin{equation}S(t)=S(-\infty) + \int_{-\infty}^{+\infty} \chi(t-t') E(t') dt'\label{eq_1} \end{equation}$$

The causality holds when imposing $\chi(u<0)=0$.

Let's consider a very simple case: voltage around an inductance. I have the law: $U=L \frac{d I}{d t}$

If I express $I$ as a function of $U$, I can write down:

$$I(t)=I(-\infty) + \int_{-\infty}^{+\infty} \frac{U(t')}{L} dt'$$

My response function will simply be $\chi(u)=\frac{\Theta(u)}{L}$ (where the heaviside is for causality).

But the system in which $I$ is the input and $U$ the output is also linear. Thus I would expect to be able to express:

$$U(t)=U(-\infty) + \int_{-\infty}^{+\infty} f(t-t') I(t')$$

However because of the derivative on $I$ in the law, I don't find how it is possible.

Am I wrong in my initial statement and \ref{eq_1} ?

Let's start with the simple system

$$y(t)=x(t)\tag{1}$$

where $y(t)$ denotes the output and $x(t)$ denotes the input. It is linear and time-invariant, and, consequently, it must be possible to represent it by the convolution integral

$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau\tag{2}$$

where $h(t)$ is the system's impulse response. In this case, it turns out that $h(t)=\delta(t)$, i.e., the impulse response is a Dirac delta impulse.

Now, for a differentiator

$$y(t)=x'(t)\tag{3}$$

we "simply" (conceptually, at least) have $h(t)=\delta'(t)$, where $\delta'(t)$ is the (generalized) derivative of $\delta(t)$. So, using the concept of a generalized derivative of a distribution, we can write the input-output relationship of an ideal differentiator as the following convolution integral:

$$y(t)=\int_{-\infty}^{\infty}x(\tau)\delta'(t-\tau)d\tau\tag{4}$$


(For c-FastICA) On covariance and pseudocovariance matrix of a complex random vector I am currently studying complex FastICA and the paper says that

  Suppose $\mathbf{s}$ is a $n\times1$ complex random vector. If $\mathbf{s}$ has zero mean, unit variance, and uncorrelated real and imaginary part of equal variances, then $E[\mathbf{s}\mathbf{s}^H]=\mathbf{I}_n$ and $E[\mathbf{s}\mathbf{s}^T]=\mathbf{0}_n$.

I don't quite get how $E[\mathbf{s}\mathbf{s}^H]=\mathbf{I}_n$ and $E[\mathbf{s}\mathbf{s}^T]=\mathbf{0}_n$ come about from the conditions.

We have the covariance matrix as
\begin{align}
\operatorname{cov}(\mathbf{s})
&= E[\mathbf{s}\mathbf{s}^H]-E[\mathbf{s}]E[\mathbf{s}^H] \\
&= E[\mathbf{s}\mathbf{s}^H]-\mathbf{0}_{n\times1}\mathbf{0}_{1\times n}\\
&= E[\mathbf{s}\mathbf{s}^H]\\
\end{align}
and the pseudocovariance
\begin{align}
\operatorname{pcov}(\mathbf{s})
&= E[\mathbf{s}\mathbf{s}^T]-E[\mathbf{s}]E[\mathbf{s}^T] \\
&= E[\mathbf{s}\mathbf{s}^T]-\mathbf{0}_{n\times1}\mathbf{0}_{1\times n}\\
&= E[\mathbf{s}\mathbf{s}^T]\\
\end{align}
I don't quite get how to equate the last line of covariance matrix to identity and the pseucovariance to zero.

If I were to write out the matrix,
\begin{align}
E[\mathbf{s}\mathbf{s}^H]
&=E\left\{\begin{bmatrix}
s_1s_1^* & s_1s_2^* &\cdots & s_1s_n^*\\
s_2s_1^* & s_2s_2^* &\cdots & s_2s_n^*\\
\vdots & \vdots &\ddots & \vdots\\
s_ns_1^* & s_ns_2^* &\cdots & s_ns_n^*\\
\end{bmatrix}\right\}
\end{align}
and
\begin{align}
E[\mathbf{s}\mathbf{s}^T]
&=E\left\{\begin{bmatrix}
s_1s_1 & s_1s_2 &\cdots & s_1s_n\\
s_2s_1 & s_2s_2 &\cdots & s_2s_n\\
\vdots & \vdots &\ddots & \vdots\\
s_ns_1 & s_ns_2 &\cdots & s_ns_n\\
\end{bmatrix}\right\}
\end{align}
I still can't quite figure how all of these eventually becomes identity and zeros.

You have

$$E[s_ks_k^*]=E[|s_k|^2]=1$$

because the complex random variables $s_k$ have zero mean and unit variance. That means that all elements of the main diagonal of $E[\mathbf{s}\mathbf{s}^H]$ equal $1$. Furthermore, with $s_k=x_k+jy_k$we have

$$\begin{align}E[s_ks^*_l]&=E[x_kx_l+y_ky_l+j(x_ly_k-x_ky_l)]\\&=E[x_kx_l]+E[y_ky_l]+jE[x_ly_k]-jE[x_ky_l]\\&=0,\qquad k\neq l\end{align}$$

because real and imaginary parts are uncorrelated. Consequently, all off-diagonal elements of $E[\mathbf{s}\mathbf{s}^H]$ are zero.

The off-diagonal elements of $E[\mathbf{s}\mathbf{s}^T]$ are zero for the same reason. (There's just a sign difference in the sum of the expectations, but since each of them is zero the result is the same). The main diagonal elements are

$$\begin{align}E[s_k^2]&=E[x_k^2-y_k^2+2jx_ky_k]\\&=E[x_k^2]-E[y_k^2]+2jE[x_ky_k]\\&=0\end{align}$$

because real and imaginary parts are uncorrelated, and they have equal variance, i.e., $E[x_k^2]=E[y_k^2]$.


How to get the signal to noise ratio from noisy data set with sinusoid? I'm new to DSP and trying to do some example problems for a class and am confused by this one.

Given a data set, and knowing it contains a signal that is a single sinusoid with random content noise, how can I calculate the signal to noise ratio? I also have the sampling rate.

I have been able to plot the autocorrelation, and can see the characteristic peak at 0 to indicate the noise. I'm just unsure if I can get the power of the signal and noise from this autocorrelation, or the power spectra density (taking the fft of the autocorrelation). Do I need to find the amplitude, frequency, and phase of the sinusoidal signal and subtract it from the original to do this? If so, how do I find those?

Any help is appreciated, I am new to DSP and learning to thanks!

One method is to pass the noisy sinewave through a notch filter tuned to the sine wave frequency. The output of the notch filter is N and the input is S + N. S/N = (S + N)/N - 1. This is the way some analog meters work.


Why does the frequency window affect the inverse fourier transform oscillation frequency? I have an oscillating pulse in the frequency domain that I would like to inverse Fourier transform. My signal looks like:

Which was coded in MATLAB using the following code:

sampleRate = 1000;
freq = 500 : (1/sampleRate) : 1500;
intensityFreq = exp(-((freq-1000).^2));
signalFreq = sqrt(intensityFreq).*exp(-1i*10*(freq-1000));
plot(freq,signalFreq)

When I inverse Fourier transform, it shows a similar oscillating pulse (in time) whose oscillations depends on the frequency range I used in the above code. For example, when:

freq = 500 : (1/sampleRate) : 1500;

is used (as in the code shown above), the inverse fourier transform looks like:

But when:

freq = 980 : (1/sampleRate) : 1020;

Is used instead, the inverse Fourier transform looks like:

This makes no sense to me, since the frequency range shouldn't matter since signalFreq goes to 0 around these points anyway. Why does the altered frequency range affect the oscillations in the inverse Fourier transform so much? How do I obtain the "actual" time domain signal that doesn't depend on frequency windows?

Any help is much appreciated. If this post is not clear I would be happy to provide more detail.

EDIT: The full code used to find the inverse Fourier transform is:

sampleRate = 1000;
freq = 500 : (1/sampleRate) : 1500;
intensityFreq = exp(-((freq-1000).^2));
signalFreq = sqrt(intensityFreq).*exp(-1i*10*(freq-1000));
Y = ifft(signalFreq);
plot(real(Y))

ifft works on the full spectrum as pointed by Hilmar.
In the discrete description of Fourier Transform theory, it means that each sample supplied as the argument to ifft is a frequency bin in the range $[0,1]$ in cycles per sample (see normalized frequency if you do not know what that means), which corresponds to $[0,2\pi]$ radians per sample or $[0, F_s]$ sample per second.

Basically what ifft saw is a spectrum ranging from (assuming $F_s=1000 Hz$, which as you see below is completely arbitrary):

0 to 1000 Hz in the first case, with a Gaussian in magnitude spectrum around 500 Hz
0 to 40 Hz in the second case, with a Gaussian in magnitude spectrum around 20 Hz

Although it gets a bit more tricky, as your Gaussian is peaking right at the Nyquist frequency, you are simply retrieving an oscillation in time-domain at the limit of your sampling rate. In both case your spectrum is the same shape but only of different length, if you assume the signal is of fixed duration that means you are transforming two signals with different sampling frequency.

It would be more appropriate to define first the frequency range as

sampleRate = 1000;
duration = 1; % 1 second
N = duration * sampleRate;
delta_f = N/sampleRate; % frequency resolution: this is fixed by all other parameters
freqs = 0:delta_f:sampleRate

% Then create your gaussian spectrum around whatever value (but below Nyquist!)
mu = 250;
sigma = 50;
intensityFreq = exp(-(freqs-mu).^2/sigma);
% etc... Play around with mu, center frequency of your gaussian, and sigma.

All in all, I would advise you to read a bit more about discrete Fourier transform and inverse Fourier transform.

Also, from your example code you will get a complex valued time-domain signal. If you were aiming for a real-valued signal, of course the magnitude spectrum supplied must be symmetric around the Nyquist frequency $Fs/2$ (and more precisely the complex valued spectrum must be complex conjugate). But this is another matter, and you should read more about Fourier transform theory if this is unknown to you (this answer is a good read on the matter).

% To transform the complex spectrum such that the second half is the complex conjugate of the first (avoiding DC 0th bin) do:
signalFreq = sqrt(intensityFreq)*exp(1i*freqs); % assuming intensityFreq is defined on one half only, 0 on range [fs/2:fs]
signalFreq(end-1:-1:end/2+1) = conj(signalFreq(2:end/2));

% You can use the option 'symmetric' to assume that the argument is conjugate symmetric:
signalTime = ifft(signalFreq, 'symmetric');
```


How to turn off calculating pyramids during image processing? I am using Erdas 2014 to run a 5x5 focal median filter on a batch of images (Raster Tab > Spatial > Focal Analysis (batch)). These images will eventually be resampled to a different pixel size, and therefore I do not want to spend the processing time calculating pyramids. How can I turn off pyramids while performing batch image processing in Erdas 2014?

Pyramids can be toggled on and off using the Erdas preferences: File > Preferences > File Formats > Image Files (General).


What mathematical function could be used to simulate this sound effect? I have this sound effect which by the way if anyone could give me its name, I would really appreciate it. And after loading it into matlab/octave (don't worry this isn't a programming question) using the following commands:
[x, fs] = audioread("Path To File");
plot(length(x(:,1)),x)
sound(x,fs) % To hear the sound

I get the following:

Could anyone tell me what mathematical function I could use to simulate this?

I guess you want to synthesize that sound, i.e. create a synthetic signal to be as close as possible to the original sound. Potentially, you can do this by creating a superposition of sine waves with varying amplitudes and phases. So you have to know the parameters for those sine functions.
From what I can see and hear, the sound has two signals, which are very similar (both visually and to ear). Moreover, if you look at the frequency spectrum, you will see, that they are quite similar as well. However, the spectrum seems quite populated with different frequencies, which will make it difficult to reproduce as a synthetic function. At least if using only the provided original sound file.
Here is the spectrum plot:

And here is a python code to make those observations. (The file with sound is downloaded to the folder of the python script.)
import matplotlib.pyplot as plt
import soundfile as sf
import sounddevice as sd

# get file
filename = "SoundEffect.wav"
x, f = sf.read(filename)
x1 = x[:,0]     # the first signal
x2 = x[:,1]     # the second signal

# play and show the original sound in time domain
sd.play(x, f)
status = sd.wait()
plt.plot(x)
plt.show()

# play the first signal
sd.play(x1, f)
status = sd.wait()
# play the second signal
sd.play(x2, f)
status = sd.wait()

# show the signals in frequency domain
plt.magnitude_spectrum(x1, Fs = f)
plt.magnitude_spectrum(x2, Fs = f)
plt.show()

Bear in mind, I am not familiar with how sound effects are done. Most probably, this sound effect is a more simple superposition of certain signal shapes. However, unless they are already known, it is hard to "reverse-engineer" them from the sampled version of the signal.


AR Modeling: Why residual is white noise? I was going through AR modeling.

The AR model of  a covariance stationary process can be expressed as:

$$x[n]=\sum\limits_{i=1}^{p} \alpha_i x[n-i] + \epsilon[n]$$

where $p$ is the model order and $\epsilon[n]$ is the residual

1) Why is it that the residual $\epsilon[n]$ is a white noise? Is it by definition or can it be shown to be a white noise?

2) Where is the assumption on the stationarity of $x[n]$ used? In other words, if $x[n]$ is not stationary, what would have changed?

Thanks,

it's a hand-wavy argument, but the idea is to derive the coefficients $\alpha_i$ so that the norm of the error $\Big\||\epsilon[n]|^2\Big\|$ is minimized.  we assume that $x[n]$ is a stationary "random" process, we know all of the previous samples $x[n-i]$ and we want to make a good guess at $x[n]$ with a linear combination of the previous samples.  we can only make some guess of the next $x[n]$ if it is not a white random process.  if it is white, knowing the previous samples is not going to help us at all.  if it is not white, then the autocorrelation of $x[n]$ is not zero for non-zero lags and from that we can compute $\alpha_i$.

now if $\epsilon[n]$ was not white, we could do further LPC to make a good guess at $\epsilon[n]$ and make a better guess for $x[n]$, but if the coefficients $\epsilon[n]$ are optimal, then we can't make a better guess for $x[n]$.  so then the error of our guess for $x[n]$ and what we really get has to be white.


Problems calculating Z-transform I am trying to solve a class exercise in which I am given the following, in Laplace domain:

$$G(s)=\dfrac{e^{-Ts}}{s+3}$$
$$H(s)=\dfrac{1}{s}$$

And I need to calculate $\dfrac{C(z)}{R(z)}$, which is defined as:
$$\dfrac{C(z)}{R(z)}=\dfrac{G(z)}{1+H(z)G(z)}$$

In order to do so, I have to transform $G(s)$ and $H(s)$ into the $Z$ domain and operate to compute $\dfrac{C(z)}{R(z)}$.

I have tried to compute $G(z)$ and $H(z)$, and I think I have failed in calculating $G(z)$, because I am not getting the right result.

I obtained that $H(z)=\dfrac{1}{1-z^{-1}}$ (that I think is correct).

Therefore, my question is how to obtain $G(z)$, I know that $\boxed{e^{sT}=z}$, but I don't know how to relate it with it, because initially I obtained: $G(z)=\dfrac{z^{-1}}{1-e^{-3}z^{-1}}$ but apparently is not correct. What am I doing wrong?

Note: in the statement of the exercise I am told to consider the period of one second (T=1s).

Solution: $$\dfrac{C(z)}{R(z)}=\dfrac{(1-z^{-1})z^{-1}}{1-(1+e^{-3})z^{-1}+(\dfrac{2}{3}e^{-3}+\dfrac{1}{3})z^{-2}}$$

We have
$$G(z)=\dfrac{z^{-1}}{1-e^{-3}z^{-1}}$$
and
$$H(z)=\dfrac{1}{1-z^{-1}}$$
and we want to calculate
$$\dfrac{C(z)}{R(z)}=\dfrac{G(z)}{1+H(z)G(z)}$$

So let's do that:
$$
\require{cancel}
\dfrac{\dfrac{z^{-1}}{1-e^{-3}z^{-1}}}{1+\dfrac{1}{1-z^{-1}}\dfrac{z^{-1}}{1-e^{-3}z^{-1}}} = \frac{z^{-1}(1-z^{-1})}{(1-z^{-1})(1-e^{-3}z^{-1}) + z^{-1}}\\
= \frac{z^{-1}(1-z^{-1})}{1 - z^{-1} - e^{-3}z^{-1} + e^{-3}z^{-2} + z^{-1}}
= \frac{z^{-1}(1-z^{-1})}{1 \cancel{- z^{-1}} - e^{-3}z^{-1} + e^{-3}z^{-2} \cancel{+ z^{-1}}}\\
= \frac{z^{-1}(1-z^{-1})}{1 - e^{-3}z^{-1} + e^{-3}z^{-2} }
$$
so it looks like either the solution you have is wrong, or that $G(z)$ is not right.

So, let's see what could be wrong. Either
$$
z \not = e^{sT}\tag{1}
$$
or
$$
s \not = 1 - z^{-1}\tag{2}
$$
or
$$
G(z) \not = \dfrac{1}{1-e^{-3}z^{-1}} \cdot {\cal T}_{s \rightarrow z}(e^{-s})\tag{3}
$$
where ${\cal T}_{s \rightarrow z}$ is the $s$ to $z$ conversion.

Comparing (1) and (2), it looks like (1) is assuming the starred transform and (2) is assuming the impulse invariance method --- which can't both be right.


z-Transform Methods: Definition vs. Integration Rule The definition of the z-transform is defined as $z = e^{sT}$ where "s" is complex frequency for continuous-time systems and "T" is the sample period.  Why are rules such as the forward rectangular rule, or Tustin's method used instead of the definition?

Forward rectangular rule:  $s \leftarrow \frac{z-1}{T}$
Tustin's rule:  $s \leftarrow \frac{2}{T} \frac{z - 1}{z + 1}$

I am asking about the characteristics of the transform methods.  If I transform G(s) to G(z), then why use Tustin's rule or forward rectangular rule instead of the definition?  I would think the integration rules are just an approximation to the definition.  What are the characteristics of these maps and how do they compare with each other?

Just to avoid a misunderstanding: the $\mathcal{Z}$-transform is a transform defined for sequences, comparable to the Laplace transform for continuous functions. What you are talking about is not the $\mathcal{Z}$-transform, but methods for converting analog to digital (actually, discrete-time) systems. [And it doesn't help that one of those conversion methods is called the matched Z-transform, which is just an unfortunately chosen name.]

There are several methods to achieve such a conversion, and all of them have certain properties. There is no ideal or "correct" mapping, and some methods may be more well-suited for certain applications than others.

The methods mentioned in your question are direct mapping methods, where the variable $s$ in the continuous-time transfer function is simply replaced by a (rational) function of $z$. Apart from the two methods mentioned in your question (1. Tustin's method - also called bilinear transform, and 2. forward Euler), there is also the backward Euler method, which replaces $s$ by $(1-z^{-1})/T$. Note that Tustin's method (bilinear transform) is equivalent to the trapezoidal integration rule. The important advantage of the bilinear transform for the conversion of analog systems is that the frequency axis ($j\omega$-axis) of the $s$-plane is transformed to the unit circle of the $z$-plane, and, consequently, the frequency response of the resulting discrete-time system is just a compressed (warped) version of the analog frequency response. One consequence of this is that analog frequency selective filters with optimal magnitudes (Chebysev, Butterworth, elliptic) are mapped to discrete-time filters satisfying the same optimality criteria. Furthermore, since the left half-plane is mapped to the inside of the unit circle in the $z$-plane, stability (and the minimum-phase property) are preserved by the transformation.

This is different for the forward and backward Euler methods. With those methods, the $j\omega$-axis is not transformed to the unit circle of the $z$-plane, and the left half-plane is not mapped to the region $|z|<1$. In case of the forward Euler method, the $j\omega$-axis is mapped to the line $\textrm{Re}\{z\}=1$, and the left half-plane is mapped to the half-plane $\textrm{Re}\{z\}<1$. This means that a stable analog prototype filter might be transformed to an unstable discrete-time filter. The backward Euler method transforms the $j\omega$-axis to the circle $|z-\frac12|=\frac12$, and the left half-plane is mapped to the region inside that circle: $|z-\frac12|<\frac12$. Consequently, stability and the minimum phase property are preserved, but the analog frequency axis is not mapped to the unit circle of the complex plane, so the resulting frequency response is not just a warped version of the original analog frequency response.

Note that there are also other very common methods for converting analog to digital systems, which are not based on direct mappings from the $s$-plane to the $z$-plane. These methods are based on time domain criteria. The two most common of these methods are the impulse-invariant method and the step-invariant method. They preserve the shape of the impulse (or step) response of the analog system by defining the desired digital system as the one whose impulse (step) response is a sampled version of the corresponding analog impulse (step) response. These methods introduce aliasing in the frequency domain, but - unlike the bilinear transform - they do not introduce frequency warping.


Properties of Spectral Transformations - Allocation (decomposition into even and odd part) I am trying to understand the Allocation property of Spectral Transformations. I can't.

I know that every function can be separated into an even part and into an odd part.

My problem is understanding proof for that:

$$s(t)=s_e(t)+s_o(t)$$
where:

$s_e(t)$ - even part of function
$s_o(t)$ - odd part of function

The problem is that in next step at every site I looked so far, it goes something like this:

$$s_e(t) =\dfrac{s(t)+s(-t)}{2}$$

Why did they do this? Where did $\frac{1}{2}$ come from, and why is $s(t)+s(-t)$

Can you please explain it to me, or point me to a good resource where I could read more about it?

Thanks!!!

Let's say that your signal is composed of two parts: even and odd:

$$s(t)=s_e(t)+s_o(t)$$

We also know following properties of this type of functions:

Even: $f(-x)=f(x)$
Odd:  $f(-x)=-f(x)$

Let's calculate the time inversion of your signal $s(-t)$ and apply above properties:

$$s(-t)=s_e(-t)+s_o(-t)=s_e(t)-s_o(t) $$

So now let's do the trick and add the: $s(t)$ and $s(-t)$:

$$\require{cancel} s(t)+s(-t) = \color{blue}{s_e(t)}+ \cancel{\color{red}{s_o(t)}} + \color{blue}{s_e(t)}-\cancel{\color{red}{s_o(t)}}=\color{blue}{2s_e(t)} $$

Solve it for $\color{blue}{s_e(t)}$, and you will get:

$$\boxed{\color{blue}{s_e(t)}=\dfrac{s(t)+s(-t)}{2}} $$

In the end let's subtract $s(t)$ and $s(-t)$:

$$s(t)-s(-t) = \cancel{\color{blue}{s_e(t)}} + \color{red}{s_o(t)} - \cancel{\color{blue}{s_e(t)}} + \color{red}{s_o(t)}=\color{red}{2s_o(t)} $$

Rearrange and you will get:

$$\boxed{\color{red}{s_o(t)}=\dfrac{s(t)-s(-t)}{2}} $$


Fit data to Gaussian distribution I want some data to fit the corresponding Gaussian distribution.

The data is meant to be Gaussian already, but for some filtering reasons, they will not perfectly match the prescribed and expected Gaussian distribution. I therefore aim to reduce the existing scatter between data and desired distribution.

For example, my data fit the Gaussian distribution as follows:

The approximation is already decent, but I really want to crunch the still tangible scatter between simulated data and expected distribution.

How can I achieve this?

The least-squares fitting is probably what your are looking for as this is maximum likelihood estimator for data points with Gaussian distribution around best fit curve.

There is already a built-in function lsqcurvefit in MATLAB:

I wrote a tutorial article on both linear and non-linear least squares methods (Gauss-Newton, Levenberg-Marquardt) if you need some deeper understanding of the topic:

Linear and Nonlinear Least-Squares with Math.NET

Please note you are asking how to use some actual piece of software (i.e. MATLAB) while this site is more about principles and practices. You can find hints for MATLAB on other forum (e.g. StackOverflow).


How can I interpolate between two FFTs? I have two FFTs that represent two impulse responses. I want to interpolate between the two of them to generate a 3rd FFT in the manner shown in the diagram below. This is to be used in real time audio rendering filtering between two FFT filters.

Linear interpolation between each of the FFT's complex coefficients doesn't work as wanted.

FFT plots of two impulse responses
(x = frequency/coefficient, y = frequency magnitude)

+                                   +
|                                   |
|        +                          |
|        |                          |
|        |                          |
|        |                          |
|        |                          |
|        |                          |                   +
|        |                          |                   |
|        |                          |                   |
+--------+------------------+       +-------------------+-------+
           Input A                              Input B

FFT plots of interpolation between A and B (Linear, and Wanted)
(x = frequency/coefficient, y = frequency magnitude)

+                                   +
|                                   |
|                                   |
|                                   |
|                                   |
|                                   |             +
|        +                          |             |
|        |                          |             |
|        |          +               |             |
|        |          |               |             |
+--------+----------+-------+       +-------------+-------------+
Linear Interpolation (A+B/2)            Wanted Interpolation

How do I do the interpolation as described in the Wanted Interpolation graph?

This is the solution for a linear interpolation between two frequencies as a "weighted average". You would just need to decide what your weighting actually is in terms of magnitude or dB scale for the vertical axis.  X is the interpolated frequency, weighted toward the stronger signal, and Y is the linear interpolation at the new frequency location of the signal strenths a and b.

Here is a simple example demonstrating the results that would be achieved with the formulas given:


How to Count the number of People in a room using CCTV Camera feed? I am trying to find out how can I count the number of people in a room using the feed coming in from a CCTV camera fixed in one of the top corners of the room. The number of cameras can be 2 at max but 1 for now. The objective is to only be able to tell how many persons are present in the room, in real time. The precision of it is not too important, for e.g. the system could suggest a range like "there are more than 5 people right now".

Upon browsing and researching, I came across a lot of very complex algorithms and terms related to Computer Vision and Image Processing that I have not understood at all like OpenCV, HOG, Blob Detection/Tracking, Foreground/Background Subtractions etc.

I am a design student. I do not have a background in programming. But, about a week back I got introduced to Processing. So, if some could help me out here with some tutorials/references/code to have it work in a language like Processing for prototyping purposes.

Thanks in advance.

This is clearly not an easy task.

The problem is, if you want a more-or-less accurate count, then you need to turn to advanced algorithms (and maybe use 2 cameras, or a kinect).

If you can't afford to take this path, then you need to try simpler options. Personally, I would try the following:

detecting skin pixels, segmenting the image with respect to that, then count the blobs (complexity medium);
detect faces, then count them. Obviously, people have to be facing the camera, but you can find any implementations of face detection with different libraries (complexity low);
use background subtraction (after some careful tuning to take into account natural illumination changes) and guess a count based on the number of pixels that belong to the foreground (complexity low, but expect poor results);
apply some existing pedestrian detection algorithm (complexity medium to high).

You don't have to implement all that by yourself. In your case, checking what OpenCV can do (or if your algorithm of choice is implemented in OpenCV) seems a good choice because:

Processing and openFrameworks have wrappers to OpenCV;
OpenCV usually has a good runtime performance;
many computer vision algorithms that will interest you are developed with OpenCV (face detection, background subtraction, probably pedestrian detection).


Find the difference equation and draw the simulation diagram Calculate the difference equation and then draw the simulation diagram of the below transfer function.
  $$
H(z) = \frac{Y(z)}{X(z)} = \frac{0.4142 + 0.4142z^{-1}}{1.4142 - 0.5858z^{-1}}
$$

I performed the normal procedure to find the difference equation, by cross multiplying and using the delay property of the $\mathcal Z$-transforms, I finally ended up with:

$$y[n] = 0.2929 x[n] + 0.2929 x[n-1] + 0.4137 y[n-1]$$

How do I draw the simulation diagram?

I want to know how to draw the block diagram of the difference equation.

I will try to draw the diagram but I am not sure if it is correct.

I don't require the answer but guidance on how to derive the diagram and final answer.

This is what I have done,

I have attached my solution.

is this correct?

I have modified my answer again.

Your equation is:
$$y[n] = 0.2929 x[n] + 0.2929 x[n−1] + 0.4137 y[n−1]$$
which has input $x[n]$ and output $y[n]$.

You also need $x[n-1]$, so there is a delay block required for that and $y[n-1]$ so a delay block is required to remember the last output.

You also need constant multipliers for the coefficients 0.2929 and 0.4137.

Read up on signal flow graphs.


Fourier coefficients of sum of two functions with different fundamental periods? If we assume $\quad x(t)\leftrightarrow a_k\:$ and it is periodic with fundamental period T.
How can we determine the fourier coefficients of the sum
$x(t-7)+x(-2t+3)$
I know that $x(t-7)\leftrightarrow e^{-jk7\frac{2\pi}{T}}a_k$
But, since these two functions have different fundamental periods we cannot apply linear combination rule. So how can we find the related coefficients? Thank you.
EDIT: After Stanley Pawlukiewicz's hints
$x(t-7)=\sum\limits_{k=-\infty}^{\infty}a_ke^{jk\dfrac{2\pi t}{T}}e^{-7jk\dfrac{2\pi}{T}}$
$x(-2t+3)=\sum\limits_{k=-\infty}^{\infty}a_{-k}e^{jk\dfrac{4\pi t}{T}}e^{-3jk\dfrac{4\pi}{T}}$
This new function is periodic with $\:T_0=T\:$ since $\:LCM(T,T/2)=T$ Thus,\begin{aligned}
\bar{a_k} = \left\{ \begin{array}{cc}
                a_k e^{-7jk\frac{2\pi}{T}}+a_{\frac{-k}{2}}e^{-3jk\frac{2\pi }{T}} & \hspace{3mm} \text{even}\;k \\
                a_k e^{-7jk\frac{2\pi}{T}} & \hspace{3mm} \text{odd}\;k \\
                \end{array} \right.
\end{aligned}

Since this looks like a homework and solving people's homework is outside the scope of this group, here are some hints

To have a Fourier Series, a function needs to be periodic.
$$ g(t)=g(t+T_1)  \quad  \Rightarrow  \sum_k \left(a_k \sin(2\pi k \frac{t}{T_1}) + b_k \cos (2\pi k \frac{t}{T_1}) \right) $$
and
$$ h(t)=h(t+T_2)  \quad  \Rightarrow  \sum_k \left(c_k \sin(2\pi k \frac{t}{T_2}) + d_k \cos (2\pi k \frac{t}{T_2}) \right) $$
If  $z(t)=g(t)+h(t)$
Under what circumstances of $T_1$ and $T_2$ would
$$ z(t)=z(t+T_3) \quad \text{?}$$
an what circumstances of $T_1$ and $T_2$  would
$$
z(t)\ne z(t+T_3) \quad \text{for any} \quad T_3
$$

and if there were circumstance where $z(t)=z(t+T_3)$,

Could we re express $g(t)$ and $h(t)$ as Fourier Series periodic in $T_3$?


How do I determine if a discrete signal is periodic or not? I want to know how I can determine whether a series of data is periodic or not.

I want to use Fourier transform/series. My data looks either aperiodic

  [111100001111000110010101010000101]

or periodic

  [11001100110011001100]

and I need to decide which it is automatically. What types of analyses or calculations can I perform to determine if a signal is periodic or not?

I would do a normalized autocorrelation to determine periodicity.  If it is periodic with period $P$ you should see peaks at every $P$ samples in the result.  A normalized result of "1" implies perfect periodicity, "0" implies no periodicity at all at that period, and values in between imply imperfect periodicity.  Subtract the data sequence's mean from the data sequence before doing the autocorrelation because it will bias the results.

The peaks will tend to diminish the farther away from the center they get simply because of having fewer overlapping samples.  You can mitigate that effect by multiplying the results by the inverse of the percentage of overlapping samples.

$$
U(n) = A(n) * \frac{N}{\lvert N-n\rvert}
$$
where $U(n)$ is the un-biased autocorrelation, $A(n)$ is the normalized autocorrelation, $n$ is the offset, and $N$ is the number of samples in the data sequence that you are checking for periodicity.

EDIT: This is an example of how to tell if the sequences are periodic.  The following is Matlab code.

s1 = [1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 1 1 0 0 1 0 1 0 1 0 1 0 0 0 0 1 0 1];
s1n = s1 - mean(s1);
plot(xcorr(s1n, 'unbiased'))

The "unbiased" parameter to the xcorr function tells it to do the scaling described in my equation above.  The auto-correlation is not normalized, though, which is why the peak at the center is around 0.25 instead of 1.  That doesn't matter, though, as long as we keep in mind that the center peak is perfect correlation.  We see that there are no other corresponding peaks except at the outermost edges.  Those do not matter because there there is only one sample overlapping, so that is not meaningful.

s2 = [1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0];
s2n = s2 - mean(s2);
plot(xcorr(s2n, 'unbiased'))

Here we see that the sequence is periodic because there are multiple unbiased autocorrelation peaks with the same magnitude as the center peak.


How to obtain the Spectrogram of a signal I am trying to use Matlab to compute the spectrogram of a signal but I can't seem to figure out how Matlab computes the spectrogram.Consider the case where I have a signal x(t) with M samples and I want to compute it's spectrogram using hann windows of length N and using an overlap of K samples.

Because there are M samples and the advance time is (N - K) samples, I would expect the resultant spectrogram to have dimensions of either $N$ by $\frac{W}{(N - K)}$ OR if matlab only provides frequencies up to $\frac{f_s}{2}$ then dimensions to be $\frac{N}{2}$ by $\frac{W}{(N-K)}$ but I am getting weird dimensions.

The code below shows the case where I have computed the spectrogram when M = 32 * 700 = 24400, N = 256 and K = 256 - 32 = 224.

%specgram(data, ChunkSize, SampleRate, Size Of Hanning Window, Overlap (in Samples)
    D = specgram( 1:1:(32*700),256, 8000 ,256 , 224);

produces a 129 * 693 matrix instead of the (256 * 700) or (128 * 700) spectrogram I was expecting.

Result should be an ($M$ by $\text{Number Of STFT Chuncks}$) matrix and
$$
  \text{Number Of Chunks} = \frac{\text{Number Of Samples in Signal}}{\text{Advance Time (in Samples)}}
$$

$$
\text{Number Of Chunks} = \frac{32 * 700}{32} = 700
$$

Have I missed out on something? Any help will be appreciated.

Since the input signal is real, half of the FFT data is "useless" because it is simply the complex conjugate of the other half. More precisely, a 256-long FFT of a real signal will give you: the DC amplitude, 127 amplitudes, the amplitude at Nyquist, and 127 conjugate amplitudes.

Among the 700 chunks, 256 / 32 - 1 = 7 extend outside the boundaries of the original signal.

Hence, 129x693.


Linear and Circular Convolution in Fourier Domain (DFT) Suppose we have two vectors A and B of length 100 and 80 obtained as a function of time. If we wish to perform convolution of the two vectors in the Fourier domain, we need to multiply the Fourier transforms of A and B. For that purpose, the length of B must be the same as A. The output would be of the same length as A or B after zero padding.

Where should we do zero padding in B, should this be at the beginning of B or at the end of B?
I was reading that convolution achieved via FFT is essentially a circular one. If our desire is to get a linear convolution, how do we ensure that the output is a linear convolution, i.e., how many end points should be rejected from both ends after inverse FFT?

the following matlab/octave code gives the linear convolution result using frequency domain :

A = ((-1).^[0:79]').*hamming(80);    % input one
B = blackman(100);   % input two

C1 = conv(A,B);     % A * B (convolution) in time domain
C2 = real( ifft( fft(A,179).*fft(B,179) ) ); % convolution using freq domain

The output will be identical of length 179 samples:


Transmission Line Information Loss Suppose we have a communication channel that has a bit-error rate of 50% and the bits are flipped at random such that the receiver has no efficient way to tell which bits have been flipped and there are no error correction bits transmitted with the signal. How do you calculate the information loss (in bits) imposed by the channel? I know there a Shannon formula for doing this but I don't know what it is! Could someone solve this for me with a specific example? Thank you.

A binary input - binary output channel which flips bits with equal probabilities is known as a Binary Symmetric Channel.

(Image Credit: Wikipedia - Binary Symmetric Channel)

The channel capacity $C$ is the maximum (actually, a supremum) mutual information between the channel inputs and the channel outputs over input distributions. You can't get information through a channel more than its capacity. Intuitively, this says that by designing the way you put inputs into the channel well, you can maximize the amount of stuff coming out of the channel that you can recover (and this limit is the capacity).

Mathematically, this is
$$C = \sup_{p(x)} I(X;Y)$$
where $p(x)$ is probability distributions on the input $X$ and $Y$ is the output of the channel.

In this case, you can write the input of the channel as $P(X=0) = q$ and $P(X=1) = 1-q$. With this distribution, you can see that $P(Y=0) = P(X=0 \text{ and no crossover}) + P(X = 1 \text { and crossover}) = q(1-p) + (1-q)p$ and $P(Y=1) = 1-P(Y=0)$. The mutual information between $X$ and $Y$ is given by $I(X;Y)=\sum_{x,y} p(x,y) \log \frac{p(x,y)}{p(x)p(y)}$ where $p(x,y)$ is the joint distribution of $x$ and $y$, and $p(x)$, $p(y)$ are the distributions of $x$ and $y$ respectively. By writing down $I(X;Y)$ and then maximizing over $q$, you will see $C = 1 - H(p)$ where $H(p)$ is the binary entropy function $-p \log_2 p - (1-p) \log_2 (1-p)$. Then, in your case, $p = \frac{1}{2}$ where $H(p) = 1$, so the capacity is $0$.

If you're interested in learning more about this, a few good books to look (at very different levels of difficulty):

Elements of Information Theory (Second Edition), by Tom Cover and Joy Thomas (Easy, Probably Best Starting Point)
Information Theory and Reliable Communication, by Robert Gallegar (A bit harder)
Information Theory: Coding Theorems for Discrete Memoryless Channels by Czisar and Korner (Pretty hard)
Information Theory by Robert B. Ash (Dover book, Easy)

And an unrelated book which I love, but couldn't keep myself from mentioning, is David Mackay's Information Theory, Inference and Learning Algorithms (free on his website), even though it doesn't go over channel capacity.


What is bit rate? How is it related to bandwidth? I am very confused about when to use bandwidth as bit rate and where not. Are these terms one and the same?

As hinted by Mr Duval, the term "bandwidth" has been adopted by networking engineers to mean something similar to "transfer rate" (just as "wideband" in networking tends to mean "high speed"). You've tagged your question "digital communications", though, where the term is very precisely defined.

Bit is the unit of information in information theory. Intuitively, a bit is the information conveyed by the answer to a yes/no question, where each answer is equally likely.

In communications we want to transfer information from a transmitter (information source) to a receiver (information sink). The information transfer rate is measured in bits per second. Informally, you can think of this as the number of zeros and ones that a system can transmit per second.

However, it turns out that information is transmitted and received using physical, analog signals (voltages or electromagnetic (radio) radiation). These signals pass through a physical medium, called a "channel", that acts like a filter: only certain frequencies make it from the channel input to its output. So, the channel is said to have a bandwidth, just like any filter.

One of the main jobs of communications engineers is to transfer information as fast and reliably as possible over a channel of given bandwidth. Basically, the problem is how to design a physical signal that:

Fits over the channel bandwidth
Conveys bits at a high rate
Does so reliably (with few bit errors at the receiver)
Does so at feasible complexity at both transmitter and receiver.

The design of such signals is, in a broad sense, covered by the term "modulation".

I hope this is clear -- please comment if you need clarification.


Phase Noise in Constellation Plot I have the following phase noise plot (-100dBc/Hz @ 1MHz offset and -130dBc/Hz noise floor):

The integrated phase jitter for this phase noise is around $4.6^\circ$. According to my understanding the spread of the constellation due to phase noise is equal to the integrated phase jitter.

But I see that the spread is much higher than the jitter. Why is it so? I am adding phase noise as follows:

modData = pskModulator(data);  % Modulate data
tx_out = modData.*exp(1i*phas_noise_signal);  % Add noise

The spread in the plot is the peak to peak phase, while the integrated phase noise is an rms quantity- The OP in his comment below this answer reported a measurement of 35° in the constellation and a 4.6° "integrated phase jitter". If this was rms phase then the measured peak as a ratio to the standard deviation would be $35\sigma/(2 x 4.6) = 3.8\sigma$ which is not an unreasonable result.

I had first assumed the OP was comparing an rms cycle-cycle jitter measurement in units of radians to the peak-peak phase deviation as shown on the constellation, hence the more detailed descriptions below that I will leave for future reference.

Determining RMS Phase

The rms phase error over the measurement duration can be determined from the SSB phase noise measurement $\mathscr{L}_\phi(f)$  using

$$\phi_{rms} = \sqrt{\int_{f=f_{low}}^{f_{BW}} S_\phi(f)df}$$

Where:

$S_\phi(f) = 2 \mathscr{L}\phi(f)$

$\mathscr{L}\phi(f)$: SSB Phase Noise (the values from the OP's plot above after converting from $\tt dBc/Hz$ to $\tt rad^2/Hz$)

$S_\phi(f)$: DSB Phase Noise (accounting for both upper and lower sidebands of phase modulation).

Th conversion from $\tt dBc/Hz$ to $\tt rad^2/Hz$ is $10^{k/10}$ where k is the value in $\tt dBc/Hz$.

With the range of integration further detailed below.

Determining RMS Jitter

A cycle-cycle jitter jitter measurement includes an effective high-pass filtering of the phase noise prior to integration, and converting to units of time.

The rms cycle to cylce jitter can be determined from the SSB phase noise measurement $\mathscr{L}_\phi(f)$  using:

$$t_{rms} = \frac{T}{\pi}  \sqrt{\int_{f=f_{low}}^{f_{BW}} (1-\cos(2\pi f)) \mathscr{L}_\phi(f)df}$$

This is arrived at by doing the following:

High pass filter $S_\phi(f)$ and integrate and square root to get rms phase:

$$\phi_{rms} = \sqrt{\int_{f=f_{low}}^{f_{BW}} |H(f)|^2S_\phi(f)df}$$

Convert units of phase in radians to units of time in seconds using:

$$t_{rms} = \frac{\phi_{rms}}{2\pi } T$$

Where

$H(f)= 1-e^{-j2\pi f T}$: Frequency response of cycle to cycle subtraction

$|H(f)|^2= (1-e^{-j2\pi f T})(1-e^{+j2\pi f T}) = 2-2\cos(2\pi f)$

$T$: Cycle time in seconds

$f_{low} = 1/T_m$ in Hz

$T_m$: total time duration of the data used for the rms jitter measurement in seconds.

$f_{BW}$: Channel bandwidth (all filtering prior to jitter measurement) in Hz.

This is further explained as follows:

The measurement of cycle to cycle jitter is done by computing the rms of the difference in time from the expected zero crossing location between adjacent symbols.

This measurement of cycle to cycle jitter is a high pass process - the jitter in one cycle is subtracted from the jitter in the previous cycle so is the filtering process given by the subtraction of the phase error (given in units of time) in one cycle from the phase error in the previous cycle.

As a discrete process this is $H(z)= 1-z^{-1}$ followed by an rms computation (compute the standard deviation of the time difference). Since the phase noise measurement represents a continuous time process, I used instead a continuous time integration with the equivalent filter $H(s) = 1 - e^{-sT}$. This is a completely valid approach since selecting the time error at each symbol location is a decimation process of the underlying continuous time function (of error vs time), which folds in all the upper frequency noise components, which then get captured in the integration.

This is actually a comb filter response, with a zero at DC forming a highpass over the measurement bandwidth given by the symbol rate.  Above cutoff and the noise is  uncorrelated so adds in power for a +3 dB gain, and below cutoff the correlation increases by the inverse of the frequency providing a $1/(f_c-f)$ rejection of the phase noise.  The easy way to see this is to consider the very very low frequency components of phase noise —- for those components the time shift in adjacent symbols is the same so is rejected when we compute the subtraction of the zero crossing offset between those two symbols. Below shows an example frequency response of the cycle to cycle measurement assuming a symbol rate of 10 MHz. Individual frequencies components above cut-off will have strong correlation every time the phase cycles through $2\pi$, hence all the nulls in the actual "comb" filter response.

Our observation time $T_m$ is given by the duration of the data set is another high pass filter (in order to see all the way down to DC you would need to observe forever!). An exact equation would add this additional high pass filter as an additional 20 dB/decade roll-off toward zero with a cutoff at $1/T_m$ and then extend the integration lower, but I have found it to not be necessary in all practical applications, getting the same result by simply starting the integral at $f_{low} = 1/T_m$ (similar to an equivalent noise bandwidth computation where the resultant filtering response is equivalent to a brickwall filter). So a brick-wall highpass with cutoff $1/T_m$ is used to model the observation duration, and similarly the channel bandwidth $f_{BW}$ is reasonably approximated as a brickwall lowpass filter since the phase noise is typically an insignificant contribution at the larger frequency offsets defined by the channel.

In summary, process the double-side-band (DSB) phase noise ($S_\phi(f)$) which is double the SSB phase noise ($\mathscr{L}\phi(f))$ with the high pass function that also has the additional noise power gain of +3 dB, and then integrate the result for total integrated power. Take the square root of this to get the standard deviation and this should match jitter once you convert between units of seconds as in jitter and units of radians as in phase noise using $t = \phi T/2\pi $ with $T$ as your symbol duration.


Hydrophone sensitivity rating If I use a hydrophone with the sensitivity given as: "-200 dB re 1V/μPa".

Question: How should I interpret this? Is it the voltage given at -200 dB SPL or??

All the Best!

Hydrophones have a weird way of specifying sensitivity. $1V/ \mu Pa$ simply means $10^6 V/Pa$ which is frankly an idiotic unit.

You scale this by -200dB and  you get $ S = 10^{-4} V/Pa$ or $S = 0.1mV/Pa$ which is a much more reasonable representation.

This has NOTHING to do with dBSPL which is referenced to the human threshold of hearing. There is no reference needed here, since you are using Pascals directly.


Discrepancy in stability conditions when calculating via RH criterion and Nyquist criteria I have the following open loop transfer function for a unity feedback system.

$$G(s)=\frac{K(s+20)^2}{s^3}$$

1.When using RH criterion it can be easily proved that the closed loop transfer function is valid for $K>10$.

2.Trying with Nyquist criteria for stability :
$$|G(jw)_{w=w_{pc}}| \leq 1$$

Here $w_{pc}$ is the phase cross over frequency ($\angle G(jw)_{w=w_{pc}}=-180^\circ$). First solving for $w_pc$, I get the value as $20$ rad/sec and  $|G(jw)_{w=20}|=\frac {K}{10}$ thus, using it with the nyquist criteria gives me $K<10$.

I don't understand why there is a difference. What am I missing?

The full Nyquist criterion can be summarized to the equation $N=Z-P$, where $N$ is equal to the number of encirclements in the clockwise direction of the minus one point, $Z$ is the number of right half plane closed loop poles and $P$ is the number of open loop unstable poles. So for closed loop stability you want $Z=0$, so $P$ counterclockwise encirclements of the minus one point. If the open loop is stable and you want the closed loop the be stable as well, then this simplifies to that you always keep the minus one point on your left, or in other words rotate around the minus one point in counterclockwise direction.

It can be noted that scaling the open loop transfer function by a constant $K$ then the Nyquist plot also gets scaled by this constant. When determining the closed loop stability this is equivalent to scaling the minus one point by $1/K$.


How do can we cancel out noises while listening voice on headphones? It's a common problem that while talking with the another person on Smart phones with the help of regular headphones,we get very annoying sounds at background . The opposite person's voice gets mixed with full of such noise signals to such an extent that sometimes it becomes literally impossible to continue the conversation.

Which kind of noises could be possible in such situation?

Also,how do can we convert a regular (ordinary) headphone into a headphone which has ability cancel out or isolate these noise signals to a maximum extent while listening voice on headphones and make conversation noise free?

Note:please consider any practical situation of your choice like opposite person(person on other) side is in silence zone area and you are walking on the road or travelling by the car or you are also in the silence zone.

here is a great website on building noise cancelling headphones http://headwize.com/?page_id=684

The concept behind noise cancelling headphones is simple. You listen to the environment and then play it back 180 degrees out of phase. When the original wave and the shifted wave meet, they cancel eachother out

The circuit is pretty simple. You have microphone go into an inverting amplifier, then you play it through the speaker (simplified version) This achieves noise cancellation.

edit

one major benefit of this solution is it is purely implemented in (analog) circuits. It doesn't require adaptive filtering, a decision engine, or any computing power. This type of circuit is what most commercial products use: Active filtering

Sometimes the solution to the problem isn't what we expected. We can get into signal estimation, adaptive filtering, DSP boards or MCUs, or any number of complex mathematically rigorous constructs but really, all that's needed is an inverting amplifier. If you want a software solution please state so (and in that even we can go back to signal estimation and adaptive filters) but since you asked for practical solutions to actually convert regular headphones I have presented a hardware implementation.

edit2

unfortunately the original website seems to have vanished, but here is the same exact content in a PDF I include the full article because it has parts lists, some helpful troubleshooting tips, and far more detail than I've included. And in case the new link breaks here is the main circuit with a brief explanation from pdf. The circuit itself looks complex but its really a 3 stage amplifier but one for each audio channel (left and right) so 6 amps total. As you can see the second stage is the inverting amp (makes signal 180 degrees out of phase) achieve our noise cancellation.

The following material is an excerpt from the noise cancellation headphone guide written by Jules Ryckebusch  and all credit goes to Jules for design and explanation.


GPS C/A Code Time Delay Using Ultrasound I am undertaking a project which requires to measure a distance from source to receiver. I am doing this using two ultrasonic transducers operating at 40kHz and using FFT search to determine the time delay and from that the distance. These are the following steps I'm taking to achieve this. All the processing is done in MATLAB.

Transmitting a C/A code at 1kHz using a carrier frequency of 40kHz
At the receiving end I sample the incoming signal at 160kHz
Generate a C/A code locally with a sampling frequency (fs) of 160kHz
Multiply the received signal with the one generated in step (3)
Multiply the signal from step (4) with sin(2*pifct) and
cos(2*pifct), where fc=40kHz and t=0:1/fs:1.023-(1/fs) to obtain the
I and Q components down-mixed to 40khz.
Then, compute the fft down-mixed signal, using fft(I+iQ), where
i=imaginary
Compute the fft of the localy generated C/A code and then take the
conjugate of that: conj(fft(C/A))
Compute the ifft,square and take the absolute magnitude of the result,
as so: abs(ifft(step7*step6)).ˆ2 to get the time delat

Are the above steps correct? Is there anything else I have to do or am I missing something? It works fine in matlab but when I use 'real' data the correlation is much weaker. I haven't actually used the ultrasonic yet I'm simply using one microcontroller to generate and send signal and another to sample it (connected with a wire). When I do the processing in matlab the results aren't as good as the pure simulation.
Note: DSP is very new to me, I'm just stuck with this project and learning as I go along.
C/A Code generator
[taken from][1]

tap=[2 6; 3 7; 4 8; 5 9; 1 9; 2 10; 1 8; 2 9; 3 10; 2 3; 3 4;
 5 6; 6 7; 7 8; 8 9; 9 10; 1 4; 2 5; 3 6; 4 7; 5 8; 6 9; 1 3;
 4 6; 5 7; 6 8; 7 9; 8 10; 1 6; 2 7; 3 8; 4 9; 5 10;4 10;
 1 7; 2 8; 4 10];

G1 = -1*ones(1,10);
G2 = G1;

s1 = tap(n,1);
s2 = tap(n,2);

tmp = 0;

for i = 1:1023
    G(i) = G2(s1)*G2(s2)*G1(10);
    tmp = G1(1);
    G1(1) = G1(3)*G1(10);
    G1(2:10) = [tmp G1(2:9)];

    tmp = G2(1);
    G2(1) = G2(2)*G2(3)*G2(6)*G2(8)*G2(9)*G2(10);
    G2(2:10) = [tmp G2(2:9)];
end

%resample
x=1:n_samples
CAa(x) = G(ceil(x*chiprate/fs));
end

Signal received from the "ADC"

  function send = s()
clear all;
%        CAa = ca(n,chiprate,fs,n_samples)
CAtest = CAa(1,1e3,160e3,(160e3*1.023)); %generate the outgoing signal
Fs=160e3;
n_data = 160e3*1.023;
t =(0:(n_data-1))/Fs
carrier = cos(2*pi*40e3*t); %carrier @40khz.
send = carrier.*CAtest %modulate

end

Code for xcorrr

close all;
modulated =  send();                          %generate incoming signal
%modulated = circshift(modulated,80000,2);    %simulate a time delay
fc = 40e3;                                    %carrier frequency
Fs = 160e3;
Ts = 1/(Fs);
n_samples = Fs*1.023;                         %sampling period
t = (0:(n_samples-1))/Fs;                     %generate time vector

b = CAa(1,1e3,160e3,(160e3*1.023));           %generate a C/A code locally

rx = modulated .* exp(-j*2*pi*fc.*t);            %downconvert to baseband

corrOut = abs(ifft(fft(rx).*conj(fft(b)))).^2;
figure
plot(corrOut);

Instead of using "square" for your carrier generation, use "cos".

Also the -1 square in your receiver will not generate a quadrature signal from + square as you have used, copied below:

This does not produce a sine and cosine:

x2 = square(2*pi*fc*t); %sin
x3 = -1*square(2*pi*fc*t);%cos

Here as well sin and cos directly would be more appropriate, but the following is how I would approach your receiver, using exp instead of sines and cosines since $e^{-j\omega t} = cos(\omega t)-jsin(\omega t)$:

modulated =  send();  %generate incoming signal at carrier fc
rx = modulated * exp(-j*2*pi*fc*t)   %downconvert to baseband
%the above will have a high frequency image at -2fc that must be filtered
coeff= firls(31, [0 .2 .5 1], [1 1 0 0]);
rxfilt= filter(coeff,1,rx)
corrOut = ifft(fft(rxfilt).*conj(fft(CA_local));

The $e^{-j\omega t}$ will shift the positve half of your modulated signal to baseband (0) and the negative half to $-2f_c$, so the example filter shown is to filter out the image. (Note there was no consideration given for actual filter design, this was quickly "pulled out of the hat" just to demonstrate what will work).


Locating the outer boundary of the iris I need help with iris localization. I have applied the imfincircle function to locate it, but it only locates the inner boundary of the pupil and does not find outer boundary.

How do I locate the outer boundary?

Hough transform is an edge based operation. For this reason, if you could extract the edges of the outer boundary (together with many noisy pixels) and adjust the radius parameter of the Hough transform, you will be finding the circle of the iris, but not anywhere else. You might want to check:

[centers,radii] = imfindcircles(A,radiusRange)

radiusRange is the parameter you should get right.

However, sometimes for some noisy images, Circular Hough Transform might not be the solution. In such cases you might try Hough ellipses or other methods for locating geometric primitives.


Low pass filtering on short int (16 bit PCM) samples I am writing software for processing audio given as 16-bit PCM samples. The first stage of the processing involves calculating the energy (or total variation) in a certain frequency range (above a certain cut-off frequency).

What I am currently doing is subtracting the energy of the low-pass filtered signal from the energy of the original signal. I found out that a lot of the processing is dedicated to converting the integer sampels to a floating point representation.

So my question is, is there a technnique for filtering the integer samples without converting them to floating point?

Yes of course, you can apply the filter directly to the integer samples, using fixed point arithmetic.

For example, if you use a FIR filter with coefficients [1 / 3, 1, 1 / 2], and an 8-bit resolution for coefficients, your output will be:

out[n] = (85 * sample[n] + 256 * sample[n - 1] + 128 * sample[n - 2]) >> 8

Two things to take be careful about:

Coefficients quantization might cause at best slight changes of filter responses, at worst cause filter instabilities. What is your filter type and the value of its coefficients?
Overflow / Data types / Truncation problems. In the example above, out can exceed the range of a 16-bit integer, so you will have to do some clipping.


Hilbert transform too large to store (out of core processing) I have a signal (saved on disk) that I would like to take the Hilbert transform of, but it's too large to fit in memory (all at once). I would like to cut it into blocks and take the transform of each block. I'm not super familiar with the Hilbert transform and I was wondering if there was some sort of padding or windowing that needed to be done / take overlapping chunks.

I would use a linear phase FIR Hilbert transformer, and use block processing, such as the overlap-add method. That means that you partition the input signal into contiguous non-overlapping blocks and compute the convolution of each block with your filter impulse response. The results are then overlapped and added. Overlap occurs because the result of the convolution is longer than the input block. E.g., if the impulse response has length $N$ and  the blocks have length $L$, the result of the convolution is $N+L-1$, and the last $N-1$ output samples must be added to the first $N-1$ samples of the next output block (assuming that $N<L$, i.e., the output block only overlaps with the next block).

In short:

$$y[n]=\sum_ky_k[n-kL],\qquad y_k[n]=(x_k\star h)[n]$$

where $h[n]$ is the impulse response, $x_k[n]$ is the $k^{th}$ input block, $y_k[n]$ is the convolution of $h[n]$ with $x_k[n]$, and, finally, $y[n]$ is the output signal.


Octave: How to attenuate all frequencies except one? I am implementing a morse code decoder in Octave. It already works, but not on very noisy signals, like the following:

I had the following idea to reduce the noise:

Perform fft to find the "beep" frequency of the morse tone
Use a bandpass filter to attenuate all other frequencies.

Unfortunately, I apparently can't do step 2. After filtering, my signal looks like this:

This shows about the first two milliseconds of the signal.

My code looks like the following. Please tell me:

whether my basic idea is promising
how to improve my code in a way that it actually does what I want

[x, f_sampling]  = wavread(filename);
t = fft(x);
l = length(t);
magnitudes = (abs(t/l))(1:l/2+1);
f = f_sampling*(0:(l/2))/l;
[peaks, locations] = findpeaks(magnitudes);
[maximum,index] = max(peaks);
f_main = f(locations(index))

f_cutoff = [0.1 0.9]*2*f_main/(f_sampling)
[b, a] = cheby1(20, 1, f_cutoff);
 y = x / max(x);

y = filter(b, a, x);

Think of the morse code as an OOK (on-off keying) modulated signal corrupted by additive noise. The optimal linear receiver that can maximize the SNR is matched filter. So "the best" filter you can choose is a filter matched to the frequency of $1$s. After doing that, you have a hypothesis test problem. So you need to choose a threshold (again, the optimal threshold is something known).

If you decided to go that way, I can give you the details.

Before any details, your code appears to be just a pseudo code. It is really not clear what specific waveform you assume. So I try to give you a general answer.

How to do it:

1 - Find the frequency of the $1$s, using FFT or any other scheme. Let's say it is $f_0$.

2 - You should know the bit duration (i.e. the time each $1$ or $0$ takes). Assume it is $T$.

3 - If the waveform of $1$'s is $s(t)$, then your matched filter is $h(t)=s(T-t)$.

4 - Calculate $y(t)=\int_{0}^{t}r(\tau)h(t-\tau)d\tau$, where $r(t)$ is your noisy signal at each interval $T$. This will give you the output of matched filter. We will look at this output at $t=T$ in each timeslot. In other words, you will have discrete values $y_k$ corresponding to each timeslot.

5 - Now you need a threshold $\gamma$ to decide which part of $y_k$ corresponds to a $0$ and which one is $1$. There are a bunch of options, depending on what you want to achieve, the statistics of noise, and the prior probabilities of $0$s and $1$s. If you are not so strict about optimality, assume equiprobable bits. In AWGN you can assume you will expect equally likely $0$s and $1$s which simplifies your choice. So an option would be just $\gamma$ equal to half of the magnitude difference between $0$s and $1$s.


Why does slowing down audio introduce aliasing? I have an 8khz mono buffer of samples, and I am trying to manipulate the playback speed of it by modifying the number of samples in it.  When slowing down (a speedScale value of 0.5 for example) then the values are repeated-- a buffer of samples: [0.1, 0.2, 0.3..] becomes [0.1, 0.1, 0.2, 0.2, 0.3, 0.3 ..].  When I do this, although the speed has slowed down correctly, I am noticing aliasing happening.  I was under the impression that repeating samples like this would not require low pass filtering to avoid aliasing.  Am I doing something wrong?

Here's the code for slowing-down/speeding-up:

function process(buffer, speedScale) {
    size = ceil(buffer.size / speedScale);

    counter = 0.0;
    ratio = 1 * speedScale;
    processed = [];

    for (i = 0; i < size; i++) {
        index = floor(counter);
        processed[i] = buffer.samples[index];
        counter += ratio;
    }
    return processed;
}

That's because you are interpolating the signal. Take a look at this figure:

This is zero interpolation, i.e. inserting zeros between samples, but applies to your case too. The replicas at 2*pi are distorting your signal. You can see them in the third graph, at -2pi/L, +2pi/L. So, you need a low pass filter with a cut off frequency pi/L, after expander.

In your case, the expander's output dtft will be X(w*2) + exp(-jw)*X(w*2), assuming input signal's dtft X(w). So, your extrapolator output will look a bit different, but the idea is the same.


Autocorrelation to diagnose faults I'm attending a very practical course on signals and i have some doubts, i hope to receive answers in layman terms.

1) My prof said i can use the autocorrelation of the output of a process to diagnose faults. I'm trying to understand how. Maybe by calculating the autocorrelation in different $\Delta t$ and notice changes ?

2) Is it possible to calculate the autocorrelation in subsequent $\Delta t$ and then compare the Power Spectral Densities of each  $\Delta t$ to notice faults ?

3) Why autocorrelation is considered a statistical approach? To me, it's simply multiplications (and a division) operations applied to signal samples...

4) In which cases i'm forced to use autocorrelation instead of FFT to diagnose faults ?

I would assume by faults you mean breaks in a signal line such that a reflection occurs. This is indeed an application of the autocorrelation: You transmit a sequence down a transmission line. If there is any change in impendance in the the line (such a break or kink etc) then a portion of your signal will be reflected back according to the reflection coefficient, given as:

$$\rho = \frac{Z_L-Z_o}{Z_L+Z_o}$$

Where $Z_L$ is the impedance of the load and $Z_o$ is the impedance of the transmission line, and $\rho$ is a complex valued reflection coefficient with magnitude ranging from -1 to +1, and is the amount of signal reflected back toward the source at the load (from 0 to 100% and with any phase angle).

Before explaining the autocorrelation further and how you would use it, I think it would be important to know some basic transmission line theory and in particular the properties that are given by this reflection coefficient.

For example, if we had a transmission line that had an impedance of 50 ohms, and if it was terminated with a 50 ohm resistor at the load, the reflection at the load would be 0 as given by the numerator of the equation, meaning none of the signal was reflected (this is the maximum power transfer condition when we match impedance of load to the line and to the source).  In such a condition, if we looked into the properly terminated transmission line from the source, since no reflection comes back it would look no different to us than if we were looking directly at a 50 ohm resistor at the source (this is one way to define the impedance of the cable; an infinitely long cable will have the same impedance as a resistor of the same value-- a shorter cable will only also look like this if it is terminated with a resistor having the value of the cable's impedance.)

Observe these other interesting conditions about the reflection coefficient and see how much it can tell us about the load.

$$Z_L = Z_o \rightarrow \rho = 0$$

$$Z_L > Z_o \rightarrow 0 > \rho \ge 1$$

$$Z_L < Z_o \rightarrow -1 \le \rho < 0$$

If the load is an open, the reflection is 1 (100%) and in phase.

If the load is a short, the reflection is -1 and out of phase.

As the reflection propagates back down the line, it constructively and destructively adds with the forward signal from the source as the line itself has phase shift proportional to the position on the line and the frequency of the signal. (A constant delay has a negative linear phase versus frequency). The velocity of this reflection is the speed of light divided by the square root of the dielectric constant of the transmission line (that sentence may be confusing, bottom line the reflection goes at the speed of light, or slower if in a transmission line with a higher dielectric constant than air).  So we can measure the reflection at the source, but we can't yet tell if the load (or anywhere along the line where a discontinuity occurs, essentially becoming our new "load") was an open or a short (if it was one of the two). If the length of the line was exactly a quarter wave-length in propagation distance of our source signal's frequency, then the source signal will have shifted 90° by the time it reached the load, and any reflection will shift yet another 90° by the time it reaches our source for a total 180° phase shift. This is exactly how a quarter wave line translates the impedance of a source to an open. Similarly a quarter wave 50 ohm line will translate the impedance of a 100 ohm load to 25 ohms (these are called impedance transformers and we use these techniques to match different impedances to minimize reflections that otherwise distort our signals).

With all that said, the great thing we can do with the autocorrelation of special sequences is to accurately resolve our distance in propogation time to the mismatched load or fault (a cable RADAR). As far as special sequences, we need to use sequences that have desirable auto-correlation properties, ideally ones that have a strong correlation to themselves only when they are completely aligned and zero at all other time offsets. White noise has this property specifically, so that is what we seek: good approximations to white noise. Such sequences are called pseudo-random sequences given they appear random yet we know exactly what the sequence is (and we need to in order to correlate to it).  If a sequence didn't have this property, meaning it also correlated to delayed versions of itself, then this would reduce our ability to accurately measure range to a reflection given multiple delays will also create a response from the one reflection.

Thus to do this, you transmit a psuedo-random sequence down the transmission line, while measuring the reflection (a "directional coupler" is a useful low cost relatively simple device for doing this). Correlate the reflection to the sequence at all delay offsets. The resulting plot of correlation vs delay offset will indicate faults AND range to fault based on size of peaks and distance from the origin in this correlation function. Opens and shorts will have the strongest peak. If you don't have an open and short but a damaging kink such as to change the line impedance (which would cause signal distortion and other problems so is a worthy fault to find) then you will see a smaller reflection at the range of this fault along with another likely even smaller reflection from the load that isn't perfectly matched. With a sensitive measurement you will often see something from the load due to imperfect matches at the physical transition from the transmission line to the load. If your source also is not well matched to the line, then you can see multiple reflections as a strong reflection bounces back and forth between source and load, but each of these will be at progressively further offsets from the origin with progressively smaller magnitudes so it will be very clear what is going on, especially if we have an estimate of the actual length of the transmission line.

I don't really understand Question 3. Autocorrelation is indeed simply multiplies and accumulation applied to signal samples, and is indeed a statistical approach (not either/or). One very important aspect of this is when you add samples that have noise values given by an independent identically distributed random process (the background noise you will get when you receive your reflected signal from the transmission line), the standard deviation of the sum will go up at $\sqrt{N}$ where $N$ is the number of samples you are adding. Coherent samples (which is what occurs when the sequence you are testing for is in perfect time alignment with the received reflected signal) will go up at $N$, so thus you achieve a $\sqrt{N}$ increase in the signal to noise magnitude ratio (in dB to be clear $10Log10(N)$).

Consider the simple Barker Code sequence that does have this nice autocorrelation property we seek, and assume you transmit a +1 when the code is 1 and an -1 when the code is 0. When you multiply this with itself in time alignment and accumulate the result, it will add to 11. The standard deviation for the IID random process of any noise on each sample will however increase by only $\sqrt{11}$. I leave it as an exercise for you to to see how it only adds to 1 for any other rotational shift. Consider what would occur with practical sequences that can have much longer lengths given by $2^{k}-1$ where k is a positive integer.

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

To answer your last question, you can do all this with an FFT since the circular autocorrelation can be computed using the FFT:

$$XCORR = ifft(fft(a)(fft(b^*)))$$

Where $a$ is the perfect sequence as transmitted and $b$ is the reflected complex signal as received. The * represents the complex conjugate. The FFT is used but you would still be doing the process I described of using the autocorrelation property of the sequence. If you had another FFT approach that could give you the range and intensity of the fault please clarify what that is and we can further compare if you don't see it now yourself.


phase mod x freq mod equivalence do any of you know how to convert a FM index to a phase mod index equivalent considering they both have the same carrier and modulator frequency? Is it even possible?

This is in respect for sound synthesis and I'm talking about one sinusoid modulating another sinusoidal oscillator.

I'm being able to get very similar sounds with these parameters:

FM: carrier freq 200Hz, modulator freq 480Hz, index 400*

this means the carrier freq oscillates 400Hz up then back to original frequency and then down 400Hz and back 480 times a second

Phase Mod: carrier freq 200Hz, modulator freq 480Hz, index pi/4 **

** this index means the phase goes up 1/8 of a cycle and then back, down the same amount and back to the original point

FM and Phase Mod are clearly two independent processes, but I've seen people mention it can be equivalent, but I've never read anything that would tell me how to get the same results with both. So I don't know if it can really be done even if it sounds practically the same thing or if I just got lucky.

thanks

phase is the integral (w.r.t. time) of frequency.  and they are definitely not two independent processes.  back in, i dunno, the 30s, Armstrong implemented FM for the first time using AM and a synchronized carrier that was 90° out of phase with the AM.  so he was first accomplishing phase modulation (PM) out of AM (for which they had the existing technology - they didn't really have VCOs back in that day), then he implemented FM out of the PM by integrating the audio modulation signal and applying that to the PM modulation input.

this is all in a basic communications systems textbook like that of AB Carlson.


Why we need Laplacian for edge detection? I am studying some trivial computer vision processing techniques and I came across edge detection algorithms. IMO sharp changes in the gradient are enough indications to detect the edges in an image but most of the modern edge detection algorithms use Laplacian with second derivation to do the same task.

Q1. Are we using second derivative because the change in flux is a much better indicator for edges?
Q2. Can we use something other than Laplacian?

"Need" or "Want"?

Q1:  The second derivative avoids the problem of a gradually changing color being "greyish".

Q2:  Yes.

Here is one:

After advice about detecting focus quality of objects in a photo detected using YoloV3

It is based on how planar the color.  Very planar means no edge.

There are many, many other variations of these techniques possible.


How to compute key of a song? I'm looking for a precise way to compute the key of a song. BPM would be cool too. Analog is ok. I tried to approach it from the DFT angle, but got stuck and can't really put a finger on what's the issue. New perspective would allow me to get out of the mental hole I'm in.

Insert a tuner into the audio chain and take down the notes coming through. Use the circle of fifths to calculate the key as best as possible.


how to extract a radio signal features with python? I have a dataset of radio signals that i want to classify. So, I decided to extract features. After some googling, I have noticed that I have to calculate moments, cumulants, Kutosis and skewness to use them as features.

for the moments I found that there is a function called scipy.stats.moment/. But then, I found another code :

def moment(sample, p, q):
    m = np.mean((sample**(p-q)) * (np.conjugate(sample)**q))
    return m

what is the difference between the 2 functions? and do I calculate the cumulants basing on the moments?

for the skewness and Kurtosis what do they really represents?
Are there any other features that could help?

Some of the more cutting edge signal classification techniques utilize deep neural networks (DNNs):

https://arxiv.org/pdf/1712.04578.pdf

This particular paper utilizes a convolutional neural network (CNN) and a residual neural network.  The convolutional "front-end" extracts its own nonlinear features prior to connection to the densely connected layers.

DNNs are typically very computationally expensive, so your approach to feature extraction approach might meet your requirements.


Comparasion of two wav files in Matlab Hi I am working on one project and I have some trouble with comparison audio files. I have read a lot of articles about this problem but I couldn't find a good one working.

In my project I need to compare two audio events like gun shot, explosion, crying, breaking of glass, ... These files have max 5sec duration but it can be different. I store these in database and there will be more examples from one audio event. Users will load his .wav file (max 5sec) into program, and program shows all files in database where is similar with the uploaded audio event. What I need is recognize which audio event is it. For example if an user uploads some explosions, program will return all files including some kind of explosion. So I am not looking for exact match. I just need find out how similar is this file with these in database. I need it in some numbers , maybe in percentage, so I can tell the files which match 90% or more.

So what is the best technique to find this kind of similarity, which way is best for my case?

Maybe you are looking for something called "audio fingerprinting". This might employ DSP techniques, but in this form the question might be percieved as somewhat broad. There is no one correct algorithm of this kind, rather a bunch of good ideas that can be combined in different ways.


What do triangles mean in frequency aliasing depiction? How edges of the triangle get on both sides of fs? What is their purpose? Why fs in the center? (from Lyons Understanding DSP)

The sampling rate is 6 kHz, so all the spectral content of what you are sampling will appear in the frequency range of -3 kHz < f < 3 kHz. The unshaded triangle from 3 kHz to 9 kHz is where the actual spectral content is.

When you sample at 6 kHz the spectral content of the unshaded triangle will appear as the triangle in the shaded region. More specifically the frequency at 4 kHz will appear at -2kHz, and the frequency at 7 kHz will appear at 1 kHz. Another way of saying this is, if you sample at 7kHz tone at 6 kHz it looks exactly the same as a 1 kHz tone sampled at 6 kHz. That's why aliasing leads to ambiguity - you can't tell whether that 1 kHz tone actually came from a 1 kHz, or a 7 kHz tone.

So the unshaded triangle represents where the original spectral content is, and the triangle in the shaded portion is where that triangle gets mapped to (when sampling at 6 kHz).


filter long data streams with python I am generating a time-stream data set that needs to be filtered.  The data-set is pretty large (i.e. too large to be filtered in one go).  I'm trying to split the filtering operation up into chunks but I obviously get a spike at each chunk start.  Is there a correct way to split a large data set up into smaller ones and have a continuous result - i.e (pseudcode).

timestream # variable with say 1e9 data points
b,a = sig.butter(1, 0.16, 'lowpass')

Split the timestream up in to say 1000 "chunks"

f_t = []
for c in chunks:
   t_data = sig.lfilter(b,a,c)
   f_t.append(t_data)

The key is to use the "state" argument of the filtering commands. The idea is that filter commands such as lfilter can accept a vector specifying the inital state of the filter, and can also output their state at the end of the filter computation. In pseudocode:

let state = vector of zeros
let filtered_data = 0
for c in chunks
    y, state_new = filter(b, a, c, state)
    append y to filtered_data
    let state = state_new
end

See the docs for more details.


What is the BER for multilevel modulated signal Assume I want to transmit a low or high symbol which can be BPSK modulated. For ordinary BPSK, bit 0 can be mapped to 1+i0 and bit 1 mapped to -1+i0. In my question, each two bits are BPSK modulated and a level dependent so bit 01 is mapped to 2+0i and bits 00 is mapped to 1+0i. Also bits 10 is mapped to -1+0i and bits 11 is mapped to -2+0i. This is BPSK (not QPSK) signal since there are BPSK demodulator and level detector at the receiver in which each of them output one bit at symbol period.I want to know what demodulation scheme to use in order to detect the symbol level and hence calculating the BER.

This is not BPSK! BPSK is binary phase shift keying with binary meaning 2 levels. What the OP is describing is Pulse Amplitude Modulation where multiple bits are mapped to multiple levels. In this case there are 4 levels so it would be 4-PAM. (As Dilip points out in the comments that the levels in 4-PAM would typically be equally spaced while here the inner levels are spaced further than the outer levels, complicating the BER computation slightly as detailed below but still PAM regardless).

The demodulation scheme would similar to BPSK demodulation in that you would multiply the received signal with the recovered carrier and low pass filter. Decision would be made on the output by setting decision thresholds half way between each of the levels (assuming equiprobable data) and then the derivation for BER follows similar to BPSK (solving for the tail probabilities under Gaussian distributions if the noise is assumed to be AWGN as is typically assumed, noting the different tail contributions involved since the levels aren't equally distributed as Dilip has pointed out - that said the thresholds at ±1.5 will likely dominate setting the error rate with significantly smaller error contribution from the threshold at 0).

With equiprobable distribution of symbols and AWGN with standard deviation $\sigma$ this would result in the following Symbol Error Rate:

$$\text{SER} = \text{erfc}(0.5/(\sqrt{2}\sigma) + \frac{1}{2}\text{erfc}(1/(\sqrt{2}\sigma) $$

For smaller $\sigma$ values the second term would likely be insignificant. For larger $\sigma$ there would be additional tail contributions from the next adjacent threshold levels which would also significantly complicate the BER computation as detailed below. Thankfully under all practical applications these tails would certainly be insignificant and can be safely ignored.

To go from Symbol Error Rate to Bit Error Rate, consider how many bit errors occur when a symbol error occurs. Because the mapping done by the OP was properly Gray-coded (only one bit change between adjacent symbols), only one bit error would occur with each symbol error.


How Does a Convolution Can Be Expressed as a Matrix Multiplication (Matrix Form)? I know this question may not be very relevant to programming, but if I don't understand the theory behind image processing I'll never be able to implement something in practice.

If I got it right Gaussian filters are convolved with an image for noise reduction since they compute a weighed average of a pixel's neighborhood and they are very useful in edge-detection, since you can apply a blur and derive the image at the same time by simply convolving with the derivative of a Gaussian function.

But can anyone explain me, or give me some references on how are they computed?

E.g. Canny's edge detector talks about a 5x5 Gaussian filter, but how did they get those particular numbers? And how did they go from a continuous convolution to a Matrix multiplication?

For this operation to work, you need to imagine that your image is reshaped as a vector.
Then, this vector is multiplied on its left by the convolution matrix in order to obtain the blurred image.
Note that the result is also a vector the same size as the input, i.e., an image of the same size.
Each row of the convolution matrix corresponds to one pixel in the input image.
It contains the weight of the contributions of all the other pixels in the image to the blurred counterpart of the considered pixel.
Let's take an example: box blur of size $3 \times 3$ pixels on an image of size $6 \times 6$ pixels.
The reshaped image is a column of 36 elects, while the blur matrix has size $36 \times 36$.

Let's initialize this matrix to 0 everywhere.
Now, consider the pixel of coordinates $(i,j)$ in the input image (not on its border for simplicity). Its blurred counterpart is obtained by applying a weight of $1/9$ to itself and each of its neighbours at the positions $(i-1,j-1); (i-1,j), (i-1,j+1),\ldots,(i+1,j+1)$.
In the column vector, the pixel $(i,j)$ has the position $6*i + j$ (assuming lexicographic ordering). we report the weight $1/9$ in the $(6i+j)$-th line of the blur matrix.
Do the same with all other pixels.

A visual illustration of a closely related process (convolution + subtraction) can be found on my personal blog Computers Don't See (Yet) - Linearized model of LBPs.


Fourier series expansion of $\exp[-j2kA\sin(\omega t)]$ I've been going through this paper where we exploit the well known Fourier expansion of the model signal as shown in the image below.

I've never come across this well known fourier expansion before. Can someone point me towards right direction to fill in the missing steps in above image and also the use of Bessel functions got me confused a bit. Any references would be of great help. Thanks.

They're talking about the Fourier series expansion of the periodic function

$$f(\tau)=e^{-jx\sin(\tau)}=\sum_{n=-\infty}^{\infty}c_ne^{jn\tau}\tag{1}$$

which sometimes appears in the literature, e.g., when computing the spectrum of a frequency modulated sine wave. The complex Fourier coefficients in $(1)$ are given by

$$c_n=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(\tau)e^{-jn\tau}d\tau=\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{j(-n\tau-x\sin(\tau))}d\tau=J_{-n}(x)\tag{2}$$

where the last equality follows from the integral representation of the Bessel function of the first kind of order $n$.

Combining $(1)$ and $(2)$ we get

$$f(\tau)=\sum_{n=-\infty}^{\infty}J_{-n}(x)e^{jn\tau}\tag{3}$$

With $x=2k_0A$ and $\tau=\omega_Dt$ you obtain the expansion used in the paper you quoted.


Achievable rate for MIMO channel with estimated channel matrix In a MIMO system, if the receiver has perfect CSI, i.e. knows that channel matrix $\boldsymbol{H}$, then the achievable rate is:
$$ R = \log \left|\boldsymbol{I}+\frac{P}{N_RN_T}\boldsymbol{H}\boldsymbol{H}^H\right|\tag1$$
where $P$ is the transmit power and $N_R,N_T$ the number of receive and transmit antennas respectively.
What happens if I don't know the channel, but estimate it?
Will it simply be
$$ R = \log \left|\boldsymbol{I}+\frac{P}{N_RN_T}\boldsymbol{H}\hat{\boldsymbol{H}}^H\right|\tag2$$
where $\hat{\boldsymbol{H}}$ is the channel estimate?

I'm interested specifically in LMMSE estimation (and Gaussian channel).

I found the following relation which answers my question

$$R=\log\det \left(\boldsymbol{C}_{MMSE}^{-1}\right)$$

where $\boldsymbol{C}_{MMSE}$ is the covariance matrix of the error in the data $\boldsymbol{e}=\boldsymbol{x}-\hat{\boldsymbol{x}}$

I found the answer in the paper "An Iteratively Weighted MMSE Approach to Distributed Sum-Utility Maximization for a MIMO Interfering Broadcast Channel" but it is said to be well known.

Therefore, we need to assume some method for estimation the data using the estimated channel, for example, match-filter or MMSE estimator, and then calculate its error (as a function of the real and estimated channels).


What is so special about this filter coef = [1,2,1]? Can anyone tell what kind of filter is this? I can see it's smoothing/averaging the input signal, but why specifically with [1,2,1] instead of something else? What is so special about this coef=[1,2,1]? Thanks

# Python code
length = 1000
signal = np.random.rand(length)
output = np.zeros(length-2)
coef= np.array([1,2,1])
for i in range(length - 2):
   output[i] = np.sum(signal[i:i+3] * coef / 4)

It is a symmetric odd-sized FIR smoothing kernel, belonging to the class of Pascal or binomial filters that somehow sample a Gaussian kernel. Plus, its coefficients are simple dyadic integers, that can be implemented as bit-shifts 1/4 1/2 1/4. The coefficients sum to one, hence it is unit gain at DC.

In simpler word: (one of) the simplest real smoother preserving symmetric feature location. It's efficiency is however limited


How do I find the transfer function of a speaker by analysing the audio produced by it? I have the original audio signal(MATLAB GENERATED) and the recorded audio signal (anechoic conditions). The setup consists of a DAC with a power amp plugged into the speaker. The audio is recorded via a condenser mic. I need to realize the transfer function of the speaker even though the properties of the rest of the setup are unknown.

Passive speaker or active speaker? If it's passive, then it's easy. If it's active, then you need to pull the transducer out and treat it as a passive speaker.
If you have a signal generator, use it as a source to generate white noise, otherwise use Matlab script.
If you are looking for impedance response, put your speaker in a voltage divider circuit (in series with a known resistor). The impedance of the speaker can be found by FFT the received signal (magnitude response and phase response are both useful). In this case you do need two devices / a tricky audio setup that generates and records at the same time.
If you are looking for surface pressure response, face your speaker at your mic (as close to each other as possible) and then play the same signal in an anechoic environment (no reverberation, ideally somewhere outdoors on a quiet sunny day) and do the same processing to received mic signal. Convert received voltages into pressure with your microphone's sensitivity data.
Impedance response tells you the physical property of the transducer and pressure response is what actually matters.