Laplace transform : integral vs poles & zeros If Laplace transform is expressed as :
$$\int_{-\infty}^{+\infty} h(t)e^{-st}dt $$
with :
$$s = \sigma + j\omega$$
and $h(t)$ an impulse response expressed as :
$$h(t) = Ae^{-\sigma_0t}\cos(\omega_0t+\phi) = e^{-\sigma_0t}\cos(\omega_0t)$$ ($A=1$ and $\phi = 0$ for simplification, $h(t)=0$ if $t<0$)
Then, each vertical line (parallel to the imaginary axis) in the $s$ plane corresponds to the Fourier transform of $f(t) = h(t)e^{-\sigma t}$ for a fixed $\sigma$.
For $\sigma = -\sigma_0$, the decaying exponential of $h(t)$ is canceled and we get the Fourier transform* of $h(t) = \cos(\omega_0t)$, that is : diracs at $\omega_0$ and $-\omega_0$ (not accurate, see (*) just below), hence two poles : $-\sigma_0 + j\omega_0$ and $-\sigma_0 - j\omega_0$ as in the following picture (illustration only, poles not located correctly) :

Indeed, we can understand that :
(*)Please, note that the following is not accurate : since $h(t) = 0$ if $t<0$, we should use the unilateral Laplace transform, not bilateral ! So here we would get the unilateral Fourier transform of a sinusoid, not the bilateral (with diracs only) one ! To see what this would be, please see the link given at the end of the accepted answer
$$\int_{-\infty}^{+\infty} h(t)e^{-j\omega t}dt $$
$$= \int_{-\infty}^{+\infty} \cos(\omega_0t)e^{-j\omega t}dt$$
$$= \int_{-\infty}^{+\infty} \frac{e^{j\omega_0t}-e^{-j\omega_0t}}{2}e^{-j\omega t}dt$$
$$= \frac{1}{2}\int_{-\infty}^{+\infty} e^{j(\omega_0-\omega)t}-e^{-j(\omega_0+\omega)t}dt$$
If $\omega = \omega_0$ or $-\omega_0$, then the integral would blow up due to the $$\int_{-\infty}^{+\infty} e^0dt $$ member, hence the poles in the s plane.
So as shown in ch.32, p.24 of The Scientist and Engineer's Guide to DSP (see figures below), with Laplace transform we multiply $h(t)$ with $e^{-st}$ = $e^{-\sigma}e^{-j\omega}$, that is we multiply $h(t)$ with sinusoids that are either :

(a) Exponentially decaying ($\sigma$ > 0)
(b) Stable ($\sigma = 0$)
(c) Exponentially growing slower than our impulse response decay ($ -\sigma_0 < \sigma < 0$)
(d) Exponentially growing, compensating our impulse response decay ($\sigma = -\sigma_0$) : OK, as studied above.
(e) Exponentially growing quicker ($\sigma < - \sigma_0$ and $\sigma < 0$)

(letters correspond to pairs of points in the s plane shown in figures below, always at a fixed $\omega$ or $-\omega$ value)

I understand case d : since we cancel the exponential part, we get only the (unilateral !!) Fourier transform of a sinusoid. That is : infinite at $\omega_0$ and $-\omega_0$ hence the poles (though I don't know why we have a continuous function of omega with infinite values at $\omega_0$ and $-\omega_0$ instead of diracs as in the original Fourier transform of a sinusoid -> Because we use unilateral Laplace hence Fourier, see end of accepted answer !).
Case a, c and e are intuitive. In case a, we multiply $h(t)$ with a decaying exponential. The integral will be some finite complex value (for all values of $\sigma > 0$. In case c, we multiply by an exponential growing slower than the decaying exponential of $h(t)$, hence some finite complex value for the integral (for all values of $-\sigma_0 < \sigma < 0$). In case e, we multiply the $h(t)$ by an exponential that grows quicker than exponential of $h(t)$ decays : hence integral does not converge (for all values of $\sigma < -\sigma_0$).
But for case b, I can't get the intuition of why the integral would be zero as shown with the area under the curve (red in the above figures) ? In other words, I understand the vertical line in the s plane at $\sigma = -\sigma_0$, it is the Fourier transform of $h(t)e^{-\sigma_0 t}$ so it is Fourier transform of $h(t)$ once its exponential component is removed, hence 2 poles due to sinusoid. We get poles whenever $e^{-st}$ is identical (compensates) to the impulse response. But what would cause Fourier transform of $h(t)e^{-\sigma t}$ to be 0 at some $\omega$ ? For which $h(t)$ and how it would impact the area under the curve (integral) ?

The definition of the Laplace transform you're using is called the bilateral Laplace transform, which is less common than the unilateral Laplace transform. The difference between the two is that the first has a lower integration limit of $-\infty$ whereas the second has a lower limit of $0$. This difference becomes irrelevant if the signal under consideration is zero for $t<0$. This is the case in the example of the book. Note, however, that the bilateral Laplace transform of the impulse response $h(t)$ that you defined in your question doesn't exist for any value of $s$. It does exist if you set $h(t)$ to zero for $t<0$ (i.e., multiply it with a unit step $u(t)$).
The figures in the book refer to a causal notch filter. Note that the region of convergence (ROC) of the Laplace transform of that impulse response is to the right of the poles. Consequently, the Laplace transform evaluated at any fixed value of $s$ to the right of the poles will be finite, simply because we're inside the ROC, i.e., the integral converges. If we choose $s$ exactly at the zero of the filter then we evaluate the filter's response at the notch frequency, which must be zero, simply by the definition of "notch frequency". The filter's response to a signal at that frequency must be zero. Finally, if we choose the value of $s$ exactly at or to the left of the filter's poles, we are outside the ROC, in which case the integral doesn't converge.


In case of using a Doppler radar, do we have to use IQ demodulation to find the direction of a target? I want to know the way to find the velocity of a target by using a doppler radar.

I read some paper about it,

In summary,

send rf signal to a target
receive the reflected signal from the target.
conduct the IQ demodulation.
compute complex FFT to find the doppler frequency.

So,

If we want to know the direction of a target, or a velocity toward or away from the radar, we must use the IQ demodulation technic to find a negative velocity, right?

If it is right,

What does The IQ demodulation mean?

The reflected radar signal is just a rf signal.

Do we need to modulate a radar signal before transmitting it to a target?

Radar designer here: It sounds like you’re talking about  pulse-Doppler (PD) radar systems. For PD radars, the process is essentially as you described:

Generate a waveform (typically at IF) and then convert to RF
Transmit the waveform at RF
Receive the waveform at RF, and eventually mix it down to IF.
Apply IQ demodulation (for digital receivers a Hilbert transformer is typically used but that’s not an important detail). In this step, the received signal is added with a version that is 90 degrees out of phase with itself multiplied by the square root of -1, or the imaginary number. In general, you’ll have the generic signal s(t) = I(t) +jQ(t), where j is the symbol for the imaginary number. You can accomplish IQ demodulation in analog hardware, or you can use a digital receiver to accomplish it as well when you sample (via a Hilbert transform)
Collect the individual pulses and “stack” them in an array
Perform pulse compression and Doppler processing. Doppler processing is accomplished by computing the DFT across each pulse.

Since each waveform has phase which is recoverable with IQ data, the cross pulse DFT basically allows you to see the rate of change of that phase information, which is naturally the Doppler frequency. By being complex, yes you can tell whether a target has a negative or positive velocity, which is naturally very useful and pretty crucial for most tracking radars.

As for your question about needing to modulate the waveform: if by modulation you mean something other than a pure sinusoid, yes generally that’s a good idea for pulsed Doppler radars, as you’ll need good pulse compression. LFMs (linear frequency modulated) signals are a popular choice due to their nice properties, though other modulation schemes exist (NLFM, PSK, FSK, etc). If you just use a pure sinusoid, you effectively get no pulse compression, and you’ll have a pretty poor estimate of where targets may be, and the performance would be awful if you have several closely spaced targets.


How can I break a filter down into second order sections? This is a follow-up question of How do I apply a Chebishev filter? I designed a filter in Scipy:

import scipy.signal as signal
fs = 240
signal.cheby1(8, 0.05, [0.1/(fs/2), 10.0/(fs/2)],btype='band', analog=0, output='ba')
/usr/lib/python2.7/dist-packages/scipy/signal/filter_design.py:268: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless
  "results may be meaningless", BadCoefficients)
(array([ -3.64482771e-08,  -3.98755170e-23,   1.27568970e-07,
     3.98755170e-23,  -2.55137939e-07,  -2.09346464e-22,
     3.18922424e-07,   1.48536301e-21,  -2.55137939e-07,
    -2.79128619e-22,   1.27568970e-07,  -3.98755170e-23,
    -3.64482771e-08,   0.00000000e+00,   4.55603463e-09]),
array([  1.00000000e+00,  -1.53699059e+01,   1.10865293e+02,
    -4.98175837e+02,   1.56087517e+03,  -3.61582779e+03,
     6.40630344e+03,  -8.85525782e+03,   9.65121227e+03,
    -8.32140081e+03,   5.65721533e+03,  -3.00062986e+03,
     1.21729735e+03,  -3.65135628e+02,   7.63720484e+01,
    -9.95191544e+00,   6.08648759e-01]))

The first array should contain the b coefficients, the second array should contain the b coefficients of the filter.

The result of this filter is oscillating and seems to be wrong. Hilmar wrote that this filter "is a numerically challenging filter, since you have poles very very close to the unit circle" and I must "break the filter down into second order sections and apply those sequentially".

Since I have no background in signal processing and found nothing about that anywhere I wanted to ask here how I can do that.

edit: The zero-pole-gain ('zpk') representation is

signal.cheby1(8, 0.05, [0.1/(fs/2), 10.0/(fs/2)],btype='band', analog=0, output='zpk')
/usr/lib/python2.7/dist-packages/scipy/signal/filter_design.py:268: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless
"results may be meaningless", BadCoefficients)
(array([ 0.95041122+0.63504447j,  0.95041122-0.63504447j,
    0.77688699+0.32179713j,  0.77688699-0.32179713j,
    0.72152481+0.14352021j,  0.72152481-0.14352021j,
    0.70710678+0.j        , -0.95041122+0.63504447j,
   -0.95041122-0.63504447j, -0.77688699+0.32179713j,
   -0.77688699-0.32179713j, -0.72152481+0.14352021j,
   -0.72152481-0.14352021j, -0.70710678+0.j        ]), array([ 1.19049737+0.j        ,  1.16303520+0.1084667j ,
    1.16303520-0.1084667j ,  1.08874047+0.19530989j,
    1.08874047-0.19530989j,  0.99220878+0.25677213j,
    0.99220878-0.25677213j,  0.93226615+0.2765259j ,
    0.93226615-0.2765259j ,  0.87852354+0.22994371j,
    0.87852354-0.22994371j,  0.83333317+0.15993377j,
    0.83333317-0.15993377j,  0.80434290+0.08130729j,
    0.80434290-0.08130729j,  0.79450811+0.j        ]), -3.6448277062436296e-08)

Below the results as sos matrix and gain. Every row in the matrix is a single second order section in the coefficient order of b0, b1, b2, a0, a1, a2. Obviously we need a0=1, and the sections are normalized so that b0=1 as well. The overall gain accumulated is in the single number "gain". The fact that this gain is actually in the order of 10^-9 shows that this is a numerically tricky beast.

Here is how this works.

Design the filter in poles and zeros (zp or zpk), NOT in polynomial form (ba). I imagine that the cheby1() function has an output option for that. If not, you can't use it.
The poles come in conjugate complex pairs, the zeros in this case are real, but you can think of them as a conjugate complex pair with a phase of 0.
Group one pole pair each with one zero pair. Start with the pole pair that has the highest magnitude and find that zero pair that's closest to it. This step is really important!!
Turn the poles and zeros into polynomial in 1/z of the shape (1-p/z)*(1-p'/z). Let's say the complex zero/pole is x+j*y. Then b0/a0 coefficients are always 1, the b1/a1 coefficients are -2*x and the b2/a2 coefficients are x^2+y^2.
Repeat over all poles. Each conjugate pole pair results in one second order section.

A word of caution, if I may: If none of things that I have listed above mean anything to you, than you probably should not be doing this. This is not trivial and many things can go wrong in non-obvious ways. Depending on what you are planning to do with the results, you really want to make sure that your code is fully vetted and properly tested. Unless you have a good idea on how to test this, I would not base anything mission critical on the results.

    sos =
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000  -1.824520683146258   0.835012167790563
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000  -1.825351261746257   0.854128332434601
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000  -1.845749373866514   0.899792317177862
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000  -1.890243419701735   0.963689345198141
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000  -1.989176997083584   0.989218209927750
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000  -1.996496707908275   0.996511870590091
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000  -1.998714416931744   0.998722644911976
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000  -1.999653030269727   0.999659262704001
gain =
    4.556034632804534e-009


Source Localization with the BLUE Estimator I'm trying to implement BLUE estimator in MATLAB for source localization and after my research I've come up with a theoretical example in Steven Kay's "Fundamentals of Statistical Signal Processing: Estimation Theory" book (Example 6.3).

In this solution, while estimating $\theta$, the Gauss-Markov theorem is used and a point called as nominal source position($R_{n_i}$) is defined rather than real source position($R_i$). There is equation which shows the relation between $R_{n_i}$ and $R_i$ found by first order Taylor series expansion.

I'm confused about finding the value of $R_{n_i}$ variable, can anyone explain me this solution?

This could be a very simple question, but I want to understand the theory and my mind is not clear now. When you look at the solution in the book, I think you will understand better than I wrote.

EDIT:

N = 8; %total number of receivers
var = 1;
c = 3*10^5; % speed of electromagnetic waves
x = [0 8 7 1 5 9 3 1]; %xi receiver location
y = [0 1 9 9 2 5 8 4]; %yi receiver location
xs = 3; ys = 6; %real source location
xn = 6; yn = 10; %nominal source location
Sxs = xs-xn; Sys = ys-yn;
for i=1:N
    R(i) = sqrt((xs-x(i))^2*(ys-y(i))^2);
    Rn(i) = sqrt((xn-x(i))^2*(yn-y(i))^2);
    co(i) = (xn-x(i))/Rn(i); %cosai
    si(i) = (yn-y(i))/Rn(i); %sinai
    noise(i) = var*randn(1,1); %e1,e2,..
end
for i=1:N-1
    e(i) = 1/c*(co(i+1)-co(i))*Sxs+1/c*(si(i+1)-si(i))*Sys+noise(i+1)-noise(i); %linear model (6.26)
end
H = 1/c*[co(2)-co(1) si(2)-si(1);
         co(3)-co(2) si(3)-si(2);
         co(4)-co(3) si(4)-si(3);
         co(5)-co(4) si(5)-si(4);
         co(6)-co(5) si(6)-si(5);
         co(7)-co(6) si(7)-si(6);
         co(8)-co(7) si(8)-si(7)];

A = [-1 1 0 0 0 0 0 0;
     0 -1 1 0 0 0 0 0;
     0 0 -1 1 0 0 0 0;
     0 0 0 -1 1 0 0 0;
     0 0 0 0 -1 1 0 0;
     0 0 0 0 0 -1 1 0;
     0 0 0 0 0 0 -1 1];
inaaT = inv(A*transpose(A));
estTheta = inv(transpose(H)*inaaT*H)*transpose(H)*inaaT*transpose(e); % eq. 6.27

This is not a modular code, I just want to be sure about not missing any point. However, when I run this code I get useless estTheta values. I think again I leave out something in the theory.

One problem with your code is that the range calculations should be

R(i) = sqrt((xs-x(i))^2$\color{red}{\bf +}$(ys-y(i))^2);

Rn(i) = sqrt((xn-x(i))^2$\color{red}{\bf +}$(yn-y(i))^2);

But that doesn't appear to help.

The thing that seems to be the issue is that your noise variance is way too high.  If I set it to var = 0.00000001 then I get good estimates.

The reason is because you're using this variance (or a standard deviation of 0.0001) to add noise to your H matrix. Looking at the values there:

you can see you're adding something with a variance of 1 (standard deviation of 1), to something of the order of $10^{-5}$.  That means your SNR (signal-to-noise ratio) is about -100dB.  I know of no algorithm that will help you at SNRs that low.

You left out that Kay defines
$$
R_i(x_s,y_s) = \sqrt{(x_s - x_i)^2 + (y_s - y_i)^2}
$$
and
$$
R_{n_i}(x_n,y_n) = \sqrt{(x_n - x_i)^2 + (y_n - y_i)^2}
$$
where $x_s = x_n + \delta x_s$ and $y_s = y_n + \delta y_s$.

Then the aim is to write $R_i$ in terms of $R_{n_i}$.

Wikipedia has an example of Taylor series in two variables, and the first order approximation is just:

$$
f(x,y) \approx f(a,b) + (x-a) \frac{df}{dx} (a,b) + (y-b) \frac{df}{dy} (a,b) + \ldots
$$

Then all you need to do is apply Wikipedia's expression to Kay's to get

$$
R_i \approx R_{n_i} + \frac{x_n-x_i}{R_{n_i}} \delta x_s + \frac{y_n - y_i}{R_{n_i}} \delta y_s.
$$


Auditory receptive field mapping The Wikipedia page on sound localisation states:

  Localization accuracy is 1 degree for sources in front of the listener
  and 15 degrees for sources to the sides. Humans can discern interaural
  time differences of 10 microseconds or less.

However, I can't find any more detailed information on localisation accuracy.

What I'm after is a plot of accuracy (Y) against azimuth (X).

Can anyone get me one step towards this information?

π

The ultimate book on the topic is probably http://www.amazon.com/Spatial-Hearing-Revised-Psychophysics-Localization/dp/0262024136 . It's a bit dated but contains and excellent review of the main psycho-acoustic research on spatial hearing.

It is a very complicated topic. For example azimuth and elevation perception work completely different and hence azimuth resolution is quite a bit better than elevation resolution. The actually resolution (or jnd for just noticeable difference) depends heavily on the actual and the properties of the signal.


Power Spectrum Density and Frequency if i have some random signals (sampling rate = 10Hz, 0.1s per data)

Using python library i transformed it to power spectral density

power spectral density forms = f, psd (using mlab.psd)

I'm really curious about f....

If f is 0.01, is it really means 0.01Hz or 0.1Hz(10Hz * 0.01Hz) in signal??

I confused about this conception

Is it relevant to sampling rate and psd frequency?

I'm really thank you for your help!

If you are asking about the grid resolution of the PSD then it is given by:

$$ \frac{{F}_{s}}{N} $$

Where $ {F}_{s} $ is the sampling rate and $ N $ is the number of samples.
Pay attention that you can increase the grid resolution by Padding the data with zeros.
Though it will increase the grid resolution it won't increase the resolution of the DFT (In the sense of resolution to see to close, in frequency, harmonic signals).


Baseband vs Passband simulation in terms of noise impact From other posts I understood that it is normally no problem to simulate a communication system in the baseband though it is carrier modulated in reality. It just speeds up the simulation.

But what if I want to add noise to the signal and the noise shape is depending on the frequency? Of course I want to superpose my signal with the noise around my carrier frequency and not the baseband noise. I guess that simply adding noise and baseband signal would falsify the simulation...

Is there a solution so that I can still simulate in the in the baseband? Maybe to transfer the noise to baseband somehow?

In fact, bandpass noise (or any kind of noise process) can be modelled as a system $h(t)$. As so, there is a lowpass equivalent $h_l(t)$ such that

\begin{equation}
h(t) = Re[h_l(t)\exp(j2\pi f_0t)]
\end{equation}

So you can work with $h_l(t)$ in the baseband.


How to find the energy of a rect signal? Given a signal
$$
x(t)=-2\operatorname{rect}\left(\frac t4\right)
$$

What is the energy of signal $-x(t + 1)$?
What are the steps needed to find the energy?

Here's a definition of the energy of a signal from Wikipedia: https://en.wikipedia.org/wiki/Energy_(signal_processing)
For your specific signal, I recommend to draw the original signal and then the transformation. Finally, just calculate the squared magnitude and the area under it.
These are plots of the steps I mentioned (from Matlab):

As you can see in the last plot, the energy is basically the area under the rectangle, which is in this case 4*4=16.


How can I run this MATLAB code snippet for the waveform synthesis? Can you give me an example for this MATLAB snippet to generate sine wav with specified phase?

function synthesize_fp(file,f,d,p,gamma)
% Matlab function synthesize_fp(file,f,d,p,gamma)
% creates a .wav audio file of a sound where all frequencies,
% amplitudes(power) and phase may be specified.
%
% file is a string which is the name of the .wav file.
% f is a length n vector of frequencies in Hz
% d is the duration in seconds
% p is a length n vector of amplitudes
% gamma is a length n vector of phase shifts, as a fraction of the
%  period of the first harmonic f1.
%
% Mark R. Petersen, U. of Colorado Boulder Applied Math Dept, Feb 2004

Fs=22050; nbits=8;              % frequency and bit rate of wav file

t = linspace(1/Fs, d, d*Fs);    % time
y = zeros(1,Fs*d);              % initialize sound data
for j=1:length(p);
  y = y + p(j)*cos(2*pi*f(j)*(t-gamma(j)/f(1))); % sythesize waveform
end
y = .5*y/max(y);                % normalize.  Coefficent controls volume.
wavwrite( y, Fs, nbits, file)

This is a programming question (won't benefit anyone apart from you) and what's more, you are using program rather different from application you described. In order to generate the sinusoidal wave with given phase you don't need to use that. Above script is used mainly for synthesis of a wave with harmonic content depending on type of sound (flute, violin, etc), and personally I don't like it as it is terribly written using a for loop.

For general answer you can refer to Phase information. Additionally usage description is very clear - did you tried to follow it? For example more complicated sound would be defined as:

file = 'output.wav';
f = [100 200 300 400];
d = 5;
p = [1 0.5 0.25 0.124];
gamma = zeros(1, 4); % Start with no phase-shift

For pure cosinusoid use the:

file = 'output.wav';
f = 100
d = 5;
p = 1;
gamma = 0;

Application of phase shift is not very straightforward as its described by relative shift with respect to fundamental frequency. Therefore I rather suggest to write it on your own and keep it clean:

fs = 400; % sampling frequency
d = 1; % duration
f = 10; % wave frequency
phi = pi/2; % phase shift - will produce cosinusoid for pi/2
t = linspace(1/fs, d, d*fs); % time vector
y = sin(2*pi*f*t + phi);


Playing samples in a loop, without clicks or pops I'm trying to build an audio player that has this specific feature: when you scroll to x sample time, [x to x + 4096] samples are played in a loop until you move the mouse to another place.

The samples are taken from the backing audio. The problem is, I fail to play them back as a continuous note: the looping is not seamless.

So far I tried to start and end the loop at zero crossings. This slightly improved it, but I can still hear the seamless looping.

One idea I had was to go to the frequency domain, remove all, but the max magnitude bin and do an inverse FFT afterwards. This, sort of works, the looping is better, not perfect, but better, however, removing all the other freq components from the loop leaves it really "dry", not resembling the original (missing the harmonics)

Is there a way to handle this in the time domain? Where should I look? Are there any filters/algorithms that could help me?

The popping comes from discontinuities in the signal or its derivative.  Even if you match the starting sample value with the ending one, you will still have a "corner" in most cases.  Imagine a pure sine wave of one frequency at the start and a pure sine wave of another at the end. Discontinuities tend to produce a lot of high frequency noise.  Your second solution using the low pass filter eliminated that noise (it smoothed out the signal), but at the cost of eliminating the signal too.

The simplest solution to this is just to crossfade the end into the start.  Admittedly, at typical playback rates (22KHz) you are looping about 200ms which gives very little time to fade in/out.  At this point, you may be closer to doing an (overlapped) window function.  Doing an abrupt crossfade would be roughly akin to a Tukey window.


How does the PSK demodulator differentiates one carrier from the other in CDMA I'm trying to understand the way in which CDMA works. Once the transmitter bit sequence is spread it is modulated for instance using a PSK modulator. Lets assume BPSK is used. I'm particularly interested in the case in which there are 2 transmitters each with its own spread sequence, orthogonal to the other, so they don't mutually interfere . Assuming they both transmit to the same receiver at the same time using the same carrier frequency and BPSK modulation. How is the receiver capable of demodulating the BPSK signal in order to recover the bit sequence from both receivers so it can be passed to the respective de-spread stages where the original messages can be recovered?

To me the transmitters' information is divided in code, but codes can only be separated once they are bits not while they are "pass band symbols", right? And before they can be successfully recovered as bits, isn't there a lot of interference at the BPSK demodulating stage where phases from 2 BPSK carriers are being received at the same time? How can the BPSK demodulator demodulate useful bits under that circumstance?

To me the transmitters' information is divided in code, but codes can only be separated once they are bits not while they are "pass band symbols", right?

Wrong!

I'm not quite sure how to continue this answer, other than saying, "please re-read the CDMA introduction you've been reading".

The whole idea is that you receive a sequence that is a superposition of the two transmit signals, and by correlation with the known spreading sequences, you can extract the bits from that.

Let's put it more mathematical:

Your spreading sequences are called $c_i$.

Orthogonality is defined via the inner product:

$$a \perp b \iff <a, b> = 0 = <b,a>\text{ .}$$

Now, let's assume perfectly orthonormal sequences, so that $<c_i, c_j> = \delta_{ij}$, i.e. the inner product is 1 only if you multiply a sequence with itself, and 0 for any other sequence.

Corollary, the $c$ span a vector space $C$ which contains all linear combinations of CDMA sequences; its dimension is the same as the number of possible sequences, since these per definition are orthogonal (and hence, especially linearly independent) and form a basis.

Now, let's apply the code. You'll take one symbol $s$, in case of BPSK, one bit, represent it as $\{-1,+1\}=S$, and multiply it to get the transmit sequence $t$. You get the set of possible transmit sequences $t_{i,s}$ for user $i$:

$$T_i = \left\{sc_i \,\forall s\in S \right\}= \left\{ -c_i, +c_i\right\}$$

Notice that $T_i$ spans a vector space that has one dimension less than $C$!

Let's now define a receive signal $r$. We don't know what symbols were sent yet, so we say that user $i$ sent symbol $s_i$, and user $j$ sent symbol $s_j$:

$$r = s_ic_i+ s_jc_j$$

Now, we simply consider $r$ to be a sequence of complex numbers.

What happens now if we multiply $r$ with the individual spreading sequences $c$? Let's do it with user $i$'s code:

$$\begin{align}
<c_i,r> &= <c_i, s_ic_i+ s_jc_j>&\text{linearity of inner product}\\
&= <c_i, s_ic_i>+ <c_i, s_jc_j>&\text{$s_i$, $s_j$ are scalars}\\
&= s_i <c_i, c_i> + s_j <c_i, c_j>&\text{orthonormality}\\
&= s_i \cdot 1 + s_j \cdot 0\\
&= s_i
\end{align}$$

How nice! The same works for every other user, too.

I was a bit pushing the vector space notion here: It's important to realize that all this simply happens due to linear algebra representing our signals. You can now also see why CDMA systems that support a high number of users get slower per user (simply because that dictates a higher number of base vectors of $C$, which dictates a higher dimensionality and thus, longer sequences). Also, it shows that CDMA on this theoretical level works with any linear modulation, since we simply need $S\subset \mathbb C$, i.e. $s$ just needs to be scalar.


Confussion on roll-off factor, sampling time and Raised Cosine pulses. I am reviewing some bandlimited communications theory and something is bothering me -

My book says increasing the roll-off factor makes the received sequence more tolerant to timing recovery errors but it's very clear that the strongest peaks are not at the sampling instances. Thus it doesn't seem that the system is very tolerant to timing errors regardless of the roll-off factor.

It looks like to make the peaks equate to the sampling times you would need to half the transmission rate.

Am I wrong? What's going on?

An easy way to see what's going on is by plotting the signal's eye diagram.

With low roll-off factors, the eye diagram will show you that deviating a small amount from the optimum sampling time will decrease the error margin. Note that the signal peaks are not the ideal sampling points. This is the eye diagram of a sinc pulse (zero roll-off):

With larger roll-off factors, the eye remains quite open for a larger period of time around the optimum sampling time. This is why it's said that these pulses are more tolerant of timing errors. This is the eye diagram when using pulses with roll-off equal to 1:

Mind you, this doesn't mean that you can be sloppy about timing synchronization. It's just that small errors in timing will result in larger variation in the sample values of low roll-off pulses compared to large roll-off pulses.


Finding impulse response $h(t)$ by impulse matching; problem understanding Example 2.3 in "Linear systems and signals", Lathi 3rd edition Before the example it has been stated that in a system described by
$$
Q(D)y(t) = P(D)x(t), \quad (1.) \iff \\
(D^N + a_1D^{N-1} + ... + a_{N-1}D + a_N)y(t) = (b_{N-M}D^M + b_{N-M+1}D^{M-1} + ...+b_{N-1}D + b_N)x(t)
$$
if $N = M$, setting $y(t) = h(t)$ = impulse response, we will get
$$
h(t) = b_0 \delta(t) + \text{characteristic modes}
$$
If $N > M$, $b_0 = 0$ so we get
$$ h(t) = \text{characteristic modes}, \quad (2.)$$
Where the characteristic modes has been defined as the exponential terms/functions in the solution to the differential equation $Q(\lambda) = 0$. Then comes an example (Example 2.3) of a method, "impulse matching", to find the unkown coefficients of the $N$ characteristic modes in $h(t)$ in equation (2.), which I do not understand.

Find the impulse response $h(t)$ for a system specified by
$$(D^2 + 5D + 6)y(t) = (D+1)x(t), \quad (2.20)$$
In this case, $b_0$ = 0. Hence, $h(t)$ consists of only the characteristic modes. The characteristic polynomial is $\lambda^2 + 5\lambda + 6 = (\lambda + 2)(\lambda + 3)$. The roots are $-2$ and $-3$. Hence, the impulse response $h(t)$ is
$$ h(t) = (c_1e^{-2t} + c_2e^{-3t})u(t) , \quad (2.21)$$
letting $x(t)$ = $\delta(t)$ and $y(t) = h(t)$ in Eq. (2.20), we obtain
$$ \ddot h(t) +5\dot h(t) + 6h(t) = \dot \delta(t) + \delta(t), \quad (2.22)$$
Recall that initial conditions $h(0^-)$ and $\dot h(0^-)$ are bot zero. But the application of an impulse at $t=0$ creates new initial conditions at $t=0^+$. Let $h(0^+) = K_1$ and $\dot h(0^+) = K_2$. These jump discontinuities in $h(t)$ and $\dot h(t)$ at $t=0$ result in impulse terms $\dot h(t) = K_1 \delta(t)$ and $\ddot h(0) = K_1 \dot \delta(t) + K_2 \delta(t)$ on the left-hand side. Matching the coefficients of impulse terms on both sides of Eq. (2.22) yields
$$ 5K_1 + K_2 = 1, \quad K_1 = 1 \implies K_1=1, K_2=-4$$
We now use these values $h(0^+) = K_1 = 1$ and $\dot h(0^+) = K_2 = -4$ in Eq. (2.21) to find $c_1$ and $c_2$. Setting $t=0^+$ in Eq. (2.21), we obtain $c_1 + c_2 = 1$. Also setting $t=0^+$ in $h(t)$, we obtain $-2c_1 - 3c_1 = -4$. These two simultaneous equations yield $c_1 = -1$ and $c_2 = 2$. Therefore
$$ H(t) = (-e^{-2t} + 2e^{-3t})u(t) $$

I do not understand the part in boldface. Why does the jump discontinuities result in impulse terms and what is the logic behind the calculation of the impulse terms from the jump discontinuities?

From Eq. $(2.21)$ you can see that the impulse response $h(t)$ jumps from zero to some non-zero value at $t=0$. The derivative of a jump is a Dirac delta impulse, so there's a Dirac delta impulse in the derivative of the impulse response. Also the derivative of $h(t)$ jumps from zero to some possibly non-zero value at $t=0$, so there's also a Dirac delta impulse in the second derivative of $h(t)$.
A slightly different way to see this is the following. Write $h(t)$ as
$$h(t)=f(t)u(t)\tag{1}$$
with
$$f(t)=c_1e^{-2t}+c_2e^{-3t}\tag{2}$$
The derivative of $h(t)$ is then
$$h'(t)=f'(t)u(t)+f(0)\delta(t)\tag{3}$$
because $u'(t)=\delta(t)$ and $f(t)\delta(t)=f(0)\delta(t)$. The second derivative is
$$h''(t)=f''(t)u(t)+f'(0)\delta(t)+f(0)\delta'(t)\tag{4}$$
From the given difference equation we can match the coefficients of $\delta(t)$ and $\delta'(t)$ on both sides of the equation:
$$\begin{align}f'(0)+5f(0)&=1\\
f(0)&=1\end{align}\tag{5}$$
Using $(2)$ this results in
$$\begin{align}3c_1+2c_2&=1\\c_1+c_2&=1\end{align}\tag{6}$$
with the solution
$$\begin{align}c_1&=-1\\c_2&=2\end{align}\tag{7}$$
From $(1)$ and $(2)$ we obtain
$$h(t)=\left(-e^{-2t}+2e^{-3t}\right)u(t)\tag{8}$$
Also take a look at this answer to a related question.


Periodically-varying phase between transmitted sound and received sound I have a function that can find the phase difference between two signals, this function uses Hilbert transform to do so.
Here is the python code for it

def phase_shift(carr_wave, rec_wave):
     assert len(carr_wave)==len(rec_wave)
     carr_comp = hilbert(carr_wave)
     rec_comp = hilbert(rec_wave)
     c = np.inner( carr_comp, np.conjugate(rec_comp) ) / math.sqrt( np.inner(carr_comp,np.conjugate(carr_comp)) * np.inner(rec_comp,np.conjugate(rec_comp)) )
     phase_diff = np.angle(c)
     phase_diff = abs(phase_diff)
     return phase_diff

I am inserting a sample code of how this works and how many samples it need to be near accurate

f = 20000
offset = 0.25
samples_in_bit = 55
t = np.linspace(0,1,44000)
signal_1 = np.sin(2*math.pi*f*t)
signal_2 = np.sin(2*math.pi*f*t + offset)
phase_difference_1 = phase_shift(signal_1,signal_2)
phase_difference_2 = phase_shift(signal_1[:55],signal_2[:55])
print phase_difference_1,phase_difference_2

output :
0.249981325271 0.249951020276

I am doing BPSK at a frequency $= 20000\textrm{ Hz}$, duration of one bit $= 0.00125\textrm{ sec}$, with a frame rate of $44000$ per second. I am modulating 1800 bits which are all zeros.

Now I plot the phase difference of this modulated wave and actual carrier wave (Both will be same as I am sending all zeros) with perfect synchronization of both waves, the plot looks like this, which is obvious.

NOTE: I am finding phase difference at each bit (from start to end sample in one bit duration) and plotting against bit number on $x$-axis and phase on $y$-axis.

Now let's say if there is some $20$-$30$ samples of miss synchronization between carrier and modulated wave there will be a an offset phase where the plot looks like this. Which is okay that its phase difference is constant through out all bits.

After modulation, when I transmit the sound and record it (Using a good mic). I found that the phase difference plot between the transmitted wave and the recorded wave looks like this,almost a periodic function.

And this one is the same plot, but I have plotted absolute of phase differences.

To exactly know the correct sync between recorded and transmitted waves I am using cross correlation techniques which give a correct synchronization correct to $30$ frames(samples). So even if synchronization instant is wrong by $30$ frames (I am having $55$ frames in on bit duration)  I should at least get a constant phase plot just like it did on my modulated wave (before transmitting).

My question is why is the phase difference between transmitted and recorded waves constantly changing ? And that too it is not a random change. It is kind of periodic.
Does this normally happen in all signal communications? or is there any problem in my transmission?
Or can I do any filtering after I receive the signal?

EDIT :

As the problem for this is found to be frequency offset between sent and recorded waves.

Is there a method or an idea to estimate the frequency from the recorded wave, which can be used to find the phase using Hilbert function.?

I am including the frequency spectrum of my sent and recorded wave as shown below.

As Matt L said, This is due to the frequency offset between the sent (same as carrier wave as all zeros have been sent) and recorded wave. Which makes the phase difference to keep on increase with time. But as the principal values of phase lie in [-pi,pi]
that we get from Hilbert function, we see this to be a periodically varying phase.

To prove that an error in frequency between sent and recorded waves does this, i have made a test on two signals with same phase but with different frequency ( an error of 0.01% ) and plotted the phase difference using Hilbert.

Python Code snippet for the test.

By increasing this frequency error( delta_f ) the periodicity of curve increases

freq = 20000
delta_f = 2                     # frquency error of 0.01%
t = np.linspace(0,1,44000)      #framerate of 44000 to match original case
sig1 = np.sin(2*math.pi*freq*t)
sig2 = np.sin(2*math.pi*(freq+delta_f)*t + math.pi/4)
phase_arr = []
for i in range (800):
    phase = phase_shift(sig1[i*55:(i+1)*55],sig2[i*55:(i+1)*55])
    phase_arr.append(phase)
x= np.linspace(1,800,800)
plt.plot(x,phase_arr,'g')
plt.xlabel('chunk number')   #chunk size is 55 i.e number of samples in one bit (original case)
plt.ylabel('phase difference')

The output plot is as shown.


GPS | Retain Signal Tracking after Noise Bound I have an IF data, 2 seconds signal sample from only one GPS satellite, I know the PRN number of the satellite.
IF Freq = 4.092e6 Hz
Sampling Freq = 16.368e6 Hz
IF data length = 2 seconds
Total sample = 32.736e6

During the first second period, there is no noise, and SoftGNSS SDR can acquire and begin tracking the signal.
However, after 1 second, noise is bonded on the signal. The tracking function loses the signal mostly. The output has distortion a lot.
Here are the figures;

IF signal

Tracking result

The tracking ends at 1.3 seconds because I only process 1.3 seconds.
I tried to adjust DLL and PLL parameters (Damp. ratio, noise bandwidth), but no handle it. Maybe my parameters were wrong. Could you redirect me about it?

Now that the OP has clarified in the comments that this is an IF file for beginners I understand the issue. The point is to not be able to demodulate or recognize the individual chips of the GPS PRN sequence, in fact in typical application the received signal on the surface of the earth can typically be -20 dB SNR (pre-correlation) meaning it is 100x lower than the power of the noise! There is no way in such a signal to recognize any such sequence. However by multiplying with the received signal once translated completely to baseband with an aligned PRN code (the GPS C/A code for a particular satellite of interest) and then summing the result over a complete 1023 chip symbol is a great demonstration of processing gain and spread-spectrum demodulation (whose intent is to "pull signals out of the noise").
I detail this in the slides copied below, starting with the sequence as transmitted from a GPS satellite, representing the phase of the 1575.42 MHz carrier going back and forth between 0 and 180 degrees according to the sign of the C/A code being transmitted. The GPS C/A code, unique to each satellite, repeats every 1023 chips as a pseudo-random sequence as a PRN symbol at a 1.023 MHz rate. After every 20 symbols the PRN can be completely inversed (multiply by -1) according to the sign of the data bit being transmitted from the satellite, meaning the satellite will send 20 symbols to send a data bit "1" or the same 20 symbols all inverted to send a data bit "0" and thus transmits a 50 bps data rate:

Below is live data that I captured (using a Lecroy O-scope as the A/D converter hence shown in the graphic), translated to baseband (all carrier frequency and phase offsets removed) and aligned using acquisition techniques to be exactly in phase with the transmitted code which is also shown in red. The received captured data in blue for all intents and purposes appears to us a additive white Gaussian noise, yet the code sequence is part of that waveform and using correlation we can extract the symbols being transmitted from the satellite (you may think that you can recognize in parts the original PRN, but trust me, you cannot possibly!):

A histogram of the actual captured waveform further confirms the appearance of Gaussian noise (and an FFT would confirm it is "white" over the sampling bandwidth used).

Here are the samples of the received signal on a complex IQ plot after being translated back to baseband and all carrier frequency and phase offsets removed. The red dots shown are all at the chip centers, which when transmitted at the satellite existing only at +1/-1 on the real axis. This plot shows the significance of the added noise to each of those samples, totally swamping out the actual chips that were transmitted.

But after doing a sliding correlation with the same reference C/A code, we are able to clearly make out the boundaries of each of the 1023 chip symbols, and their overall phase associated with the 50 bps data transmitted by the satellite. Each sudden peak below occurs when the received waveform (after having been translated to baseband with all carrier frequency and phase removed) aligns in time with the reference C/A code, which is slid past the received waveform and after every offset a new correlation is computed (cross correlation function of the received sequence with the C/A code). We can see how there are at least 20 symbols in a row before the entire sequence is inverted (to change the data in the satellite between 1 and 0).

Hopefully this makes it all clearer that the intent of demodulation with GPS is nor to be able to extract what would appear to be a noise-free or even low noise C/A code itself, but operate at the symbol level after correlating with a received waveform where the C/A code is unrecognizably buried in the noise.  The extracted C/A code symbols after correlation become the cleaner signal (at a much lower rate) with which to them perform typical receiver processing on.  So the intent of tracking (and acquisition!) is to operate on post-correlation symbols in the latter low noise condition; the actual signals as received in practical application will never have the benefit of the high SNR as presented in the first half of the OP's waveform. Any solution that depends on such high SNR will therefore be flawed in its approach.
For further details of GPS and GPS acquisition in negative SNR condition, where some of these same figures are repeated, please see:
Auto-correlation of C/A code with 30dB noise
Phase transitions in GPS C/A signal
GPS Coarse Acquisition PRN Codes
GPS signal acquisition
GPS Receiver in Frequency Domain
Understanding time delay estimation for GPS?
Isn't Gardner's algorithm and Early-Late gate the same thing?


What method to use for interpolation and extrapolation? Below are discrete samples { t1, f(t1) }, { t2, f(t2) }, ... ,{ tn, f(tn) } using Mathematica syntax.

{{7.0,0.354887404925574},{7.3,0.4003399403324751},
{7.6,0.5849632195845474},{7.9,0.8785638270289906},
{8.2,1.235637591701261},{8.5,1.6025351028031707},
{8.8,1.925797159944803},{9.1,2.160385431197419},
{9.4,2.276553177467603},{9.7,2.2643104784151262},
{10.0,2.1348061618705962},{10.3,1.9184189440030162},
{10.6,1.6598520161702173},{10.9,1.4109826017641867},
{11.2,1.2225616685298717},{11.5,1.1360369427321668}}

Now I plot the samples.

The samples are from a bandlimited function f(t) that I defined. However, suppose they are discrete samples of an unkown function, but we know the Fourier spectrum of f(t) is band-limited to frequencies 0.25 Hz and below. What would be a good method to use for interpolating between the samples and extrapolate beyond the samples?  Here  I mentioned that polynomial interpolation has less error than cubic spline interpolation in a case like this? What method would have even less approximation error in a case like this?

@Maximilian Matthé, The discrete samples above were made using the function

$$ \mathrm{Truth}(t) = 1.4 + \frac{\sin(13-t)}{13-t} - \sin\big((1.3)(13-t)\big) $$
In this case, the sample rate is several times the Nyquist frequncy of the signal. When that is the case, we can get a very good interpolation over the samples by finding a linear combination of eight sin(f t) terms and eight cos(f t) terms that interpolate the sixteen samples. Results using eight frequencies equally spaced from ($\tfrac{\pi}{16}$ = 0.19635  radians/second) to ($\tfrac{\pi}{2}$ = 1.5708 radians/second) are shown below.

This method is far better for extrapolation (i.e. beyond the samples) than any other method I know of, and the interpolation error (i.e. between samples) is less than $4 \times 10^{-12}$.

This example was not selected to make sure the results look very good, but it is possible to find examples where the approach doesn't excel. If the same frequencies are used and the function is nearly linear over the same sample times, then cubic spline interpolation is better than an interpolating trigonometric series. However, that issue can be overcome by selecting the range of frequencies based on how much oscillation there is between the samples. If the samples are nearly linear, then a set of frequencies closer to zero could be used. Then an interpolating trigonometric series would still work very well if the range of frequencies is appropriate. Also, the samples used in this example are equally spaced, but that is not necessary for an interpolating trigonometric series to work well.


What Is Adaptive Compressive Sampling? I have just started my work in compressive sensing. the measurement vectors are obtain by multiplying the sensing matrix with input signal. the thing i cant figure out adaptive compressive sampling.

can anyone please explain it as practical ?

thank you

In "ordinary" compressive sensing you take measurements $y = Ax$ and try to reconstruct $x$ from $y$ where $A$ - of dimension $(M \times N)$ - has been generated with no knowledge of x, except that it is sparse. $A$ is typically generated randomly.

In adaptive compressed sensing, you in principle take each individual measurement $y_m,\ m = 0, 1,\ldots M$, trying to adapt you next measurement $y_m$ cleverly so you get the best possible measurement for reconstruction of $x$, given what you have already observed in $y_0, y_1,\ldots, y_{m-1}$. Effectively, this means determining the rows $A_m$ of $A$ one at a time, taking into account the measurements that previous rows $A_0, A_1,\ldots, A_{m-1}$ gave you.

One example of adaptive compressive sensing could be Malloy & Nowak, 2013.


How to emulate human sound recognition A person can recognise a large variety of sounds and fairly subtle differences in sound, even if they have only heard the sound once.  Is it possible to recreate this with software, and how?

Consider these two problems:
(1)  You have a recording of a dog bark and a cat meow (same volume).  You then record another bark, how would you go about comparing the sounds in order to tell them apart?
(2)  A person says "Yay" and "Tray".  Similar, but we can all tell the difference.

I have already investigated the following:
(1)  Fingerprinting:  Does not work unless the sounds you are comparing are identical.  So it does not work if a different dog barks.  --- Dejavu Python Package and its GitHub repo
(2)  Generic Voice recognition:  Very inaccurate if the recorded word does not come with context.  Also, dogs and cat's don't speak English.  --- Python Speech recognition package, CMU speech recognition
(3)  Phonetics/ phonemes:  The breakup of a (speech) sound into constituent pieces.  E.g: "word" becomes w ə r d.  It returns rubbish if a sound is just taken out of context. --- Python pocketsphinx library and example code

I do not think that there is a standard approach to this problem, because as far as I know it has not yet been solved satisfactorily. I believe that you can find a solution to a (very) simplified version of the problem by restricting it for example to the recognition of a fixed set of animal sounds. The way to go would probably be to use methods from speech recognition (Hidden Markov Models, HMM), because the problem is similar to voice control systems, where commands from a fixed set of predefined commands are identified. With some good engineering, these systems can work relatively well.

In the case of animal sounds you would need to define a fixed set of animal sounds that you want to be able to recognize, collect a large database of such sounds, and use the data to train HMM models, which can then be used to score unknown sounds. Of course, such a system does not even come close to human capability of distinguishing sounds, but for such a restricted application it may work well enough.

So in that sense I do not agree that voice control systems are "very inaccurate" if there is no context available, because that's the usual scenario for simple voice control systems, which can in fact be made to work in practice.


Extracting original signal from overlapping area Considering this image, where I have a spectral decomposition of a light bulb.

Due to the non point shaped light source, the spectral decomposition using a diffraction grating in front of my camera leads to overlapping spectral components. As a human though I can still perceive the structure of the main four or five spectral components because I know they should have the same shape as the lamp. Knowing the structure of the lamp, is it possible to isolate each of those spectral components with their correct amplitude, and what algorithm would be used in that case?

I'd calculate the horizontal point spread function for each color channel separately and use the spectral sensitivity ratios of the channels for frequency calibration; the resulting sensitivity-weighted spectra should obey the same ratios.

Deconvolution can be implemented as frequency domain division of each diffracted row by the corresponding direct transmission row. The division results can be averaged over all rows, excluding from statistics any noisy bins that came from division by a value too close to zero. Hopefully the rows differ enough to cover all frequency bins. Combining the color channels at the end further reduces noise.

For this to work properly, the image should be corrected for any distortion and rotated so that the point spread function due to diffraction is fully horizonal. Chromatic aberrations will be hard to fix.

A completely different approach would be kind of a matching pursuit, which should work quite fast thanks to the spiky spectra of fluorescent phosphors. There you'd try to pattern-recognize what position in the diffracted image is most similar to an intensity-scaled copy of the direct transmission image (using a color appropriate to that position, according to a calibration). Then you'd subtract the match from that position in the diffracted image and simultaneously accumulate an estimate of the spectrum. You'd repeat the process until you'd get good enough an estimate, or would no longer be able to find a good match. It should be possible to do image geometry correction in the pattern recognition part of this approach, which would also be able to avoid problems due to chromatic aberration.


How to get sensor covariance matrix from CSI values? I try to use musicdoa() function in Matlab to estimate the direction of arrival of the paths of signals received by the receiver in a wireless link. The function needs sensor covariance matrix as one of the arguments.

I have a matrix of Channel State Information (CSI) that is obtained from a WiFi card. The CSI values come in a matrix shown below:

$C=\begin{bmatrix}c_{1,1} & c_{1,2} & ... & c_{1,30}\\
 c_{2,1} & c_{2,2} & ... & c_{2,30}\\
 c_{3,1} & c_{3,2} & ... & c_{3,30}\end{bmatrix}$

In which, $c_{i,j}$ is the CSI value for $i^{th}$ antenna on $j^{th}$ subcarrier.

How can I get the sensor covariance matrix using the C matrix?

The Channel State Information (CSI) matrix is calculated over a packet of known signals that is sent as part of the communications protocol.

The CSI matrix contains two things (per carrier). The bearing of the "target" and the "strength" of transmission along the line of that bearing. So basically, this can tell you where is the transmitter with respect to the receiver array.

The CSI matrix is calculated via cross correlation between the signals received by each antenna and correlation is normalised covariance.

So, strictly speaking, it is impossible to derive covariance from the CSI because you need access to the received signals themselves. But, the angle (bearing) component of the CSI could be considered as being proportional to the real covariance for the pilot tones that are used to generate them. You could even obtain some theoretical covariance from the pilot tones by delaying them progressively but that would not include the covariance introduced due to noise.

For more information please see here and here.

Hope this helps


Efficient calculation of correlation function every $N^{\rm th}$ point I would like to calculate a long correlation function of length say, 1e6 points.
I have a prior knowledge that the correlation peak will be in point k*1000.
Is there an efficient way to apply this prior knowledge and perform a shorter calculation.

I don't want to calculate correlation by definition and implement a 1000 correlators. This solution is not sufficient.
Batch algorithm (FFT) is always possible but it does not make use of the mentioned structure.

If you have a signal of length $N$ you can calculate auto correlation function of length $K=2N-1$.You'll need perform $N$ multiplications for each value of correlation function $r_{xx}$ or in other words you need one correlator. So you need as many correlators as your $r_{xx}$ function length is. If you want to find correlation in the lag $0$ use one correlator. You can can calculate correlation of any lag $i \in (-N+1...N-1)$ via one correlator with signal properly shifted by $i$ times. From the other side $FFT$ method allows you to find $N$ taps of cyclic auto correlation function via only one $FFT$ operation. So it is not optimal if you need only a few correlation values.


Impulse response of an LTI system given the input and output signals I have been given the input and output signals of an LTI system as:

$x[n] = (\frac{1}{2})^nu[n] + 2^nu[-n-1]$

$y[n] = 6(\frac{1}{2})^nu[n] - 6(\frac{3}{4})^nu[n]$

I have found the system function $H(z)$ by using $H(z) = \frac{Y(z)}{X(z)}$ and using the standard transform tables to get:

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

$Y(z) = 6*\frac{1}{1-\frac{1}{2} {z^{-1}}} - 6*\frac{1}{1-\frac{3}{4}z^{-1}} = \frac{-\frac{3}{2}z^{-1}}{(1-\frac{1}{2}z^{-1})(1-\frac{3}{4}z^{-1})}$

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

Then this can be further simplified to:

$H(z) = 1 - \frac{\frac{5}{4}z^{-1}}{1-\frac{3}{4}z^{-1}}$

However from here I am stuck on how to get the form of the second term back into the time domain as it does not match anything from the standard tables.

I am just wondering if I have missed anything or made a mistake in my calculations, or if there is some way to get it into time domain form from here that I am not seeing, so that i can get the impule response $h[n]$.

Thank you for any help you can give me.


Finding polynomial approximations of a sine wave I want to  approximate the sine wave given by $\sin\left(\pi x\right)$ by applying a polynomial waveshaper to a simple triangle wave, generated by the function

$$T\left(x\right)=1-4\left|\tfrac{1}{2}-\operatorname{mod}(\tfrac{1}{2}x+\tfrac{1}{4},\ 1)\right|$$

where $\operatorname{mod}(x, 1)$ is the fractional part of $x$:

$$ \operatorname{mod}(x, y) \triangleq  y \cdot \left( \left\lfloor \frac{x}{y}\right\rfloor - \frac{x}{y} \right) $$

A Taylor series could be used as a waveshaper.

$$S_1\left(x\right)=\frac{\pi x}{2}-\frac{\frac{\pi x}{2}^3}{3!}+\frac{\frac{\pi x}{2}^5}{5!}-\frac{\frac{\pi x}{2}^7}{7!}$$

Given the above functions, $S_1(T(x))$ will get us a decent approximation of a sine wave. But we need to go up to the 7th power of the series to get a reasonable result, and the peaks are a little low and won't have a slope of exactly zero.

Instead of the Taylor series, we could use a polynomial waveshaper following a few rules.

Must pass through -1,-1 and +1,+1.
Slope at -1,-1 and +1,+1 must be zero.
Must be symmetrical.

A simple function meeting our requirements:

$$S_2\left(x\right)=\frac{3x}{2}-\frac{x^3}{2}$$

The graphs of $S_2(T(x))$ and $\sin\left(\pi x\right)$ are pretty close, but not as close as the Taylor series. In between the peaks and zero-crossings they visibly deviate a bit. A heavier and more accurate function meeting our requirements:

$$S_3\left(x\right)=\frac{x(x^2-5)^2}{16}$$

This is probably close enough for my purposes, but I'm left wondering if another function exists which approximates the sine wave more closely, and is computationally cheaper. I have a pretty good grasp on how to find functions meeting the three requirements above, but I'm not sure how go about finding functions which meet those requirements and also most closely match a sine wave.

What methods exist for finding polynomials which mimic a sine wave (when applied to a triangle wave)?

To clarify, I'm not necessarily looking only for odd-symmetrical polynomials, although those are the most straightforward choice.

Something like the following function could also fit my needs:

$$S_4\left(x\right)=\frac{3x}{2}+\frac{x^2}{4}+\frac{x^4}{4}$$

This meets the requirements in the negative range, and a piecewise solution could be used to apply it to the positive range as well; for example

$$\frac{3x}{2}-\frac{P\left(x,2\right)}{4}-\frac{P\left(x,4\right)}{4}$$

where $P$ is the signed power function.

I'd also be interested in solutions using the signed power function to support fractional exponents, as this gives us another "knob to twist" without adding another coefficient.

$$a_0x\ +a_1P\left(x,\ p_1\right)$$

Given the right constants, this could potentially get very good accuracy without the heaviness of fifth- or seventh-order polynomials. Here's an example meeting the requirements described here using some hand-picked constants: $a_0=1.\overline{666}, a_1=-0.\overline{666}, p_1=2.5$.

$$\frac{5x-2P\left(x,\ \frac{5}{2}\right)}{3}$$

In fact, those constants are very close to $\frac \pi 2$, and $1 - \frac \pi 2$, and $e$. Plugging those in gives something that looks extremely close to a sine wave.

$$\frac{\pi}{2}x\ +\left(1-\frac{\pi}{2}\right)P\left(x,e\right)$$

To put it another way, $x-\frac{x^e}{6}$ looks very close to $\sin(x)$ between 0,0 and π/2,1.  Any thoughts on the significance of this? Maybe a tool like Octave can help discover the "best" constants for this approach.

about a decade ago i did this for an unnamed music synthesizer company who had R&D not too far from my condo in Waltham MA.  (can't imagine who they are.)  i don't have the coefficients.

but try this:

$$\begin{align}
f(x) &\approx \sin\left(\tfrac{\pi}{2} x \right) \qquad \qquad \text{for }-1 \le x \le +1 \\
\\
     &= \tfrac{\pi}{2} x (a_0 + a_1 x^2 + a_2 x^4) \\
\end{align}$$

this guarantees that $f(-x) = -f(x)$.

To guarantee that $f'(x) \Big|_{x=\pm1} =0$ then

$$ f'(x) = \tfrac{\pi}{2}( a_0 + 3 a_1 x^2 + 5 a_2 x^4) $$

$$ a_0 + 3 a_1 + 5 a_2 = 0 \tag{1}$$

That's the first constraint.  To guarantee that $|f(\pm 1)| = 1$, then

$$ a_0 + a_1 + a_2 = \tfrac{2}{\pi} \tag{2}$$

That's the second constraint.  Eliminating $a_0$ and solving Eqs. (1) and (2) for $a_2$ in terms of $a_1$ (which is left to adjust):

$$ a_0 =  \tfrac{5}{2 \pi} - \tfrac{1}{2} a_1 $$

$$ a_2 = -\tfrac{1}{2 \pi} - \tfrac{1}{2} a_1  $$

Now you have only one coefficient, $a_1$, left to twiddle for best performance:

$$ f(x) = \tfrac{\pi}{2} x \left((\tfrac{5}{2 \pi} - \tfrac{1}{2} a_1) + a_1 x^2 -  (\tfrac{1}{2 \pi} + \tfrac{1}{2} a_1)x^4 \right) $$

This is the way I would twiddle $a_1$ for best performance for a sine wave oscillator.  I would adjust use the above and the symmetry of the sine wave about $x=1$ and place exactly one entire cycle in a buffer with a power of two number of points (say 128, i don't care) and run the FFT on that perfect cycle.

The FFT result bin 1 will be the strength of the sine and should be about $N/2$.  Now you can adjust $a_1$ to bring your 3rd harmonic distortion up and down.  I would start with $a_1 \approx \tfrac{5}{\pi} - 2$ so that $a_0 \approx 1$.  That's in bin 3 of the FFT results  But the 5th harmonic distortion (value in bin 5) will be consequential (it will go up as the 3rd harmonic goes down).  I would adjust $a_1$ so that the strength of the 5th harmonic level is equal to the 3rd harmonic level.  It will be around -70 dB from the 1st harmonic (as I recall).  That will be the nicest-sounding sine wave from a cheap, 3-coefficient, 5th-order, odd-symmetrical polynomial.

Someone else can write the MATLAB code.  How does that sound to you?


GPS CA Signal Acquisition The coarse acquisition C/A GPS signal has a symbol rate or chip rate of 1.023 Mcps with code period of 1 ms. The code bandwidth is +/-1.023MHz (to 1st sinc nulls). I'd like to record this signal using a USRP.

What's the minimum complex sample rate that I'd need to use so I could acquire this signal in software? 1 sample/chip = 1.023e6? 2 samples/chip = 2.046e6?

Since the signal is coming from a satellite, there will be an unknown carrier frequency offset due to Doppler. From what I've read, for GPS the Doppler is +/- 10 KHz. What's the maximum residual frequency offset allowable when doing the despreading of the CA code? 500 Hz? Or, what's the SNR loss due as a function of frequency missmatch due to residual Doppler?

What's the maximum residual frequency offset allowable when doing the despreading of the CA code? 500 Hz? Or, what's the SNR loss due as a function of frequency missmatch due to residual Doppler?

The maximum frequency offset is dependent on SNR required for acquisition as the roll-off of correlation versus frequency offset is due entirely to the duration of the correlation.

The magnitude of correlation vs frequency offset is a Sinc function with the first null at 1/T Hz where T is the duration of the correlation interval in seconds.  So if your correlate over 1 full PRN period, which is 1 ms, the correlation will go to zero at a 1 KHz Doppler offset.  Note the the PRN is repeated 20 times before a data transition (the data rate if 50 Hz), so you could actually correlate over 20 ms for 13 dB more processing gain (Processing gain is 10Log(N)), at the expense of reducing the Doppler tolerance (the first null would appear with a 50 Hz offset in this case!).

A good demonstration of this is the DFT itself for those familiar with it. Each bin in the DFT is a correlation over the time interval of the sequence to the particular frequency in that bin. The DFT can be described as a bank of filters with each filter approaching a Sinc magnitude response as N gets larger (when a rectangular window is used), with the first null appearing in the adjacent bin and then additional nulls in every other bin (hence we show that each bin is uncorrelated to the others). Given the time domain sequence is N samples long in time, and there are N bins also in frequency, then the spacing between bins is 1/N, as we would expect by the description above.

This is demonstrated in the generic graphic below showing the correlation to a 3 Hz sine wave (red waveform) buried in noise over a 2 second correlation interval.

From this we also see how the post correlation magnitude is reduced as a function of frequency offset, while the noise, if white, would not be affected- from which we can determine SNR vs frequency offset.

We also see the trade space involved between processing gain and Doppler offset tolerance. If the received SNR is strong, we can reduce the duration for the correlation and then benefit from a wider frequency acquisition window. For example with GPS we could correlate over half a PRN sequence (0.5 ms) at the expense of 3 dB in processing gain, with the benefit of having the first null out at 2 KHz instead of 1 KHz.  In practice I have typically used the metric of half the main lobe as my Doppler range for the correlation (so for GPS would be 1 KHz wide or +/- 500 Hz when correlating over 1 ms), but in weaker signal conditions a smaller search range may be necessary in order to maximize the SNR required for acquisition.   (Meaning how far I would step the frequency to repeat the correlation again while searching for signals; typically with a 1ms correlation, I would step 1 KHz but if I felt more SNR may be necessary I may search in finer steps, I could also increase the correlation time which would also require finer frequency steps in the search but here I was just referring to being less sensitive to the roll-off within a set correlation time by not stepping as far in frequency before repeating the correlation again).

Also see this post related to GPS acquisition and a joint delay and frequency acquisition approach that can acquire in one PRN sequence (at the expense of significant processing): GPS signal acquisition

Also to note, once within a Doppler bin, the I and Q correlated outputs can then used to determine the precise Doppler offset and do carrier tracking with implementations described here (where we see the algorithm to determine the precise frequency offset from two adjacent I Q correlation outputs is quite simple):
High modulation index PSK - carrier recovery


Is a negative Es/No for an SCPC RF signal possible, and what would such a value signify? I apologize if this isn't the correct Stack on which to ask this question, but it seemed like the best fit.  I've tried to Google an answer, but I haven't been able to find anything that applies directly to Es/No, or which is written in a way I can understand (ie, the mathematical formulas on Wikipedia are beyond me at the moment).

I have worked with legacy satellite communications equipment for years, and I am typically used to seeing receive signal performance measured as a positive Eb/No or Es/No; any carrier visible on a spectrum analyzer would have--in my experience--a positive Es/No value.

However, I was recently shown an SCPC system where the carrier's Es/No was reported as a negative value, but I can't get an explanation as to why.  Is a negative Es/No value even possible?  Specifically, is such a value possible for a carrier which appears to be above the noise floor (as seen on a spectrum analyzer)?

If it's possible, then what would such negative values signify?

I appreciate any insight anyone can offer, even if it's just to point me to another Stack.

Assuming the number is reported in dB, a negative value of $E_s/N_0$ means is that the energy per symbol $E_s$ is less than the noise spectral density $N_0$.

Even then, a negative $E_s/N_0$ seems hard to obtain if the signal's power is larger than the noise power (as measured on the spectrum analyzer). Assuming that the maximum symbol rate is used ($E_s=2B$, where $B$ is the bandwidth), and that the noise power $\sigma^2=N_0B$, then $$\frac{E_s}{N_0} = \frac{P}{2\sigma^2},$$ where $P$ is the signal power.

So, under these assumptions, a negative $E_s/N_0$ (in dB) means that the signal's power is less than twice the noise power.


Error probability for orthogonal signaling with unequal symbol probabilities I'm trying to derive the equation of error probability in orthogonal signaling. Here what I've tried:
Let $\vec{s_1} = (\sqrt{E_b} , 0)$ and $\vec{s_2} = (0, \sqrt{E_b} )$ be vector representation of $s_1$ and $s_2$. We know that $$f(\vec{y}|\vec{s_m}) = (\frac{1}{\sqrt{\pi N_0}})^2e^{-\frac{|\vec{y} - \vec{s_m}|^2}{N_0}} , \ m = 1,2$$ According to the Bayes' law, we can write $$P(\vec{s_m}|\vec{y}) = \frac{f(\vec{y}|\vec{s_m})P(\vec{s_m})}{f(\vec{y})}, m = 1,2$$Let $P(\vec{s_1}) = p_1$ and $P(\vec{s_1}) = p_2$. Using MAP rule, we want to find $m$ such that $P(\vec{s_m}|\vec{y})$ be maximized. So $\vec{s_1}$ is chosen if $$P(\vec{s_1}|\vec{y}) \gt P(\vec{s_2}|\vec{y})$$ and otherwise $\vec{s_2}$. After simplification we get $$\frac{N_0}{\sqrt{E_b}} \ln(\frac{p_1}{p_2}) \gt -2(y_1 - y_2)$$ Let $a = \frac{-N_0 \ln(\frac{p_1}{p_2})}{2\sqrt{E_b}}$ and note that $y_1 - y_2$ is a Gaussian random variable with mean $\sqrt{E_b}$ and variance $\sigma^2 = N_0$. So the probability of error is $$P_e = p_1\int_{-\infty}^a \frac{1}{\sqrt{2\pi N_0}}e^{-\frac{(z-\sqrt{E_b})^2}{2N_0}}dz + p_2\int_{a}^{+\infty} \frac{1}{\sqrt{2\pi N_0}}e^{-\frac{(z-\sqrt{E_b})^2}{2N_0}}dz$$ which simplifies to $$P_e = p_1\int_{-\infty}^{\frac{a - \sqrt{E_b}}{\sqrt{N_0}}} \frac{1}{\sqrt{2\pi}}e^{-\frac{u^2}{2}}du + p_2\int_{\frac{a - \sqrt{E_b}}{\sqrt{N_0}}}^{+\infty} \frac{1}{\sqrt{2\pi}}e^{-\frac{u^2}{2}}du = p_1Q(\frac{\sqrt{E_b} - a}{\sqrt{N_0}}) + p_2Q(\frac{a - \sqrt{E_b} }{\sqrt{N_0}})$$ Is this result correct? I couldn't find a source which includes error probability for the orthogonal signaling when prior probabilities are unequal.
Edit: Thanks to Dilip Sarwate's comment, the second exponent should be
$$-\frac{(z+\sqrt{E_b})^2}{2N_0}$$ So the final answer is $$P_e = p_1Q(\frac{\sqrt{E_b} - a}{\sqrt{N_0}}) + p_2Q(\frac{\sqrt{E_b} + a }{\sqrt{N_0}})$$ Where $a = \frac{-N_0 \ln(\frac{p_1}{p_2})}{2\sqrt{E_b}}$. We can simplify this formula and get $$P_e = p_1Q(\sqrt{\frac{E_b}{N_0}} - \frac{\sqrt{N_0}}{2\sqrt{E_b}}\ln(\frac{p_2}{p_1})) + p_2Q(\sqrt{\frac{E_b}{N_0}} +\frac{\sqrt{N_0}}{2\sqrt{E_b}}\ln(\frac{p_2}{p_1}))$$ which agrees with Dilip Sarwate's answer.

I haven't checked the details fully but the OP's revised solution is correct. Once upon a time, I had posted a to a class webpage from 1998 where I had some lecture notes giving the answer
but the link is now defunct because my University seems to have deleted that webpage from its servers.  Anyway, the answer that I came with to the OP's question is
$$P_e = \pi_0 Q\left(\mathsf{SNR} -(2\cdot\mathsf{SNR})^{-1}\ln\left(\frac{\pi_1}{\pi_0}\right)\right) + \pi_1 Q\left(\mathsf{SNR} +(2\cdot\mathsf{SNR})^{-1}\ln\left(\frac{\pi_1}{\pi_0}\right)\right) \tag{1}$$
where $\mathsf{SNR} = \dfrac{\Vert s_0-s_1\Vert}{\sqrt{2N_0}} = \sqrt{\dfrac{E_b}{N_0}}$ in this case. Note that when $\pi_1 > \pi_0$, $\ln\left(\frac{\pi_1}{\pi_0}\right) > 0$, and so $$Q\left(\mathsf{SNR} +(2\cdot\mathsf{SNR})^{-1}\ln\left(\frac{\pi_1}{\pi_0}\right)\right) < Q\left(\mathsf{SNR} -(2\cdot\mathsf{SNR})^{-1}\ln\left(\frac{\pi_1}{\pi_0}\right)\right),$$
that is, in the weighted sum of two error probabilities in $(1)$, the larger weight $\pi_1$ weights the smaller error probability, which is exactly how it should be.


Orthogonal signal generator using integer arithmetic I have a problem with implementing an orthogonal signal generator (OSG) algorithm on a microcontroller using integer arithmetic. I use this algorithm for a single-phase phase-locked loop (PLL) algorithm, for which I need an orthogonal component of a grid voltage.

The OSG algorithm is defined as follows:

$$\frac{d}{dt}v_x = \hat{\omega} \cdot \bigl((v_g-v_x)-v_y\bigr)$$
$$\frac{d}{dt}v_y = \hat{\omega} \cdot v_x$$

where $v_g$ is the measured grid voltage, $\hat{\omega}$ is the estimated grid frequency, and $v_x$ and $v_y$ are estimated components, with $v_x$ being equal to $v_g$ for ideal estimation. For this purpose, let us assume that the grid frequency is known.

The numerical integrator is implemented as follows:

$$y_k = \frac{T_s}{12} \bigr( 23u_{k-1} - 16u_{k-2} + 5u_{k-3} \bigl) + y_{k-1}$$

where $T_s=50~\mu\text{s}$ is the sample time.

Now, this algorithm works fine in floating point implementation, but is poor in integer arithmetic implementation. Here I give both implementations:

Floating point implementation

Code declaration.

float w = (2*PI)*50;
float Ts = 50e-6;
float i1u1, i1u2, i1u3, i1y1;
float i2u1, i2u2, i2u3, i2y1;

float NumInt3rd(float u1, float u2, float u3, float y1) {
    return (Ts/12)*(23*u1-16*u2+5*u3)+y1;
}

Main function.

float vg = floor(Input(0));

float vg_x = NumInt3rd(i1u1,i1u2,i1u3,i1y1);
float vg_y = NumInt3rd(i2u1,i2u2,i2u3,i2y1);

i1u3 = i1u2;
i1u2 = i1u1;
i1u1 = ((vg-vg_x)-vg_y)*w;
i1y1 = vg_x;

i2u3 = i2u2;
i2u2 = i2u1;
i2u1 = vg_x*w;
i2y1 = vg_y;

The Input(0) is a macro to get an input signal (sine wave with an amplitude of $2048$).

Integer arithmetic implementation

Code declaration.

int w = 643398L; // (2*PI)*50*2048
int i1u1, i1u2, i1u3, i1y1;
int i2u1, i2u2, i2u3, i2y1;

int NumInt3rd(int u1, int u2, int u3, int y1) {
    int iu = 23*u1-16*u2+5*u3;
    int iy = 240000L*y1;
    return (iu+iy)/240000L;
}

Main function.

int vg = (int) Input(0);

int vg_x = NumInt3rd(i1u1,i1u2,i1u3,i1y1);
int vg_y = NumInt3rd(i2u1,i2u2,i2u3,i2y1);

i1u3 = i1u2;
i1u2 = i1u1;
i1u1 = ((vg-vg_x)-vg_y)*w/2048;
i1y1 = vg_x;

i2u3 = i2u2;
i2u2 = i2u1;
i2u1 = vg_x*w/2048;
i2y1 = vg_y;

Note that I've checked for possible overflows, it never occurs. Also, the interesting thing is that the same algorithm works fine for $T_s=300~\mu\text{s}$.

I'm not that experienced with integer arithmetic implementations. Can you please give me an advice how to possibly fix this? Thanks!

Here is the estimation error for both implementations (on y-axis: percentage of the estimation error). The estimation error in case of integer artihmetic implementation is around $\pm 5\%$.

I managed to find an answer to my problem.

The problem is with rounding in integer division. For example, -4/3 will be rounded to -1, just as -5/3. Because of this, the integration error is constantly accumulated. Instead of explaining, here I give a code how to fix this.

Code declaration.

// Global variables
int w = 643398L; // (2*PI)*50*2048
int i1u1=0, i1u2=0, i1u3=0;
int i2u1=0, i2u2=0, i2u3=0;
short i1y1=0, i2y1=0;

// Numerical integrator implementation
short NumInt3rd(int u1, int u2, int u3, short y1) {
    int iu = 23*u1-16*u2+5*u3;
    int iy = 240000L*y1;
    short y0 = (((iu+iy)>>12)*2237+65536)>>17;
    return y0;
}

Main function.

// Get voltage measurements (-2048 to +2048)
short vg = (short) Input(0);

// Numerical integrators, 3rd order
short vg_x = NumInt3rd(i1u1,i1u2,i1u3,i1y1);
short vg_y = NumInt3rd(i2u1,i2u2,i2u3,i2y1);

// Downsample frequency to prevent overflow
short wb = w>>5;

// Update integrator #1 states
i1u3 = i1u2;
i1u2 = i1u1;
i1u1 = ((int)((vg-vg_x)-vg_y)*wb)>>6;
i1y1 = vg_x;

// Update integrator #2 states
i2u3 = i2u2;
i2u2 = i2u1;
i2u1 = ((int)vg_x*wb)>>6;
i2y1 = vg_y;

It should be noted that instead of using integer division, which is very expensive in terms of required number of instruction cycles, I rather use bit shift operation. For example, 1/240000 can be well approximated as 2237/2^29 - the approximation accuracy is 0.00169277%.

And here is the sine wave estimation error. As you can see, the results are much better now.


How do I know what is causing the noise in my signal I have a noisy audio file, available here (It is safe to download if you wish to)
https://ufile.io/pstrt

Now, I inspect it in matlab:

And then i plot it on normalised frequency:

Now I need to filter this signal to remove the noise. However, I have no idea what is causing the noise in the first place. Perhaps it is obvious even from the two figures I have shown but unfortunately not for me. I could play a guessing game and apply different configurations of lowpass, highpass and bandpass filters, but I have quite a few of such audio files to analyse. Thank you for any help.

Since it's home work, I can't give the full solution.

The noise seem to have been artificially added. This type of noise doesn't happen normally
The noise is basically two sine waves that really don't belong there
The best filter type to eliminate these are notch filters.


Implementation of Block LMS In the implementation of block LMS, i need one clarification.

In the 3rd step as shown in the figure attached, the summation over a product of input, $\mathbf{u}$ and error, $e$, associated with each sample value in a block, gives rise to a single value (some integer).

Every first statement (formulae) in a line belongs to Block LMS algorithm and second statement (formulae) belongs to simple LMS Algorithm. What I want is only to consider first formulas in every line.

If all the initial weights $\hat{w}(k)$ are the same (all are zeros), then this constant value will be added to that at the end of each block which results in updates weights having the same value. i.e;

  [0 0 0 0 0] + x = [x x x x x].

If this process continues till the end of input, all the filter weights will be of the same value.

Is this correct?

As applesoup says in the comments the term
$$
\mathbf{u}(kL+i)e(kL+i)
$$
is a vector, not a single value (some integer).

Why do you think it's a scalar?

To answer your question: no, it's incorrect to state that all the algorithm does is add a constant to all elements of $\hat{w}(k)$.

As per this slide from here

the block LMS error is still a scalar. It's just that the errors across multiple blocks are used to form the updates.

Another way to look at it is that $\Phi(k)$ can be written as:
$$
\mathbf{u}_k \mathbf{e}_k
$$
where $$\mathbf{u}_k =
\left[
\begin{eqnarray*}
&\mathbf{u}(kL+0)
&\mathbf{u}(kL+1)
&\mathbf{u}(kL+2)
&\ldots
&\mathbf{u}(kL+L-1)
\end{eqnarray*}
\right]$$
is an $N_u \times L$ array and
$$\mathbf{e}_k =
\left[
\begin{eqnarray*}
&e(kL+0)\\
&e(kL+1)\\
&e(kL+2)\\
&\vdots\\
&e(kL+L-1)\\
\end{eqnarray*}
\right]$$
is an $L\times 1 $ vector.

As a result $\mathbf{u}_k \mathbf{e}_k$ is an $N_u \times 1$ vector.


Determining power spectral density from a sample To start, I have absolutely no background in electrical engineering, I am a computer scientist by trade.

I'm currently working on a project where I am attempting to code a function that uses the power spectrum density of a particular signal at a given frequency.

Specifically, I want to calculate the following function

$$F(x,f) = \int [P(f) * X(x,f)] df$$

where $F$ is the function I want to calculate, $f$ is the frequency, $X$ is an  function that doesn't pertain to my question, $x$ is a variable that also doesn't pertain to the question, and $P$ is the power spectral density of a signal.

I don't have a function describing the signal, instead I have a vector of samples of the signal.  What I'm interested in calculating is $P$ using my sample vector where $P$ has been normalized such that

$$ \int [P(f)] df = 1$$

Questions:

How might I go about this in MATLAB?  I have looked at the pwelch function (Welch's power spectrum density estimate) and I think it does what I want.  However, the values seem to be scaled to $[0,\pi]$ which makes me a bit wary.  In addition, just calling the function will not scale the values so that $ \int [P(f)] df = 1$.  Is using pwelch an acceptable approach for what I'm trying to do and how can I get things scaled correctly?
In addition to Computer Science, I also have a background in math.  However, I'm a bit thrown off by the integrals that I have been provided since the result has been "evaluated" despite no upper and lower limits ($F$ should evaluate to an actual value).  Is it common notation to integrate from negative to positive infinity in signal processing using this notation?  I understand that this question is a bit vague since I haven't really explained much of what I am trying to do.  This question is meant more as a general notation question regarding convention in the signal processing / electrical engineering community.  Specifically, is there a common convention as to what limits to integrate over when the limits have been omitted?

In Matlab, you are dealing with sampled signals, so your frequency variable goes from $0$ to Nyquist. This is why the x-axis is labeled in normalized frequency that goes from $0$ to $\pi$. If you like you can use [P, f] = pwelch(___, fs) where fs is your sampling rate, and then the function will return the power spectral density (PSD) estimate P as a function of analog frequencies in the vector f. See example here.

pwelch is not the only method for estimating PSD. The algorithm choice will depend on your application. A nice summary of other competing algorithms can be found in the answer to this question.

As for these integrals, the limits should go from $-\infty$ to $\infty$ for continuous signals. But if you know that the signal is bandlimited to some range of frequencies $(f_1, f_2)$ then you can get away with integrating from $f_1$ to $f_2$. Recall the definition of PSD is inspired by continuous time signals $x(t)$ with finite energy:
$$
E = \int_{-\infty}^{\infty} |x(t)|^2 dt < \infty
$$
and denoting $X(f)$ to be the Fourier transform of $x(t)$, by Parseval's theorem we have:
$$
\int_{-\infty}^{\infty} |x(t)|^2 dt = \int_{-\infty}^{\infty} |X(f)|^2 df
$$
This leads to the intuition that $|X(f)|^2$ is the "density" of "spectral energy" around the frequency $f.$

Requiring that the PSD integrate to 1 implies that your signals of interest are unit energy. For a non-unit-energy finite energy signal, you can simply divide the signal $x(t)$ by $\sqrt{E}$ to transform it into a unit energy signal.

For implementation purposes, all your signals and transforms are discrete. So you will have to approximate these integrals with sums.


Double size inverse fourier transform This question is about 1-d signals.

I have an array of size X with the signal on which I preform Fourier transformation using FFTW.

Is there a way to use the frequency domain data from this transformation to get an array of size X * 2 that contains the time-domain signal resampled to twice the original sampling frequency?

Would that be the same as doing linear interpolation on the original array?

Copy the spectrum to an array that is long enough to contain the spectrum of the double-length array. Put all the frequencies in their right places so that the zero frequency bin, the negative frequency bins (do not exist if you used real-input-FFT), and the positive frequency bins fall in their correct places. Zero the remaining (high-frequency) bins. If your original array was of even length, depending on the flavor of FFT, you may need to split the Nyquist frequency (highest original frequency, sampling frequency / 2) bin to its identical negative and positive frequency component and divide their value by two. You may further need to scale all bins by a constant (2 in the example) to account for the way the particular FFT and IFFT implementations are normalized. Then at the end use IFFT to get back to time domain (example in Octave):

octave:1> x = [1 2 3 4]
x =

   1   2   3   4

octave:2> X = fft(x)
X =

   10 +  0i   -2 +  2i   -2 +  0i   -2 -  2i

octave:4> Y = [10 -2+2i -1 0 0 0 -1 -2-2i]*2
Y =

 Columns 1 through 6:

   20 +  0i   -4 +  4i   -2 +  0i    0 +  0i    0 +  0i    0 +  0i

 Columns 7 and 8:

   -2 +  0i   -4 -  4i

octave:5> y = ifft(Y)
y =

   1.0000   1.0858   2.0000   2.5000   3.0000   3.9142   4.0000   2.5000

The process is not equivalent to linear interpolation but to ideal sinc interpolation with wraparound of the sinc function tails to account for the assumed periodicity of the sequence.


Looking for Fastest 2D Convolution in Python on a CPU Convolutions are essential components of many algorithms in neural networks, image processing, computer vision ... but these are also a bottleneck in terms of computations...  In the python ecosystem, there are different existing solutions using numpy, scipy or tensorflow, but which is the fastest?

Just to set the problem, the convolution should operate on two 2-D matrices. We will here always consider the case which is most typical in computer vision:

a first matrix $A$ is the input and is typically large ($N \times N$ where $N$ is typically larger than $2^{10}=1024$),
a second matrix $B$ is the template and is typically smaller (say $M=128$),
the result of the convolution $C = A \ast B$ is padded such that it is of the same size as $A$.

Often, you need to do that with many images on many kernels, so a method that does it on many has a bonus point.

Thanks for any hint!

PS: the goal is at some point to sum up these results in this blog post which already contains some examples...

To provide a more quantitative approach, I have made a jupyter notebook (which can be seen as a web page here.

The results can be summarized in the following plot:

In practice, I have found that numpy was always faster by a significant amount. However, it seems new alternatives are emerging and I have to include (not exhaustive): architecture specific optimizations (IPP), simplified graphs in deep-learning algorithms (such as pyTorch), and so on. Comments welcome!


Do the DCT coefficients correspond to particular pixels? Getting the DCT coefficient based on what I've read, is usually done through a matrix, usually 8x8 or 16x16 using these formulas:

$$B_{pq}=\alpha_p\alpha_q\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}A_{mn}\cos\left(\frac{\pi(2m+1)p}{2M}\right)\cos\left(\frac{\pi(2n+1)q}{2N}\right),\quad

\begin{align}
&0\leq p\leq M-1\\
&0\leq q\leq N-1
\end{align}
$$

$$\alpha_p=\left\{
\begin{array}{ll}
1/\sqrt{M},& p=0\\
\sqrt{2/M}, & 1\leq p\leq M-1
\end{array}\right.
\quad
\alpha_q=\left\{
\begin{array}{ll}
1/\sqrt{N},& q=0\\
\sqrt{2/N}, & 1\leq q\leq N-1
\end{array}\right.
$$

However, I would like to know if it is possible to determine the corresponding DCT coefficient of a pixel or of the edges of an image? Is there a way to detect edges through DCT coefficients?

What I'm planning to do is to apply edge detection to obtain the
  locations/indices of the pixels and then get their corresponding DCT
  Coefficient

As others stated in the comment - any pixel in the given block is related to all co-efficient in the block and vice-versa.

However, understand that Edge detection is a process of finding the gradient. Linear convolution essentially does the job. What you can do is to apply equivalent of that convolution in the DCT domain which will give you an idea about the edge orientation in the DCT block domain.

While, this is relatively narrow area - there are papers which allows you to find edges in DCT blocks. See the references below:

A Compressed Domain Scheme for Classifying Block Edge Patterns by Hyun Sung Chang, and Kyeongok Kang IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 14, NO. 2, FEBRUARY 2005 145
Direct feature extraction from compressed images by Bo Shen et. al.
SPIE vol. 2670, Storage & Retrieval for Image and Video Databases IV, 1996


Low pass filter output showing instability I am fairly new to DSP, and I'm trying to implement the Pan Tompkins algorithm for QRS detection of ECG signals in MATLAB. The first stage of the algorithm consists of a second-order low-pass filter. The difference equation of the filter is

$$
y(n)=2y(n-1)-y(n-2)+x(n)-2x(n-6)+x(n-12)
$$

Here is the MATLAB code I have written.

load 108m.mat; sig_108 = val(1,:); clear val;
lowPassResult = lowPassFilter(sig_108);

figure;
subplot(2,1,1)
plot(lowPassResult);
title('Lowpass')
subplot(2,1,2)
plot(sig_108);
title('Original')

function y = lowPassFilter(x)
    %y(n)=2y(n-1)-y(n-2)+x(n)-2x(n-6)+x(n-12)
    y = zeros(1,21600);
    y(12) = 0;
    y(11) = 0;
    for n = 13:length(x)
        y(n) = 2*y(n-1)-y(n-2)+x(n)-2*x(n-6)+x(n-12);
    end
end

I am using ECG record 108 from the MIT-BIH arrhythmia database available at PhysioBank ATM. One problem I am having is that the output of the lowpass filter appears to grow without bound towards $+\infty$ as can be seen in the figure above. Where have I gone wrong in my code that could be causing this issue?

I am going to explain Robert's answer below.

Apply Z-transform to both sides of the equation. You get the following:
\begin{align}
Y(z) [1 - 2z^{-1} + z^{-2}] = X(z)[1 -2z^{-6} + z^{-12}]
\end{align}

The transfer function is given as
\begin{align}
H(z) = \frac{Y(z)}{X(z)} = \frac{1 - 2z^{-1} + z^{}-2}{1 - 2z^{-6} + z^{-12}} = \frac{(1-z^{-6})^{2}}{(1-z^{-1})^{2}}
\end{align}

The denominator represents the poles of the filter $H(z)$. In this case, there exists two poles, both of which lie on the unit circle at $z=1$. For a filter to be stable, all the poles should lie inside the unit circle, i.e., the roots of the denominator polynomial of the transfer function $H(z)$ should lie inside the unit circle.

Perhaps, you can also try a different MATLAB implementation of the above filter:

a1 = [ 1 −2 1 ];
b1 = [ 1 0 0 0 0 0 −2 0 0 0 0 0 1 ];
lp_sig = filter( b1, a1, x );


Taking the Moving Average of STFT Results I realize that this seems redundant, but I'm getting consistently better results when performing machine learning techniques that use spectral features extracted using the moving average process below, rather than using STFT results alone.  I'm a machine learning guy by trade, and I don't have a background in signal processing.  I'm hoping you guys might be able to guide me in the right direction here.

STFT process:
I use a 10 second STFT window with 9 second overlap. The result is 1 second spaced observations that contain 10 seconds of information.

STFT Moving Average process:
I use a 2.5 second STFT window with a 1/480 second increment (480Hz sample rate), then take a 10 second moving average with a 1 second increment.  The result is 1 second spaced observations that contain 12.5 seconds of information.

My reasoning for using the moving average process was that I get a smoother signal to work with, but I'm not sure that it makes sense to signal processing folks. Is there anything fundamentally wrong with doing things this way besides having to deal with the additional computational cost?

edit: It seems like the reason for the improved results is related to the theory behind using Welch's method for power estimation. Maybe?

As rb-j says in his comment, this is often how the (magnitude-squared) STFT is smoothed $-$ using a moving average window.

The only concern might be that you're using too much smoothing. A 10 second moving average window (resulting in a 12.5 second "effective" window) may mean that some artifacts take a while to show or do not show up at all.


How to intuitively understand the state space formulation of discrete time system? The SS formulation of DT system is given by

$$x[(k+1)T] = Ax(kT) + Bu(kT)$$

$$y(kT) = Cx(kT) + Du(kT)$$

Note: T is the sampling period and often omited

Can someone explain to me why the state equation $x[(k+1)T] = Ax(kT) + Bu(kT)$ is only between the currently sampling period $kT$ and the next sampling period $(k+1)T$?

How can this model be used to describe the state occurring at all times?

Let me start by answering your last question

  How can this model be used to describe the state occurring at all times?

If you know the state at one single time instant $k$ (e.g. by a given initial condition when the system starts operating), then, given an input signal for times $\ge k$, you can use the state-space equations to compute the output signal for $k+1$, $k+2$, etc. to infinity. The state at a time $k$ contains all information about the past history of the system necessary to compute its current and future output given the input signal.

As for the other question why the state equations only relate time $k+1$ to time $k$ and not also to times $k-1$, $k-2$, etc., the answer is again that the state summarizes all past history of the system. Note that in general the state $x$ is a vector. E.g. if you have a discrete-time linear time-invariant system specified by a linear difference equation with constant coefficients, the output is given by a linear combination of past output values, and current and past input values. Such a system can be implemented by a tapped delay line as shown in the figure below for a second-order system. Now if you define the state vector $x$ to be the content of all delay elements at a given time, you kind of exchange the explicit formulation of time dependency spanning several past time indices with a vector-valued equation spanning only one time interval. The more memory a system has, the higher will be the dimension of the state vector. So by using a vector difference equation (i.e. a system of equations) instead of a scalar equation, you can always obtain a first order equation relating time $k+1$ to $k$ and to no other past time indices.

"Biquad filter DF-II" by Akilaa - Own work. Licensed under GFDL via Wikimedia Commons


How many data points are required for one revolution? I have a question about an FFT. I have recorded the vibration of a gearbox with a sampling rate of 48 kHz. The speed of the input is 50 Hz, the speed of the output is 12.5 Hz. At a sampling rate of 48 kHz I have 48000 data points per second and after the output rotates 12.5 times per second I have (48000/12.5) 3840 points per rotation.

Does it make a difference if I use only 3840 points for the FFT (i.e. 1 revolution) or if I use 48000 points and thus have 12.5 revolutions?

If someone has literature on this or can explain to me whether there are rules for this, I would be very grateful.

@Hustler. Hi. I suggest you use all 48000 vibration samples. As a general rule: In mathematical analysis of measured data it's preferred that you use all available data to estimate some physical quantity. If you perform 48000-point DFTs (discrete Fourier transforms) your first DFT bin frequency will be zero Hz and your DFT bin spacing will be (as jithin said) 1 Hz. Your full spectrum analysis width will be zero –to- 24000 Hz.
When you see a narrow spectral peak in your DFT spectral magnitude samples you can relate the frequency (vibrations/second) of that spectral peak to some gear in your gear box by knowing the gear's rotational velocity and how many teeth the gear has.
If you are forced to use radix-2 FFT (fast Fourier transform) software, I suggest you zero-pad your 48000 samples out to a length of 65536 samples (2^16 samples) and then perform a 65536-point FFT. In that situation the FFT spectral bin spacing will be 48000/65536 = 0.7324 Hz and your full spectrum analysis width will again be zero –to- 24000 Hz.


Applying fourier transform twice (DSP course) I stumbled upon a question in a DSP course (coursera) which I don't understand, shown below is a screenshot of the question and answer.

The part which I don't understand is circled. Why is it equal to 0? For multiples of N I understand that each element in the summation is 1 because the exponent is 0. But why does the summation equal 0 when (i+n) is not a multiple of N?

Just use the formula for the geometric sum:
$$\sum_{k=0}^{N-1}e^{-j\frac{2\pi}{N}kl}=\frac{1-e^{-j2\pi l}}{1-e^{-j\frac{2\pi}{N}l}},\qquad l\neq mN,\;m\in\mathbb{Z}$$
Now realize that $e^{-j2\pi l}=1$ for $l\in\mathbb{Z}$ and the result follows.


Measuring Audio Signal Similarities I have recorded samples of lungs sound (person breath) signals from electronic stethoscope.
How can I measure (In percentage) the similarity between two audio signals (lung sound), the general idea is to take a clean recorded signal(without environmental noises, Laboratory conditions) and compare it with the regular recorded lung sound. If the percent is high enough then it is valid, else it is with too much environmental noises.
I thought to do it by cross-correlation but after a research I found out that it is more suitable for time lag between signals. Any good idea to start with? I want to implement it in MATLAB.
** the recorded files are in AAC format.

You can use the Normalized Cross Correlation for that.
http://en.wikipedia.org/wiki/Correlation_and_dependence

Basically, representing each recorded sound as a vector, this gives the angle between them.

Another approach is dealing with features of the sound instead of the signal itself.
You can start by reading how Shazam works:

https://www.google.com/webhp?sourceid=chrome-instant&ion=1&espv=2&ie=UTF-8#q=how%20shazam%20works

http://www.reddit.com/r/askscience/comments/1mdb08/how_do_apps_like_soundhound_and_shazam_work_is/

http://gizmodo.com/5647458/how-shazam-works-to-identify-nearly-every-song-you-throw-at-it

http://www.slate.com/articles/technology/technology/2009/10/that_tune_named.2.html

http://en.wikipedia.org/wiki/SoundHound

http://www.quora.com/How-does-SoundHound-work

http://everythingelsematterstoo.blogspot.co.il/2010/11/how-shazam-works.html


Can a digital image be considered as a tensor ? or have a tensor like structure? I am wondering if one can assume that a digital image has a tensor structure. I mean for example when we talk about isotropy in image processing, we are talking about the intensity arranged in a tensor like structure; having coordinates and values as opposed  to a matrix which carries only information about values.

Why is this important ? because I would like to explain some processes in comparison with tensor functions. For example be able to explain isotropic operations

Thank you,

An isotropic tensor is one invariant under rotations. An isotropic kernel is also invariant under rotations, but the way the rotation acts on the object is very different.

Tensors transform as linear forms, whereas images are just rotated like a scalar function of coordinates and the coordinates are rotated. So if you come from continuum mechanics, the image/kernel transforms like the pressure field.