FIR Filter Design - obtaining FIR coefficients from frequency response I would like to obtain FIR filter (low-pass) coefficients from a given frequency response by using C++ to develop algorithms. I have read at http://iowahills.com/B2PolynomialFIRFilters.html, that, if done properly, almost any filter response can be implemented as an FIR. In principle, one simply does an Inverse Fourier Transform on the filter's frequency response. Is there a C or C++ code that will do this? Any links or detailed guidelines how to solve this would be appreciated! I want to thank member hotpow2 for suggesting gradient descent routine, but I need more details on that as well.

here's some very old source in C to do an FFT.  it's nowhere near as good as FFTW, but it's free and self-contained.

if you're doing this inverse DFT and windowing (i would recommend the Kaiser window) to design an FIR, probably it doesn't need to be an FFT.  you can write a very simple DFT (no Cooley-Tukey radix-2 optimizations) and run it.  i don't think you'll care too much if it takes a couple of seconds longer to run.

/*
A set of utility programs to compute the Fast Fourier Transform (FFT):

                       N-1
                X[k] = SUM { x[n]exp(-j2(pi)nk/N) }
                       n=0

and inverse Fast Fourier Transform (iFFT):

                           N-1
                x[n] = 1/N SUM { X[k]exp(+j2(pi)nk/N) }
                           k=0

To speed things up, multiplication by 1 and j are avoided.  The input, x[],
is an array of complex numbers (pairs of doubles) of length N = 2^p.  The
calling program supplies p = log2(N) not the array length, N.  The output
is placed in BIT REVERSED order in x[].  Call bitreverse(x, p) to swap back
to normal order, if needed.  However, iFFT(X, p, trigtbl) requires its input,
X[], to be in bit reversed order making bit reversing unnecessary in some
cases, such as fast convolution.  trigtbl is an array of doubles of length
N/4 containing the sin function from 0 to pi/2 used to compute the FFT.
Call sintable(trigtbl, p) ONLY ONCE before either FFT(x, p, trigtbl) or
iFFT(X, p, trigtbl) to set up a sin table for FFT computation.

Written in Megamax C (for the Mac) by Robert Bristow-Johnson (1985).
*/

/*
#include <complex.h>
*/

typedef struct {
double real;
double imag;
} complex;

#define Re(z) (z).real
#define Im(z) (z).imag

#define PI (3.14159265359)

double sin();

FFT (x, p, trigtbl)
complex x[];
int p;
register double *trigtbl;
{
register long length, step, stepsize, size;
register complex *top, *bottom, *end;               /* top & bottom of FFT butterfly */
complex temp;

size = 1L<<p;
end = x + size;

length = size>>1;
size >>= 2;
stepsize = 1L;
while ( length > 1L) {
    top = x;
    while (top < end) {
        bottom = top + length;

        Re(temp) = Re(*top) - Re(*bottom);          /* butterfly: twiddle = 1 */
        Im(temp) = Im(*top) - Im(*bottom);
        Re(*top) += Re(*bottom);
        Im(*top) += Im(*bottom);
        *bottom = temp;
        top++;
        bottom++;

        for (step = stepsize; step < size; step += stepsize) {
            Re(temp) = Re(*top) - Re(*bottom);      /* butterfly: twiddle = exp(-jÌ) */
            Im(temp) = Im(*top) - Im(*bottom);
            Re(*top) += Re(*bottom);
            Im(*top) += Im(*bottom);
            Re(*bottom) = Re(temp)*trigtbl[size - step] + Im(temp)*trigtbl[step];
            Im(*bottom) = Im(temp)*trigtbl[size - step] - Re(temp)*trigtbl[step];
            top++;
            bottom++;
        }

        Re(temp) = Im(*top) - Im(*bottom);          /* butterfly: twiddle = -j */
        Im(temp) = Re(*bottom) - Re(*top);
        Re(*top) += Re(*bottom);
        Im(*top) += Im(*bottom);
        *bottom = temp;
        top++;
        bottom++;

        for (step = stepsize; step < size; step += stepsize) {
            Re(temp) = Im(*top) - Im(*bottom);      /* butterfly: twiddle = -j*exp(-jÌ) */
            Im(temp) = Re(*bottom) - Re(*top);
            Re(*top) += Re(*bottom);
            Im(*top) += Im(*bottom);
            Re(*bottom) = Re(temp)*trigtbl[size - step] + Im(temp)*trigtbl[step];
            Im(*bottom) = Im(temp)*trigtbl[size - step] - Re(temp)*trigtbl[step];
            top++;
            bottom++;
        }
        top = bottom;
    }
    stepsize <<= 1;
    length >>= 1;
}

top = x;
bottom = x + 1;
while (top <  end) {
    Re(temp) = Re(*top) - Re(*bottom);              /* butterfly: twiddle = 1 */
    Im(temp) = Im(*top) - Im(*bottom);
    Re(*top) += Re(*bottom);
    Im(*top) += Im(*bottom);
    *bottom = temp;
    top += 2;
    bottom += 2;
}
}

iFFT (X, p, trigtbl)
complex X[];
int p;
register double *trigtbl;
{
register long length, step, stepsize, size;
double scale;
register complex *top, *bottom, *end;               /* top & bottom of FFT butterfly */
complex temp;

size = 1L<<p;
end = X + size;

scale = 1.0/(double)size;
top = X;
bottom = X + 1;
while (top <  end) {
    Re(temp) = (Re(*top) - Re(*bottom))*scale;      /* butterfly: twiddle = 1/N */
    Im(temp) = (Im(*top) - Im(*bottom))*scale;
    Re(*top) = (Re(*top) + Re(*bottom))*scale;
    Im(*top) = (Im(*top) + Im(*bottom))*scale;
    *bottom = temp;
    top += 2;
    bottom += 2;
}

length = 1L;
size >>= 2;
stepsize = size;
while ( stepsize >= 1L) {
    length <<= 1;
    top = X;
    while (top < end) {
        bottom = top + length;

        temp = *bottom;                             /* butterfly: twiddle = 1 */
        Re(*bottom) = Re(*top) - Re(temp);
        Im(*bottom) = Im(*top) - Im(temp);
        Re(*top) += Re(temp);
        Im(*top) += Im(temp);
        top++;
        bottom++;

        for (step = stepsize; step < size; step += stepsize) {
                                                    /* butterfly: twiddle = exp(+jÌ) */
            Re(temp) = Re(*bottom)*trigtbl[size - step] - Im(*bottom)*trigtbl[step];
            Im(temp) = Im(*bottom)*trigtbl[size - step] + Re(*bottom)*trigtbl[step];
            Re(*bottom) = Re(*top) - Re(temp);
            Im(*bottom) = Im(*top) - Im(temp);
            Re(*top) += Re(temp);
            Im(*top) += Im(temp);
            top++;
            bottom++;
        }

        Re(temp) = -Im(*bottom);                    /* butterfly: twiddle = +j */
        Im(temp) = Re(*bottom);
        Re(*bottom) = Re(*top) - Re(temp);
        Im(*bottom) = Im(*top) - Im(temp);
        Re(*top) += Re(temp);
        Im(*top) += Im(temp);
        top++;
        bottom++;

        for (step = stepsize; step < size; step += stepsize) {
                                                    /* butterfly: twiddle = +j*exp(+jÌ) */
            Re(temp) = -Im(*bottom)*trigtbl[size - step] - Re(*bottom)*trigtbl[step];
            Im(temp) = Re(*bottom)*trigtbl[size - step] - Im(*bottom)*trigtbl[step];
            Re(*bottom) = Re(*top) - Re(temp);
            Im(*bottom) = Im(*top) - Im(temp);
            Re(*top) += Re(temp);
            Im(*top) += Im(temp);
            top++;
            bottom++;
        }
        top = bottom;
    }
    stepsize >>= 1;
}
}

sintable(trigtbl, p)
register double *trigtbl;
int p;
{
register long size, i;
double theta;

size = 1L<<(p-2);
theta = (PI/2.0)/(double)size;

for (i = 0; i < size; i++)
    trigtbl[i] = sin(theta*(double)i);
}

bitreverse (x, p)
register complex *x;
int p;
{
complex temp;
register long k, i, r, size, count;

size = (1L<<p) - 1L;
for (k = 1L; k < size; k++) {
    i = k;
    r = 0;
    for (count = size; count > 0; count >>= 1) {
        r <<= 1;
        r += (i & 0x00000001L);
        i >>= 1;
    }
    if (r > k) {
        temp = x[r];
        x[r] = x[k];
        x[k] = temp;
    }
}
}


Problem with simulation of the desired shape of the signal I am trying to generate a signal with specific frequency and shape. I have accomplished to get a time domain signal with the specified frequency, however the shape does not match. Is there any method to get the desired shape. Will phase be any helpful?

For a square wave as used in your example phase don't do anything.
Since the fourrier series is just based on a sine wave without phase.
The fourrier series are:
$$
S_n (x) = \frac{a_0}{2}+ \sum_{n=1}^{N} A_n \cdot \sin (\frac{2 \pi \cdot n \cdot x}{P}) \qquad \text{for integer } N \geq 1
$$

Since the square wave have only the odd frequencies. you can choose for the sinewaves: 1, 3, 5, 7, 9, 11, 13 ,...

For other functions it is different. In that case you have both sine and cosine signals. Combination of both give as result the phase.
if you have an series of 7 or 8 frequencies it is allready close to a square wave.


What Is the Difference between Parseval's Theorem and Plancherel Theorem? Wikipedia gives Parseval's theorem as follows. Parseval's Theorem

  ,     is the Fourier transform of x(t).

And Plancheral's theorem is given as

,

where  is the fourier analog.

Plancheral theorem wiki link

Both look similar to me. What exactly is the difference between the two?

According to the links you shared, and using Wikipedia's notation, Parseval assumes $A(x),B(x)$ are two square integrable (with respect to the Lebesgue measure), complex-valued functions on $\mathbb {R}$  of period $ 2\pi $  with Fourier series

$$A(x)=\sum_{n=-\infty}^\infty a_ne^{inx}$$
and
$$B(x)=\sum_{n=-\infty}^\infty b_ne^{inx}$$
Then
$$\sum_{n=-\infty}^\infty a_n\overline{b_n} = \frac{1}{2\pi} \int_{-\pi}^\pi A(x)\overline{B(x)} \, \mathrm{d}x$$
Plancherel provides a continuous version of Parsevals theorem (in a sense that $a_n$ and $b_n$ are now the $f(x)$ and $g(x)$, respectively and $A(x),B(x)$ are replaced by $ \widehat{f}(\xi), \widehat{g}(\xi)$)


Calculating covariance matrix for MVDR beamforming I am trying to calculate the covariance matrix that is required for the calculation of an MVDR beamformer. I am getting confused as to how to calculate it. I have an array of 3 microphones each with a vector of speech signals with N samples in the time domain. I apply an STFT to the signals giving a N x M matrix, where N represents the time domain or a single frame and M is the frequency bins. The covariance matrix can be described as:

$$
R_{xx} =E[xx^{H}]
$$

Thus, for a 3 microphone system we get:

$$      R_{xx} = E\left(
        \begin{matrix}
        x_1x_1^* & x_1x_2^* & x_1x_3^* \\
        x_2x_1^* & x_2x_2^* & x_2x_3^* \\
        x_2x_1^* & x_3x_2^* & x_3x_3^* \\
        \end{matrix} \right)
$$

Keeping in mind that the system will be done in real-time, would the expectation function apply a time-average across past and current signals (i.e. for the first frame divide by one, second divide by two and so on) or just current signals (i.e. divided by one)? And so, would each frequency bin in each frame have its own $R_{xx}$ giving N x M $R_{xx}$ values?

Thanks!

MVDR is a narrowband beamformer. For broadband signals it is usually applied for each frequency bin. That means that $\mathbf{R}_{xx}$ is frequency dependent. In other words, for each time you should have $M$ matrices, each one is $3\times 3$.

Now, since you usually cannot compute $\mathbf{R}_{xx}$ exactly, you perform covariance estimation $\tilde{\mathbf{R}}_{xx}$. For more information you can look here for example. In short, as you mention you perform averaging of the covariance matrix in order to improve the estimation. In order to adapt you estimation to time-varying statistics, usually exponential averaging is performed, i.e.

$$
\tilde{\mathbf{R}}_{xx}[n] = \alpha\tilde{\mathbf{R}}_{xx}[n-1] + (1-\alpha)\mathbf{x}[n]\mathbf{x}^H[n]
$$

for some $0<\alpha<1$ (usually closer to 1).


Oversampling, decimation and Nyquist I know this is a bit of a fundamental DSP question but I want to finally understand this and I think I need some help to finally put it all together.
Let's say I have an RF complex sampling ADC with a sample rate of 4GSPS.  The device has built in decimation with a minimum decimation of 4X so my maximum capture bandwidth is 1GHz (I know slightly less to allow for aliasing guard bands etc.).
So I would put a 1GHz analog anti alias bandpass filter in front of the ADC to avoid any aliasing, and capture the full 1GHz bandwidth with the ADC. At this point I effectively have 4X oversampling and my data is output to the FPGA at 1000MSPS IQ samples.
The question I have is if I now want to break up this bandwidth in the digital domain into smaller pieces for simultaneous processing of different channels, lets say 100MHz chunks for simplicity.
Is that possible with decimation without breaking Nyquist rules?  My analog anti alias filter is still 1GHz, it would seem like if I decimate by 40x to get the 100MHz bandwidth channels from the original samples that I would now have an effective sample rate of 100MSPS.
Does this work because the original data was taken with a higher sample rate and I am decimating and digital filtering non-aliased data?
I think that is the root of what i was getting at, would those filters need to be analog and limit the BW to half the final sample rate after decimation? The real question ends up being: is it possible to sample a wide-bandwidth (1GHz) and break it up digitally by digital filtering and decimating?

Let's say I have an RF sampling ADC with a sample rate of 4GSPS. The device has built in decimation with a minimum decimation of 4X (and I/Q output, per comment below) so my maximum capture bandwidth is 1GHz (I know slightly less to allow for aliasing guard bands etc.).

Ultimately the ADC is sampling at 1Gsps with I/Q output, so it's producing $2\cdot10^9$ independent samples per second.  The Nyquist rate is twice the bandwidth, so between your analog filter and the ADC's, you'd need a final bandwidth of less than 1GHz.
You're asking if you can separate the signal digitally into a number of sub-bands and then decimate.  The answer is -- yes.  There's no reason that your filtering has to be in the analog domain -- you can just as easily filter in the digital domain.
So -- make a filter bank in digital-land, and decimate the filter outputs.  Done properly, this'll work just fine.
As an example, say you separate that 1GHz signal into 8 signals, spaced 125MHz apart and filtered to only include the middle 100MHz of each signal.  Then you decimate down to 125MHz.  Now you have 16 samples (8 channels, I & Q) on a 125MHz interval.  That is the same number of independent samples per second that you started with, and you've reduced the signal content.  So -- you have less information, and just as many samples. You're fine.


Why is beginning of this FIR filtered signal attenuated? I apply an FIR filter passing the band around 50/1000 and attenuating the other three frequency components 30/1000, 70/1000 and 110/1000.

Both filters are generated using the Remez algorithm. The first filter is of length 300 and the second is of length 100.

Why is the beginning of the filtered signal attenuated and why is this attenuation effect stronger for a longer filter?

library(signal)
par(mfrow=c(3,2))

ch <-
  sin(2*pi*1:1000/floor(1000/30)) +
  sin(2*pi*1:1000/floor(1000/50)) +
  sin(2*pi*1:1000/floor(1000/70)) +
  sin(2*pi*1:1000/floor(1000/110))

# FFT of signal
barplot(abs(fft(ch)[1:120]))

# unfiltered signal
plot(ch,type="l")

filter <- function(c0,d1,d2,n) {
  fir <- remez(n=n,f=c(0,c0-d2,c0-d1,c0+d1,c0+d2,1),a=c(0,0,1,1,0,0))

  freq <- freqz(fir,n=n)
  y <- signal::filter(as.vector(fir), 1, x=ch)

  # frequency response of filter
  barplot(abs(fft(y))[1:120])
  return(list(freq = freq, fir = fir, y = y))
}

f <- filter(2*50/1000,1/80,1/40,n=300)

# first filtered signal
plot(f$y, type="l")

f <- filter(2*50/1000,1/80,1/40,n=100)

# second filtered signal
plot(f$y, type="l")

The filter you designed is a linear phase filter, therefore, response of the filter will start at $N/2$th sample, where $N $ is the filter order. As $N$ increases so does the delay of the response. You can find this delay value by computing group delay of the filter.

In case you want to reduce the delay, you can design a minimum phase filter.


Online interpolation of temporal signal I have a signal generated from an accelerometer, what I want to do is to receive the incoming signal and check if it is in a certain interval, if it is then I can apply interpolation on it(but interpolation is not the main issue). The basic problem is how to divide a signal into chunks Online and by online I mean that I don't already have an upper time limit of the signal.

A graph to describe further:

Green Signals are the actual samples I am receiving. Let suppose I have received first 2 Green Bars then I'll wait for  3rd green point and use interpolation to find a point at 10ms after that the 4th point will arrive and I want to use it to improve my interpolation by including 3rd and also 4th point until I reach 20ms where I will use the same method used on 10 ms point.

Until now what I have done is divided whole time sequence into chunks and performed interpolation but that is not very practical as in practice I won't have complete sequence already. I'm using python but it is not important. Any suggestion in this regard is highly appreciated.

If this question is in the wrong place please let me know which platform should I use.

The idea of applying interpolation works well when you have a very high sampling rate. In some sense, what you are trying to do here is fit the missing data between the second and third green bars.

Approach

This approach is motivated from the results discussed in [1].

Construct gapped IMFs (Intrinsic Mode Functions)
Fill in gaps in each mode, based on geometrical constraints
Add up all gap-filled IMFs

Though, I am not sure if this can be implemented in real-time. You can make assumption that you are processing 100 ms data in every batch.

Reference

[1] A. Moghtaderi, P. Borgnat and P. Flandrin, "Gap-filling by the empirical mode decomposition," 2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Kyoto, 2012, pp. 3821-3824.


Wrong estimation of derivatives with an extended Kalman filter I am trying to implement an extended Kalman filter (EKF) in MATLAB for the estimation of joint trajectories (angular position, angular velocity and angular acceleration) from noisy motion capture measurements using a constant angular jerk process model. This is my process model:

$x_{prior}^k = \begin{bmatrix} 1 & \Delta t & \frac{\Delta t^2}{2} & \frac{\Delta t^3}{6} \\ 0 & 1 & \Delta t & \frac{\Delta t^2}{2} \\ 0 & 0 & 1 & \Delta t \\ 0 & 0 & 0 & 1 \end{bmatrix} \cdot x_{postiror}^{k-1}$

and process noise covariance:

$Q = 2 \cdot \sigma_p^2 \cdot \begin{bmatrix} \frac{\Delta t^7}{252} & \frac{\Delta t^6}{72} & \frac{\Delta t^5}{30} & \frac{\Delta t^4}{24} \\
\frac{\Delta t^6}{72} & \frac{\Delta t^5}{20} & \frac{\Delta t^4}{8} & \frac{\Delta t^3}{6} \\
\frac{\Delta t^5}{30} & \frac{\Delta t^4}{8} & \frac{\Delta t^3}{3} & \frac{\Delta t^2}{2} \\
\frac{\Delta t^4}{24} & \frac{\Delta t^3}{6} & \frac{\Delta t^2}{2} & \Delta t \end{bmatrix}$

My estimates for the angular positions look ok, but my estiamtes for the derivatives are quite wrong. Here is a graph from a simple example application with a comparisson between the angular velocity estimated by my implementation (EKF) and a quite noisy central difference (CD).

Somehow, I cannot improve the estimation of the derivatives by tuning the covariances, initial values or (almost) anything else. So I think there must be something wrong with my implementation. But I cannot find anything, maybe I suffer from tunnel vision...

I adopted the EKF equations from the paper by Yu et al. 2.

% Apply extended Kalman filter equations
for k = 1:T

% Apply time update
if k == 1

    % Equation (15) in [Yu2004]
    X_prior(k, :) = f(x0);
    A = dfdx(x0);
    W = dfdw(x0);

    % Equation (12) in [Yu2004]
    P_prior(:, :, k) = A * P0 * A' + W * Q * W';

else

    % Equation (15) in [Yu2004]
    X_prior(k, :) = f(X_posterior((k - 1), :)');
    A = dfdx(X_posterior((k - 1), :)');
    W = dfdw(X_posterior((k - 1), :)');

    % Equation (12) in [Yu2004]
    P_prior(:, :, k) = A * P_prior(:, :, (k - 1)) * A' + W * Q * W';

end

% Apply measurement update
H = dhdx(X_prior(k, :)');
V = dhdv(X_prior(k, :)');

% Equation (11) in [Yu2004] using the pseudo inverse to handle
% badly scaled matrices
K = P_prior(:, :, k) * H' / (V * R * V' + H * P_prior(:, :, k) * H');

% Equation (14) in [Yu2004]
X_posterior(k, :) = X_prior(k, :) + (K * (Z(k, :)' - h(X_prior(k, :)')))';

% Equation (10) in [Yu2004]
P_posterior(:, :, k) = P_prior(:, :, k) - K * H * P_prior(:, :, k);

end

There seems to be one difference in the computation of the prior error covariance compared to the equations given in the paper by Welch and Bishop 3.

Yu et al.:

$P_{prior}^k = A \cdot P_{prior}^{k - 1} \cdot A^T + W \cdot Q \cdot W^T$

Welch and Bishop:

$P_{prior}^k = A \cdot P_{posterior}^{k - 1} \cdot A^T + W \cdot Q \cdot W^T$

The equation by Yu et al. produces the wrong derivatives and the equation by Welch and Bishop seems to result in instable estimations.

I compiled a simple example application with my implementation 4. It should work out of the box, if you run "TestKalmanFilter.m" in MATLAB and you should see the given graph and a small animation of a double pendulum.

Any help is greatly appreciated.
Janis

Finally, I found my mistake. It was a combination of unfavorable parameters and the wrong computation of the prior error covariance. Yu et al. and Welch and Bishop actually give the same equation. I just got confused by the indices used by Yu et al. So the correct listing for the EKF should be:

% Apply extended Kalman filter equations
for k = 1:T

% Apply time update
if k == 1

    % Equation (15) in [Yu2004]
    X_prior(k, :) = f(x0);
    A = dfdx(x0);
    W = dfdw(x0);

    % Equation (12) in [Yu2004]
    P_prior(:, :, k) = A * P0 * A' + W * Q * W';

else

    % Equation (15) in [Yu2004]
    X_prior(k, :) = f(X_posterior((k - 1), :)');
    A = dfdx(X_posterior((k - 1), :)');
    W = dfdw(X_posterior((k - 1), :)');

    % Equation (12) in [Yu2004]
    P_prior(:, :, k) = A * P_posterior(:, :, (k - 1)) * A' + W * Q * W';

end

% Apply measurement update
H = dhdx(X_prior(k, :)');
V = dhdv(X_prior(k, :)');

% Equation (11) in [Yu2004] using the pseudo inverse to handle
% badly scaled matrices
K = P_prior(:, :, k) * H' / (V * R * V' + H * P_prior(:, :, k) * H');

% Equation (14) in [Yu2004]
X_posterior(k, :) = X_prior(k, :) + (K * (Z(k, :)' - h(X_prior(k, :)')))';

% Equation (10) in [Yu2004]
P_posterior(:, :, k) = P_prior(:, :, k) - K * H * P_prior(:, :, k);

end

Thank you for your comments and votes!


Update: Fourier Transform of a shifted and scaled $\operatorname{sinc}$ signal Let $x_N$ be the function given by
$$x_N(t)=A\frac{\sin(M\pi(t-N))}{\pi(t-N)}$$
The Fourier Transform of $x_N$ is
$$\begin{align}
X_N(j\omega)&=\mathscr{F}\{x_N\}(j\omega)\\\\
&=\int_{-\infty}^\infty x_N(t)e^{-j\omega t}\,dt\\\\
&=Ae^{-jN\omega }\int_{-\infty}^\infty \frac{\sin(M\pi t)}{\pi t}e^{-j\omega t}\,dt\tag1
\end{align}$$
Enforcing the substitution $t \rightarrow t/M\pi$ reveals
$$\begin{align}
X_N(j\omega )&=\frac{Ae^{-jN\omega }}\pi\int_{-\infty}^\infty \frac{\sin(t)}{t}e^{-j(\omega/M\pi) t}\,dt\\\\
&=\frac{Ae^{-jN\omega }}2\left(\text{sgn}(-\omega/M\pi +1)-\text{sgn}(-\omega/M\pi -1)\right)\\\\
&=\begin{cases}Ae^{-jN\omega}&,|\omega|<M\pi\\\\0&,\text{elsewhere}\end{cases}
\end{align}$$
where we used the Fourier Transform of the sinc function, $\operatorname{sinc}(t)=\frac{\sin(t)}{t}$
$$\mathscr{\operatorname{rect}}(\omega)=\begin{cases}\pi&,|\omega|<1\\\\0&,\text{elsewhere}\end{cases}$$

Let $\displaystyle x(t):=\frac{\sin(2\pi(t-1))}{\pi(t-1)}$. Then we have $A=1,M=2, N=1$ so we get :
$$
X(j\omega)=\begin{cases} e^{-j\omega}&\text{if $|\omega|<2\pi$}\\
0&\text{if otherwise}\end{cases}
$$

Is this example correct? because in my book it says the answer is $e^{2\omega}$

I would hope for someone to assist me in determining if this is correct or not. I have spent entire day on this problem. I would much appreciate any help and thank you :)

HINT:
Instead of having someone else confirm your result, it would be more instructive to confirm it yourself by computing the inverse Fourier transform of the expression for $X(j\omega)$ that you obtained:
$$x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(j\omega)e^{j\omega t}d\omega\tag{1}$$
In the definition $(1)$ it is assumed that in the forward transform you use a negative exponent and no scaling constant, i.e.,
$$X(j\omega)=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\tag{2}$$
After having derived and confirmed the result yourself, it would be useful to remember the formula for the impulse response of an ideal low pass filter with cut-off $\omega_c$. Combined with the time shifting property of the Fourier transform it is straightforward to immediately write down the transform of your original function $x_N(t)$.


YIN Pitch detection Algortithm ( how do I improve my results ) I am using YIN algorithm in a school project of mine which uses pitch detection on guitar sound. I when I play a note I get random frequencies at the beginning until they stabilize. I am thinking those are probably from action of pick on the strings.

I am going through the original paper:

  Cheveigne A, Kawahara H. - YIN, a fundamental frequency estimator
  for speech and
  music

trying to reverse engineer the library and improve my results.
I am studying computer science and my knowledge of signal processing is limited.
A summary:

Step 1)
Auto-correlation :- We try to find the correlation of the signal with itself shifted by a lag within a window.

Function can possibly have infinite values. We chose the highest peak with non-zero lag. Within a lag range.(Why the highest peak? Does it mean the loudest frequency?). The paper says if upper limit for $\tau$ is high. The algorithm may chose higher order peak (what are higher order peak?)

The following steps are to improve the accuracy of the results

Step 2)
Difference Function :- Model the signal in form of a difference function.
$$ d(\tau) = \sum_{j=1}^{j=W} (x_j - x_{j+\tau})^2 $$. Which gives :
$$ d(\tau) = \sum_{j=1}^{j=W} (r_t(0) + r_{t+ \tau}(0) -2r_t(\tau)$$

So we're basically using amplitude as bias.

step 3)
Cumulative mean normalization  Replace the difference function by cumulative normalized difference to avoid selecting value with zero lag:

$$d_t(\tau')=\frac{d_t(\tau)}{(1/\tau)\sum_{j=1}^{\tau}d_t(j)}$$

step 4)
Absolute Threshold ( Could Anyone explain this section?)

step 5)
Parabolic Interpolation : Fit the $d(\tau)$ estimates to a parabolic curve.

step 6)
Choose The Best Local Estimate : Self explanatory

I am trying to compare guitar sound with a monophonic midi.

I think the parameters I should be thinking about tinkering with are window size and threshold to improve my results or I could discard first few frames.
Could anyone point me in the right direction?

The parameters I am using :

SAMPLE RATE : 44100

WINDOW SIZE : 1024

HOP SIZE : 512

THRESHOLD : 0.1

If you are not doing this in low-latency real-time, you can work backwards from the stable portion of the pitch estimate to the transient attack portion of the waveform at the beginning.

The sound of a plucked guitar string evolves in a possibly predictable pattern over time (e.g. more so than voice).  If you can estimate the onset time and/or have neighboring pitch estimates, you can adaptively set window sizes and threshold levels over time to more optimal values, as potentially determined by experimentation on some data set of guitar notes.  You can also use statistical decision theory to determine if any local pitch estimate fits the history of a reasonable spectral evolution of any guitar note, and reject outliers as noise, transients, harmonics or octave errors (harmonic and octave errors potentially being correctable errors).  This is especially useful working backwards in time, as the attack is usually noisier than the sustain/decay portion of a note's evolution.

Some improvement in steps 4) and 6) can also be acquired from psycho-acoustic experimentation with human listeners.  For instance: For an octave difference to be perceived above a certain error rate by typical human listeners, how much normalized difference in octave pitch estimate peaks over what amount of time is required?  Any difference over a smaller amount of time might be imperceptible.

ADDED: A window size of 1024 at a sample rate of 44100 only allows correlating/using a little more than 2 periods of the pitch of the lowest string of a guitar (E4 = 82.40Hz).  A 3X or 4X longer window might be more reliable for the lowest guitar notes, but shorter windows would be more responsive or provide better time locality for the higher fretted notes.


Frequency domain multiplication I am confused why zero-padding in the frequency domain does not result in linear convolution in the time-domain using the relationship:

$$corr(a,b) = ifft(fft(a) fft(b))$$

See my process below for more details:

1.

$$a = \{a_0, a_1,\dots,a_{N-1}\},\;\;\;\; w = \mathrm{rect\{N\}}$$
$$\boxed{y_c = \mathrm{conv}(a, w)}$$

2.
$$M=2N-1$$
$$X_a = \mathrm{FFT}(a, M),\;\;\;\;W_a = \mathrm{FFT}(w, M)$$
$$\boxed{y_f=\mathrm{IFFT}\left(X_a * W_a\right)}$$

3.
$$y_c=y_f$$

I am aware that zero-padding in frequency domain is interpolation in time-
domain.

4.

$$X_b = \mathrm{FFT}(a, N),\;\;\;\;W_b = \mathrm{FFT}(w, N)$$

Now zeropad the $X_b$ and $W_b$ the spectrum to size $M$
Perform frequency domain multiplication
Take inverse transform
$$\boxed{y_{zf}=\mathrm{IFFT}(X_{zb} * W_{zb})}$$

where, $X_{zb}$ and $W_{zb}$ are zero-padded spectrum of $X_b$ and $W_b$

Results: $y_{zf} \neq \left(y_c\;|\; y_f\right)$

Attached MATLAB example below

clear all;

clear;

M = 7;
N = 4;

A = 1:4;
window = ones(1, 4);

y_c = conv(A, window)

A3 = fft(A, M);
W3 = fft(window, M);

y_t = ifft(A3.*W3)        % y_t = y_c

A1 = fft(A, N);
middle_a = A1(N/2+1)/2;
A2 = [A1(1:N/2) middle_a zeros(1, M-length(A1)-1) middle_a A1(N/2+2:end)];

W1 = fft(window, N);
middle_w = W1(N/2+1)/2;
W2 = [W1(1:N/2) middle_w zeros(1, M-length(W1)-1) middle_w W1(N/2+2:end)];

y_f = ifft(A2.*W2, M)     % Expected y_f = y_t = y_c but y_f gives bad results

The result is expected as you are simply appending zeros to the product

fft(A1).*fft(A2)

Which does the time domain interpolation as you describe. However the ifft of the above product results in a series of constant values in the time domain since it is a circular convolution. Appending additional zeroes simply interpolates this for more constant values (in addition to a scaling factor).

To achieve the result you desire, which is the linear correlation of the two vectors, you must interpolate in the frequency domain prior to the multiplication (not insert zeros). This is exactly what occurs if you review the result of fft(A,M) which appends zeros in the time domain prior to taking the FFT, resulting in frequency domain interpolation.

Essentially, the process is to zero-pad the time domain samples prior to circular convolution, which results in linear convolution. Zero padding in frequency is not actually introduced in the process ifft(fft(A,M).*fft(B,M)).


Find a specific line in a picture I want to find the line between color blue and color gray in the following picture:

More exactly, the end points of that line.

Is the Hough transform a choice? If yes how should I set up the picture before using Hough?

Any way of solving this would be helpful.

In general, you want an edge detector, like the Canny edge detector, for this kind of problem. The Hough transform is useful for extracting lines rather than edges.

However, in this specific case, you're better off doing something like:

imdata = imread('grayblue.jpg');
colidx = find(diff(imdata(1,:) ~= 0);

since it's such a simple image.


Using Chebyshev polynomials in wavetables How are Chebyshev polynomials used in wavetables for audio synthesis? I dont understand how to store these inside wavetables?

How does the synthesized audio result in matching to the correct index in the wavetable and what part of the Chebyshev polynomials are actually stored in the wavetable?

Chebyshev polynomials are mainly used for waveshaping, where an input signal is distorted by applying a memory-less non-linear function to it. This shaping function is stored in a look-up table (wavetable). If you use Chebyshev polynomials, then the polynomial just defines the shaping function, and the wavetable is generated just like for any other shaping function.

If the output is given by $y=f(x)$, where $x$ is the input, then the wavetable contains values $y$ for a given range of input values $x$. The function $f(x)$ can be defined in any desired way, and one way to define it is to use Chebyshev polynomials, which for a pure sinusoidal input generate a given harmonic.


Difference between texture and illumination in image I am new in image processing.

What are texture and illumination in image?
How can we define these terms?

Very basically illumination refers to falling of light onto the objects and making them visible as they are illuminated, whereas texture is essentially geometric patterns on the surface of objects such as dots, lines, crosses, letters, symbols, figures etc. Hence they are two quite unrelated aspects of optical images, assuming the image is composed of objects and alike.

The subject of illumination is dealt within the science of photometry, in which light source characterisation, light radiation, interaction of objects with light that falls on them, and finally how humans perceive the resulting optical stimulus that comes into our eyes as reflections from our environment is studied primarily...


Energy-preserving convolution kernel? Is there a function, in continuous or (primarily) discrete time, such that convolving with it preserves the input's energy? For $x$ that is finite-valued,  finitely supported in time, and:

Real or complex?
Real?
Real and non-negative?

where "energy" refers to L2^2:
$$
\text{E}(f(x)) = \int_{a}^{b} |f(x)|^2 dx \ \ \text{or}\
 \ \sum_{n=a}^{b} |f[n]|^2
$$
and following must hold
$$
\text{E}(h \star x) = \text{E}(x)
$$
where $\star$ is linear or circular convolution (either works for me). If  there's no such kernel, can this be proven?

I'm aware of the trivial case of the delta; $h$ should be such that $x$ is transformed non-trivially  (not just a copy, conjugate, rescaling, sign flip, etc).
Idea: we can satisfy Plancherel's formulation per Envidia's answer with a distribution, example being the Dirac delta, but unsure there'd be something one could implement finitely.

I think any allpass filter should do that. Since they only change phase and not amplitude in the frequency domain, the time domain energy remains unchanged per Parseval's theorem.
$$\sum x[n]^2 = \sum |X(k)|^2$$
An all pass only changes the phase of $X(k)$ so total energy is maintained.
You could argue that an allpass is an IIR filter, so it's not suitable to be a convolution kernel (which implies FIR). But some allpass filters are FIR (delays) and you can always truncate an allpass IIR impulse response to the "desired precision".


Why do some radio stations still use amplitude instead of frequency modulation? Why is frequency modulation (FM) not used on the current amplitude modulation (AM) radio stations? The Police have FM 2 way radios that use only 5 kilocycles of modulation. Instead of the current AM radio stations that use the same bandwidth of 5 kilocycles. The limiter & discriminator or ratio detector filter out the static & loud pops that are on AM; Thus the FM would have much less static & noise than AM. FM channels would occupy the same bandwidth as AM channels. For example if WABC in New York used frequency modulation over the current frequency of 770 they would need to have the oscillator changed but the power amp would stay the same. But the receiving radios would need to be replaced with a radio to demodulate frequency modulation at 770.

Because (as you note) all the receivers would have to be replaced.
Because international agreements allot that part of the spectrum to AM.
Because an FM transmitter with 5kHz modulation needs 25kHz channel spacing to avoid problems with adjacent channels.
Because an FM transmitter with 5kHz modulation only has an audio bandwidth of 3kHz - which is far less than AM (depending on the standard used 5khz, 10kHz or 15kHz.)


Convolution of even functions is even Prove that the convolution of two even functions is an even function.

I have my own proof which I have included as an answer, but it assumes a linear time-invariant system. I want to know if there is a general way to prove this (or if it even holds for a general convolution) ?

Just do a change of variables:

\begin{align}
z(t) &= x\star y \big|_{t}\\
&=\int_{\tau=-\infty}^{\tau = \infty} x(\tau)y(t-\tau) \,\mathrm d\tau
&\scriptstyle{\text{Set}~\tau=-\lambda, \mathrm d\tau=-\mathrm d\lambda} \tag{1}\\
&=-\int_{-\lambda=-\infty}^{-\lambda = \infty} x(-\lambda)y(t+\lambda) \,\mathrm d\lambda &\scriptstyle{\text{Now remember that}~ x~\text{and}~y
~\text{are even functions}}\\
&=\int_{\lambda=-\infty}^{\lambda=\infty} x(\lambda)y((-t)-\lambda) \,\mathrm d\lambda &\scriptstyle{\text{Same as (1) except with }\lambda~\text{instead of }\tau ~\text{and with}~-t}\\
&= z(-t)\end{align}
with nary a mention of LTI systems.

For a more DSPish proof, consider that the Fourier transform of a real-valued even function of time $t$ is a real-valued even function of frequency $f$, and thus if $X(f)$ and $Y(f)$ both are real-valued even functions of $f$, then so is $Z(f) =X(f)Y(f)$ also a real-valued even function of $f$, implying in turn that $z = x\star y$ is a real-valued even function of $t$ whenever $x(t)$ and $y(t)$ are real-valued even functions of $t$.


What do I with the sampling period when I inverse $\mathcal Z$-transform? Given the input
  $$e(kT_s) = {0.4, 0.8,1.2, - 0.9}, \quad k = 0,1,2,3$$

  and transfer function
  $$T\left(z^{-1}\right)=\frac {z^{-1}-0.8z^{-2}}{1-1.1z^{-1}+0.3}$$
  Find the output $Y(kT_s),\quad k=0,1,2,3$

(I don't know why $0.3$ and $1$ haven't been added).

Starting with the inverse transform I found $$T(k)=\frac{1}{1.3}\left(\frac{1.1}{1.3}\right)^{k-1}u(k-1)-\frac{0.8}{1.3}\left(\frac{1.1}{1.3}\right)^{k-2}u(k-2)$$

I'm not sure if that is correct as I'm not so familiar with $\mathcal Z$-transform. Anyway, I don't have a sampling period in $T(k)$. Do I just put a $T_s$ next to $k$? Do I ignore $T_s$ and do the following?
$$Y(0)=T(0)E(0)$$

$T(0)=0$ so $Y(0)=0$, but why am I given $0.4$ for $E(0)$? This has to be wrong.
$$Y(1)=T(1)E(1)=\frac{1}{1.3}0.8=10.4 $$
etc etc.

There should probably be a $z^{-2}$ term with the $0.3$. It looks like a typo.

Once you've transformed to the $z$ domain (or discrete-time), and are not interested in transforming back to the continuous-time domain, there is no real need to carry the $T_s$. Just refer to $e[k], y[k]$, etc.


closed-form analytical expression for DFT coefficients I want to calculate the closed-form analytical expression for DFT coefficients of the following problem

  $$x[n]=\left\{
 \begin{array}{ll}
  1  & \mbox{if } 0\leq n \leq M-1 \\
  0 & \mbox{if }M \leq n \leq L-1
 \end{array}
\right.$$
  Write out the closed-form analityical expression for its DFT coefficients $X[k]$.

I always get the wrong solution... My attempted solution is: $$\frac{1 - e^{-j\pi k \frac{M}{L}} }{ 1 - e^{-j\pi k \frac{1}{L}} }$$

What is the correct solution and why?

The DFT of a vector of size L is L samples of the Fourier transform of the signal:

$$
X[k] = X(e^{j\theta})|_{\theta = 2\pi k/L}
$$

So now we need to calculate the FT of $x[n]$:

$$
X(e^{j\theta}) = \sum_{n=0}^{n=M-1} e^{-j n \theta}
$$

This is a geometric series which results in:

$$
X(e^{j\theta}) = \frac{1 - e^{-jM\theta}}{1 - e^{-j\theta}}
$$

Evaluate on the DFT samples to give:

$$
X[k] = \frac{1 - e^{-j M 2 \pi k/L}}{1 - e^{-j 2 \pi k / L}}
$$

You can extract half phase from top and bottom to give:

$$
X[k] = e^{-j(M-1)\pi k/L} \frac{\sin(M \pi k/L)}{\sin(\pi k/L)}
$$


Phase difference between two wave files I am doing comparison between two wave files, How to calculate phase difference between those files?
Second wave is generated by making some changes in the first one .I have tried calculating phase like this:-

y, sr = librosa.load('first.wav')
D = librosa.stft(y)
magnitude, phase1 = librosa.magphase(D)

y1, sr1 = librosa.load('second.wav')
D1 = librosa.stft(y1)
magnitude1, phase2 = librosa.magphase(D1)

then calculated the difference as (phase1-phase2)

Is this the correct way??

You need to cross-correlate those two wavfiles : from wikipedia "cross-correlation is a measure of similarity of two series as a function of the displacement of one relative to the other"

The result obtained will give you a measure of the similarity between those two audio signals depending on the time.

Therefore you will need to find the time for which the cross-correlation is at its maximum (ie : maximum similarity), it will indicates you the time delay, phase between those two signals.

A couple of advices :

1/ If your signals are periodic, you only need to compute cross-correlation over one period (or more if you want to average the noise).

2/ Don't forget to take the absolute value of the cross-correlation while searching for the maximum.


Cross-correlation, sharp peak at 0? First of all, I have to stress that I am not a professional of coding, no more than a professional of signal processing. I am a chemist that happen to be working on a project involving both. So in advance, please forgive me if my questions are not well worded, or seem trivial to you (in addition, I am not a native english speaker). I have several questions that are not directly related to each other, so I might open several threads as I got answers.

So, I am working on a project involving the interference patern of a laser, and more importantly, shifting of this pattern alongside the horizontal (x) axis. Here is the typical diffraction pattern I observe through my camera :

My goal is to precisely measure the shift of that pattern that happens during my experiment. For a better understanding of what I mean, I attach an exemple with a signal shifted from the other.Keep in mind that the shift presented here, that is of about 40 pixels, is considered HUGE, and will never really happen during my work, I am typically observing shift of about 0-5 pixels.

To achieve this goal, I wrote a C# code that computes the cross-correlation between two arrays (I reduced my image to a single 1024 pixels line, by calculating the mean of the signal alongside each 768 pixels row). The cross-correlation has been computed using this forumula (using the Accord Framework):

corr(a, b) = ifft(fft(a_and_zeros) * conj(fft(b_and_zeros))).

The code seems to work, here is the kind of cross-correlation figure I get :

This particular figure is obtained when cross-correlating two identical images (I guess I could call it auto-correlation, but I am not sure, and, as expected, I find a shift of 0 pixels). Keep in mind that this figure is the product of a modification of the cross-correlation results array, as applying Fast Fourier Transform results in the peak being "cut in half" at the edge of the array.

My problem arise when I zoom on the peak :

As you see, the peak is sharp, and I think it should be smooth. This is a problem when I want to observe low shift, as the "sharpness" of the peak hides the real maximum :

The problem disappears when the shift becomes bigger, but then the presence of a really sharp peak I can't explain is exposed :

Note that this peak doesn't come from me modifying the array to center the peak, it is present at the edge if I don't modify the results.

This ends up with my results not being realistic for low shifts (AKA the shifts I care about) as shown here :

As you can see, things are not exactly on point (this is the very start of the project, so the experimental setup is far from perfect, as is the code), but the tendancy is there, and I would really be able to solve this peak issue before going further.

So I come to you in hope of finding out why my experiment are behaving this way, and what I could do about it ! Is there a reason linked to the way cross-correlation works, or something else ? I couldn't find any issue in the code, as of now.

Thanks a lot !

Since your shift is going to be in the 0-5 pixel range, I am assuming you are going to want a method that gives you sub-sample accuracy.

I have a suggestion for you to try.  Read my article Exponential Smoothing with a Wrinkle and apply the technique to your 1D signal.  Look for some matching zero crossings in the difference case near the center of your data frame.  You should be able to determine the zero crossing by linear interpolation of the two nearest points, or do linear regression on the nearest four or five.  The linear interpolation should be sufficient.  You should be able to measure the shift at several zero crossings. If it is a true shift, all the values you calculate should be the same.  Average them for noise mitigation.

This method is easy to code and takes way fewer calculations than the method you are using.

BTW, that is a very good presentation of your question.

Ced


Intersymbol Interference Currently I have generated a ISDN 2B1Q signal and added a random signal to it,then I have designed a lowpass butterworth filter to remove the noise as much possible in MATLAB.

My next task is to add Intersymbol Interference to this 2B1Q signal and implement DFE using LMS algortihm to reduce this ISI.Can anybody help me as to how I can add ISI to this signal in MATLAB?

You can simulate an ISI (inter symbol interference) effect, by passing the message signal pulses through a dispersive communication channel model filter.

The dispersive channel will spread the pulses and make them overlap, hence creating ISI.

You can simulate a dispersive channel via nonlinear phase filtering or multi path effects.


STFT multiplied with Gaussian window vs STFT multiplied by rectangular window I am currently analysing my signal by looking at its spectrogram to determine the sinusoidal frequency content at the local sections of my signal overtime. I can do my spectrogram in 2 ways:

Directly get the FFT of each local sections. Same as I multiply each section by a rectangular window with a height 1 and get the FFT.
Multiply each section by a Gaussian window and get the FFT.

Apparently, I get 2 different pictures of my spectrogram. Number 2 is known as the Gabor Transform.

My question here is, what is the motivation for doing the Gabor Transform vs directly doing the FFT at each section of a signal?

Using a Gaussian window has one advantage over all other windows-- it yields the most compact time-frequency view of the spectrum. In other words, the product of the uncertainty in the time and frequency domains is minimized. This is due to the uncertainty principle of the Fourier transform.

Incidentally, in physics this property of the Fourier transform is the mathematical underpinning behind the Heisenberg uncertainty principle, which states that the position and momentum cannot simultaneously be known to infinite precision. Replace "position" with "time" and "momentum" with "frequency" in that statement and you can see how this applies to spectral analysis.

(There are other reasons not to use a Gaussian window as well, depending on what your application requires.)

Also, incidentally, Gabor was a physicist who first applied many mathematical tools from physics to signal processing problems, such as the analytic signal representation. It would not surprise me if his use of the Gaussian window in signal processing was carried over from quantum mechanics.


get amplitude difference via cross correlation? im currently comparing 2 signals:

A = some music
B = A + noise

and ive implemented cross correlation in frequency domain via complex multiplication to determine the time offset, which works perfectly so far.

Now what i need in addition to the time offset is the amplitude offset, means a faktor by which A should be scaled to match the amplitude of A contained in B.

Ive been given advice to do this i have to multiply my A by the cross-correlation peak but heres my problem: my peak is NOT as expected between -1 and 1 but much higher (e.g. ~300). (even though it shows the time offset perfectly)

So here are my questions:
am i missing some normalization of my CC values?
what exactly is the difference between normalized CC and "normal" CC.
what is the best way to get the amplitude factor?

Thanks!

You divide the peak of the correlation by the energy of $A$.

$B$ is a scaled version of $A$, with some noise added to it.  Thus, the equation for the cross-correlation peak is-

$ A * B' = A * (kA' + n') = k(A*A') + A*n'$

Let's ignore the noise term since we cannot predict its value.  The remaining term, $k(A*A')$, is the energy of $A$ times the scaling factor.  Thus, you simply need to divide the correlation peak by the energy of $A$ (which you know), and you have the scaling factor.  The noise term will introduce some error into the result, but there isn't really anything you can do about that other than reducing the noise or increasing the signal power in $B$.

An example- if $A$ is [1 2 1], then the energy is $1^2 + 2^2 + 1^2 = 6$.  If $B$ is twice $A$ ([2 4 2], then the cross correlation peak would occur when they were lined up, which would give us $1*2 + 2*4 + 1*2 = 12$.  Dividing the peak (12), by the energy of $A$ (6), gives a factor of 2, which is the correct result.


How are line pairs and pixels related in image resolution? I hope this is the correct place to ask this. If not, please refer me to a better place, thanks. Now to my question,

I am learning about digital image processing, and the book I am reading is discussing spatial resolution. They discuss spatial resolution in terms of line pairs. They talk about constructing a chart with pairs of dark and light lines, and how the width of a line pair is the width of a light and dark line. This is all of the information the book gives about line pairs, and I don't really understand what role these play in digital image formation and resolution.

What does a line pair actually correspond to or respresent in the context of image formation, specifically in a camera? There is a question in my book about how many line pairs/mm a camera will be able to resolve when taking an image of a subject. It gives the size of the CCD in the camera, as well as how many elements there are vertically and horizontally in the CCD. What does line pairs/mm of the CCD mean? Is this something like the number of pixels that will be resolved per CCD element? It just isn't clear to me how line pairs are realized in the camera.

I hope this question makes sense. I am a CS student so don't do much of anything with signal processing, but I would like to have a better understanding of how images are formed with respect to their spatial resolution capabilities.

The spatial resolving capability of an optical system refers to its ability to distingusih between closest (and possibly tiniest) details. The distance between a pair of lines, yields a measure of how small the distance between them can be while they are still sperated from each other. Resolution depends on several factors such as brightness level, color, and environment.

Those resolution (line) charts are designed to test optical equipment and systems, including your eyes too. Consider two black lines (with some thickness) viewd at some length and let them come closer under your control. It's known that at some small distance, you can no longer tell that there are two lines, but they will seems like a single one. That distance is the limit of your optical resolution. To get rid of the viewing distance dependence, the resolution can often be stated in a spatial angle unit.

That's also true for a camera system. It's spatial resolution is limited by two factors; first the resolution limit of the lens (its sharpness), and then the resolution limit of the CMOS / CCD sensor grid (its sampling density).


BPSK signal that has not been pulse shaped by any filter Can one consider that a BPSK signal that has not gone through any pulse shaping filter to be equivalent to a signal that has gone through a rectangular filter of which the pass BW is equal to the symbol ( bit in that case) rate?

To answer your question, NO!  A BPSK signal with a random data bit pattern that has not gone through any (additional) pulse shaping, meaning the pulse shape is rectangular in time equal to the bit duration, has a $Sinc^2$ power spectral density centered on the carrier frequency (as the Fourier Transform of a rectangular pulse is the Sinc function).

This signal has infinite bandwidth, but the power spectral density is dropping within an envelope of $1/f^2$ where f is the offset from the carrier (as expected from the $Sinc^2$ function). Any reduction of this bandwidth by filtering in the frequency domain, results in a finite transition time between symbols, so therefore induces a pulse shape in the time domain that is non-rectangular.

Another way to see this simply, a symbol that can transition from one state to the next in an infinitely short time (which is what the rectangular pulse shape is required to do) would require infinite bandwidth in order to transition that quickly. We slow down the transition rate on purpose in order to limit the bandwidth occupied since spectrum is such a valuable (and regulated) resource. We just do that very strategically to minimize inter-symbol interference, and optimize the SNR we can receive within the bandwidth used.


How Many Samples Are Identical In Linear and Circular Convolution Result? Suppose we have two sequences $x[n]$ of length $L$ , and $y[n]$ of length $P$. If we are doing a $L$ length circular convolution between the two sequences ($L>P$), how many of the resultant samples will be uncorrupted, i.e. same as linear convolution values ?

I think by intuition it should be $P-1$, but I'm not sure.

assuming $y[n]$ is padded with $L-P$ zeros (to make it the same length as $x[n]$ and your circular convolution length), then the number uncorrupted samples is $L-P+1$.

you might want to Google "Overlap Save" (OLS, a better name i have read is "Overlap-Scrap"), which is directly about this issue of uncorrupted samples in circular convolution.  the alternative is "Overlap Add" (OLA) which has no corrupted samples, but you have to overlap the tails of each block and add them.

this is well covered in books like O&S.


Processes/Transforms involved to get brainwave data from raw EEG? (Autocorrelation confusion) Not clear on what the autocorrelation function of raw EEG means physically
why can't you take the FT of a the EEG itself and get frequency data?
With BCI & basic electrode setups you can understand your brain activity:
neuron-firing --> EEG voltage reading --> brainwave readings
How I understand that this is done:
Fourier Transform: Operation bringing you from TIME domain to FREQUENCY domain

Time Domain: EEG raw ... x= time , y= volts
Frequency Domain: brainwave data... x= frequency, y=intensity

Fourier Theorem ... basically that every line/wave can be broken down into pure components (single frequency sinusoids)
First you need the autocorrelation of the raw EEG
- Correlate EEG signal with itself - compare EEG signal at one time pt -- to another time pt... shifted tau
- Get the nature of the EEG signal through time ... how voltages are shifting
FT of Autocorrelation --> PSD
We treat the signal as stationary, even though in reality it is moving in time and apply a windowing method like the Welch's Method .. and slide this window along the entire EEG signal.
Outcome is a PSD plot <-- which gives intensity of each frequency band
Thanks!

And the power spectral density is the Fourier transform of the autocorrelation.
I don't know about EEG data, but you can compute autocorrelation either in time domain or frequency domain. For instance if you need only a few samples of the autocorrelation it may be computed more efficiently using inner products than it would by $FFT^{-1}(|FFT(x)|^2)$. You know that the autocorrelation and the power spectral density are a time-frequency pair, so if the nature of information harmonic responses, oscillation modes, (power spectral density). But if the nature of the signal is chaotic with some dispersion (memory effect) you use the time domain (autocorrelation). autocorrelation are very good to analyze channel response in communication, and the signals in the brain are communication signals, so to me it makes sense to use that instead of frequency domain that could be used for identifying vibration in an engine for instance.


Multi-band audio delay filter? I'll be processing stereo music signal. What I need is to first split the frequency range into 2 bands at 50-100 Hz. In the higher band, I need to delay one of the channels by a certain number of milliseconds. The lower band has to be delayed as a whole (left+right channel) by a different value. Is this doable? Is this doable in real-time (think VST plugin for an audio player)? If it's possible, can you please give me pointers (algorithm-wise) on how to implement it?

P. S. I need to preserve the sound quality. If such a filter would ruin it, then my idea goes to the trash can.

Not sure about audio plugins but what you want to do is certaininly doable by using filters (delay line in your case in the subbands. Split the band using using a Perfect Reconstruction Analysis filter bank. Now plugin the delay line and then combine two subbands using synthesis filter bank. Note that you need to be extra careful while plugging in the filter between analysis and synthesis filter bank. You can check some of these papers for ideas on how to do subband filtering.

https://dl.acm.org/citation.cfm?id=2817121

http://ieeexplore.ieee.org/abstract/document/1661283/


noise of zero mean and variance? It might be a very silly question, however I'm facing problem with it. Now, the Gaussian noise is a signal with zero mean and variance. I understand what the mean and the variance mean however, I'm facing problem with zero mean. What does that mean? Does it mean that it is only a signal with variance which represents the magnitude of the signal? And if this is the case, why we bother ourself with variance since the magnitude of this signal is random?  I want to understand that because I want build a sensor (a function in C/C++) and add some noise to the data. In MatLab, there is randn to generate Gaussian noise, but I have no idea how to do same concept in C++.

Edit: Let's say I have a sensor with $\pm$ 0.2 as an error. What I understand that the noise should be a Gaussian noise with zero mean and variance. Now, I have a real data which represents the position in 1D. My sensor now should give the position with some uncertainty. This is my code in C++11

#include <iostream>
#include <chrono>
#include <random>

int main(int argc, const char * argv[])
{

    double noise; // Gaussian noise with zero mean and variance 1
    double data_sensor;
    double position(2.0);

    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
    std::default_random_engine engine(seed);
    std::normal_distribution<double> dist(0, 1);

    for (int i(0); i < 10; ++i)
    {
        noise = dist(engine);
        data_sensor = position + 0.2*noise;
        std::cout << "Real Data: " << position << "  Sensor Data: " << data_sensor << std::endl;

    }
    return 0;
}

The output

Real Data: 2  Sensor Data: 2.19025
Real Data: 2  Sensor Data: 1.96754
Real Data: 2  Sensor Data: 2.03551
Real Data: 2  Sensor Data: 1.58706   << Why this out of the range
Real Data: 2  Sensor Data: 2.11237
Real Data: 2  Sensor Data: 1.98482
Real Data: 2  Sensor Data: 1.94809
Real Data: 2  Sensor Data: 1.8264
Real Data: 2  Sensor Data: 1.9776
Real Data: 2  Sensor Data: 1.70071   << Why this out of the range

The Box–Muller method uses two independent random numbers distributed uniformly on $(0,1)$, then uses such two uniformly distributed numbers to produce values with normal distribution.

Zero-mean indicates the average value of the detected signal at each time point is zero. In digital signal processing, zero often refers to $2^{n-1}$ where $n$ is the bit number of digitizer. For example, if you have a 16-bit signal ($65536$), then $32768$ represent the zero in double, and those numbers that less than $32768$ are negative.


Simulating Range Bearing Sensor with MATLAB with Gaussian Noise (Generating Gaussian Colored Random Vector) I would like to simulate a sensor that provides range and direction of a beacon. This is for EKF localization, so the noise must be Gaussian (i.e. $\mathcal{N}(0, \sigma^{2})$. Also, I would like to compute the measurement noise matrix $Q_{k}$ which is
$$
Q_{k} = \begin{bmatrix} \sigma^{2}_{r} & 0 \\ 0 & \sigma^{2}_{\theta} \end{bmatrix}
$$

This is what I've done

sigma_range = 1;
sigma_angle = degtorad(5);
Q = [sigma_range^2,               0;
                 0, sigma_angle^2];

In the sensor, this is I've done

% add Gaussian noise with zero mean and some variance
Z(1) = Z(1) + sigma_range*randn(); % Z(1) the distance btw robot and beacon
Z(2) = Z(2) + sigma_angle*randn(); % Z(2) the angle btw robot and beacon

Is this correct?

Few notes:

First you should multiply the noise by the Standard Deviation (Root of the Variance for zero mean noise).
You can do that by multiplying the Lower Cholesky Decomposition of matrix by a column vector of Gaussian noise. Yet since you assume no correlation, you can do that by independent multiplication.

Something like:

vZ = vZ + [sigma_range; sigma_angle] .* randn([2, 1]);


Scipy - Audio Processing I am looking for good tools for Audio signal processing. e.g Speech & music analysis, automatic language identification, etc.

Does the Scipy library provide functions for audio processing? Does it have a good tool for Audio Signal processing?

Can you please suggest a tool for this?

Personally I find Python one of the best choices out there and did myself some work in area of audio identification. You are welcomed to check for instance my software for automatic identification of birds from noisy audio recordings: Ornithokrites. The program is used by Department of Conservation of New Zealand and they are happy about it. Based on this example I would like to point out several advantages of using Python:

Huge, fast developing community providing tons of libraries. SciPy provides plethora of methods for signal processing (granted, not that many and mature as Matlab). Mind though that SciPy, although one of the most important, is only one of hundreds that can help you in your endeavours. I found Aubio best for music analysis. For speech and music recognition for sure you will enjoy great number of audio features Yaafe can extract.
It's free! Once out of academia, you quickly find out that Matlab is rather expensive. And even if you can afford it, then your perspective users will not be happy about this dependency. For instance, mentioned Department of Conservation would not accept proprietary software.
Identification often requires machine learning and Python has great toolkit for it: sklearn. It is state of the art library - and easy to use. Have a look at Kaggle competitions  (machine learning) and check how many top programmers are using Python and sklearn.
You can manage "big data". If you want to run analysis against huge networked databases of recordings, then Python has well established set of tools. I don't think Matlab / Octave  interface easily with e.g. Hadoop, although please do correct me if I am wrong. R does better this area.
Speaking of interfacing, you can easily interface your program with a web site. This is the way I manage Ornithokrites (bird recognition): the program runs on Amazon Web Services cloud computing service. Great if you want to provide your software to other people who do not necessarily want to go through installation procedure of all required libraries.

My second choice would be R. Although not that feature-rich as Python, it has great number of useful libraries (check e.g. seewave for your applications). Installation of those on both Windows and Linux is piece of cake, which is important if you would like others to use your program. However, to my experience high-performance computing in R is more difficult - an important thing to notice if you need to do A LOT of processing and identification.

Examples of music classification in Python:

Lecture
Audio fingerprinting
Jack Minardi blog

Book Building Machine Learning Systems with Python has a chapter on music classification

Other tools (list by no means complete): Python in Music


Why a pair of complex-conjugate zeros provide a nulling filter? (FIR filter case) FIR filter having a pair of complex-conjugate zeros that lie on the unit circle, with zeros of the form:
$$
z_1 = e^{j\omega_i}\qquad\text{and}\qquad z_2 = e^{-j\omega_i}
$$
And transfer function:
$$
H(z) = \left(1 - z_1z^{-1}\right)\left(1 - z_2z^{-1}\right) = 1 - 2\cos\left(\hat{\omega}_i\right)z^{-1} + z^{-2}
$$

Has the following pole-zero map, and magnitude & frequency response graph:

But what is the intuitive idea behind that? And how the locations of zeros will affect the result? Thank you!

The unit circle on the z-plane represents the frequency axis, similar to the imaginary axis $j\Omega$ on the s-plane for the Laplace Transform in the continuous time case. So the frequency response of the system is given by $H(z)$ when $z= e^{j\omega}$ with $\omega$ going from $0$ to $2\pi$ representing the normalized fractional radian frequency (which is the continuous time radian frequency $2\pi f$ divided by the sampling rate $f_s$.

That said, any zero on the unit circle will create a null in the frequency response. With the OP's case of complex conjugate zeros (resulting in a real response), two nulls would result as shown. The location of the zero if on the unit circle, is the fractional radian frequency where $H(z) = 0$, thus called a "zero".

If the zero is not on the unit circle, the null will not be zero but will be lower the closer the zero is to the unit circle for that frequency.

This may be clearer from the plot below showing the frequency response for a 2 point moving average filter, which has a zero at $z= -1$. The frequency response is $H(z)$ as z sweeps over the unit circle, thus giving a numerator magnitude as the difference between z at any point on the unit circle and the zero location: $z-q_z$ (or the multiplication of multiple such magnitudes if there is more than one zero), and a denominator magnitude given by the same for the pole locations: $z-q_p$. In this case the pole is at the origin, so $z-q_p=1$ for all $z=e^{j\omega}$. What also should be now clear is how the resulting phase response is formed since the net phase will be the difference between the phase of the numerator and the phase of the denominator (phases subtract in the division of complex numbers).

This type of nulling filter (zero-only) is not very effective given the gradual roll-off in frequency. To achieve very sharp nulls, place a pole very close to the zero; the closer the pole, the sharper the response! Given all poles must be inside the unit circle for a stable causal linear time invariant system, the magnitude of the pole would therefore be less than but close to 1.

This IIR approach is further detailed here: Transfer function of second order notch filter

Also this is an excellent write-up by Richard Lyons on linear phase nulling (or notch) filters that do provide a sharp notch with an FIR appraoch. This could similarly be translated to provide a notch at any frequency: https://www.dsprelated.com/showarticle/58.php


FIR Linear Phase Condition This is the linear phase condition for FIR filters as expressed by my prof:

I don't understand why $G(f)$ can be negative. Isn't the Fourier transform expressed in polar form ? So the magnitude is always positive.

You can describe a frequency response in terms of its real-valued amplitude and its phase:

$$H(f)=A(f)e^{j\phi(f)}\tag{1}$$

Note that $A(f)$ is not the magnitude, but a bipolar amplitude function. Equivalently, you can express $H(f)$ in terms of its magnitude and its phase:

$$H(f)=M(f)e^{j\tilde{\phi}(f)}\tag{2}$$

Now we have $M(f)=|A(f)|\ge 0$, and, consequently, the phase $\tilde{\phi}(f)$ in $(2)$ is generally different from the phase $\phi(f)$ in $(1)$. For frequencies for which $A(f)<0$ is satisfied, we have to add $\pi$ (or $-\pi$) to the phase to compensate for the sign change:

$$\tilde{\phi}(f)=\begin{cases}\phi(f),& A(f)\ge 0\\\phi(f)\pm\pi,&A(f)<0\end{cases}\tag{3}$$

For a linear phase filter of type I or type II (even symmetry), the phase $\phi(f)$ is a linear function of $f$, where $\phi(f)$ is defined as in $(1)$, i.e., with a bipolar amplitude function. The phase $\tilde{\phi}(f)$ as defined in $(2)$ has jumps wherever $A(f)$ changes its sign.

You can extend the definition to type III and type IV linear phase filters (odd symmetry) by expressing $H(f)$ as

$$H(f)=jA(f)e^{j\phi(f)}\tag{4}$$

where $A(f)$ is again a real-valued bipolar amplitude function, and $\phi(f)$ is a linear function of $f$.

Note that by linear function I mean a strictly linear function without an additive term, i.e., $\phi(f)=a
\cdot f$, with $a=-2\pi\tau_g$, where $\tau_g$ is the filter's group delay.


What word can I use to describe time-varying spectral content? I have two devices that generate different classes of signals, one simple and one complex. I would like to find a single word or short phrase that describes the increased richness in one class compared to the other. I'm wondering if such a word or phrase exists.

Class #1 is limited to pure tones and colored noise. These sounds are time-invariant, and they are described as having attributes of "frequency and duration," where frequency is taken to mean any number of frequencies together, so long as they are invariant over time.

Class #2 may be any kind of sound. Typically the sounds used are more complex than those in class #1, with time-varying features. They usually have phase and amplitude changes affecting various frequencies over time. Often they contain semantic content such as recognizable real-world sounds or spoken words. This class of sound has been loosely described as having an attribute of "spectral content" as opposed to an attributes of "frequency and duration," but that description doesn't seem to be very good because people may consider spectral content to mean a distribution of time-invariant frequencies.

Perhaps "time-varying spectral content" gets closer? That might be a problem, because the rich class is a superset of the simple class, so there's no requirement for a signal to vary over time. In any event, I'm not sure having time-varying spectral content rather than a frequency and duration really captures the qualitative difference between white noise and a song fragment. Spectro-temporal variations seem to be part of it, and the temporal envelope is certainly important .

Rich sounds would be things like a child laughing, a dog barking, or a human speaking your name, while simple sounds would be things like an 800 Hz sine wave, or white or pink noise.

Is there a pre-existing word or short phrase that meaningfully describes the kind of sonic richness that distinguishes one from the other?

The notions of harmonicity or non-stationary seem connoted.

Stochastic White noise or pink noise, or deterministic pure tones, are (strongly) structured in the spectral domain. So,  strongly structured, or strongly  spectrally-structured. So, in a soft manner, I would  describe the rich class as weakly structured, or weakly spectrally-structured.


taylor expansion of scale space function I see the following expression from http://en.wikipedia.org/wiki/Scale-invariant_feature_transform

The quadratic Taylor expansion of the Difference-of-Gaussian scale-space function,  with the candidate keypoint as the origin is

$$D(\textbf{x}) = D + \frac{\partial D^T}{\partial \textbf{x}}\textbf{x} + \frac{1}{2}\textbf{x}^T \frac{\partial^2 D}{\partial \textbf{x}^2} \textbf{x}$$

where D and its derivatives are evaluated at the candidate keypoint and $\textbf x = (x,y,\sigma)$  is the offset from this point.

I'm very confused about this expression. I don't know why $D^T$ appears and don't understand everything. Is there someone to help me?

The quantity $\frac{\partial D}{\partial \textbf{x}}$ is a vector, since it is the derivative of the scalar function $D(\textbf{x})$ w.r.t. all the elements of $\textbf{x}$. In the formula it is assumed that all vectors are column vectors, so in order to compute the dot product of the derivative $\frac{\partial D}{\partial \textbf{x}}$ and the vector $\textbf{x}$, you need to transpose one of them, which gives you a matrix product row times column (= scalar). The expression $\frac{\partial D^T}{\partial \textbf{x}}$ is simply the transpose of the vector of derivatives. Maybe it would have been clearer to write it as $\left(\frac{\partial D}{\partial \textbf{x}}\right)^T$.

Similarly, for the last term you have a row vector times a matrix (second derivative of $D$ w.r.t. $\textbf{x}$) times a column vector, which again results in a scalar value.


how to design LPH and HPF from APF I came across a concept that says we can design LPH and HPF from APF when the transfer function of the two APFs given as

$$\begin{align}A_0(z)&=1\\
A_1(z)&=\frac{-a+z^{-1}}{1-az^{-1}}\\\\
H_{LP}(z)&=A_0(z)+A_1(z)\\
H_{HP}(z)&=A_0(z)-A_1(z)\end{align}$$

i tried to implement this in matlab but i did not get the desired frequency response..kindly help me with the matlab code for the same..

Since you don't say what response you're actually getting from Matlab, this is all guesswork.

So let's do the basics: A little formula juggling.

$$\begin{align}A_0(z)&=1\\
A_1(z)&=\frac{-a+z^{-1}}{1-az^{-1}}\\[4em]
H_{LP}(z)&=A_0(z)+A_1(z)\\
   &= 1 + \frac{-a+z^{-1}}{1-az^{-1}}\\
   &= \frac{{1-az^{-1}}-a+z^{-1}}{1-az^{-1}}\\
   &= \frac{(1-a)\quad +\quad (1-a)z^{-1}}{1-az^{-1}}\\[4em]
H_{HP}(z)&=A_0(z)-A_1(z)\\
   &= 1 - \frac{-a+z^{-1}}{1-az^{-1}}\\
   &= \frac{{1-az^{-1}}+a-z^{-1}}{1-az^{-1}}\\
   &= \frac{(1+a)\quad +\quad (-1-a)z^{-1}}{1-az^{-1}}
\end{align}$$

So, taken from these formulas, the filter coefficients for $H_{LP}$ are $b=[1-a,1-a]$, $a=[1,-a]$, and  for $H_{HP}$ it's $b=[1+a,-1-a]$, $a=[1,-a]$ for the canonical form of recursive filters

$$H= \frac{\sum\limits_{i=0}^N b_i z^{-i}}{\sum\limits_{i=0}^M a_i z^{-i}}\text.$$

Plugging that into scipy's scipy.signal.freqz for different $a$ yielded (abscissa: normalized frequency $\in[0,\pi[$, ordinate: magnitude):

$H_{LP,a}$ for $a\in\{0,0.1,\dots,0.9,1\}$:

$H_{HP,a}$ for $ a\in\{0,0.1,\dots,0.9,1\}$:

which both very nicely match the terms low and high pass filters, respectively.


conversion of 1d signal to 2d signal for CNN operation I have a 1D signal data sampled at 25kHz and I want to modify in to 2D signal so as to modify it as the input to ConvolutionalNets as they work only on 2 dimensional data.
What are possible ways to convert the signal from 1D to 2D for such calculations. I have seen some papers but I couldnt find any relevant understanding how the conversion of signal is modified as to apply to CNN operations.
For example if I have 25 data samples in the input and if i rearrange in to a 2d signal as 5X5 sized(inbuilt function in matlab: "reshape").If I am applying 3x3 filter for convolution operation, smallest filter and the output of that 5x5 size input will be 3x3 size.
Can I apply such input to any convnets? Is that a valid type of conversion to apply for CNN ? In such kind of inputs I understand there will be relation in row wise but not column wise.
Any explanation or suggestion would be helpful.

Updated for the comment

I want to convert 1d signal in to 2d image signal. In an image a pixel will have relation with the adjacent 8 pixels. If I am rearranging a 1D signal to 2D signal as quoted in the example , the row 1 will have 5 data samples and next row starts with 6th sample and so on. If I take a sample in row 2, 2nd sample, I have relation first and third sample in the same row but nothing with samples from row 1 and row 3. So If Im converting such way is that appropriate ? In the paper mentioned ( Fig 1) I wanted to know how the 1D signal is converted and applied to CNN link

There are a couple things with that question.

1) Why do you need 2D convolutions exactly? You could think of a 1D signal as an $Nx1$ 2d signal if you like ;) Then you can design your convolutional filters to be $7x1$, $5x1$, $3x1$ etc. The pooling you only apply in the x-direction. And there you go. It's a 2D convolutional neural network, which works on 1D signal. Just kidding, it just a CNN.

2) Simple fully-connected deep learning neural networks are better for handling 1D data, as far as I know. Of course it's best if your data is evenly sampled. CNN comes handy when the fully connected structure explodes the computation. In your case, that should not be a big issue.

If your 1d data vector is too large, just try subsampling, instead of a convolutional architecture.

3) Converting a 1d data to 2d is probably valid only if you know in advance that this 1d manifold carries non-uniform neighborhood information, which could be represented with a 2D matrix with nearby connections. The paper probably inserts a layer to the network which re-arranges the indices into a 32x32 matrix. The weights or the indices where to put each value maybe learnt by the network itself. You would need to write the backprop for this though.


Pi/4 dqpsk carrier recovery I think we don't have to do carrier phase recovery. Do i need to do carrier frequency recovery. My symbol rate is high around 2 MHz  so i am thinking carrier freq offset is less compared the symbol rate. Or does my gardner timing recovery can take care of this carrier freq offset recovery. I have another question should i do differential decoding then timing recovery or timing recovery followed by differential decoding . Because the pi/4 dqps will have valid symbols at zero level also ie symbol with zero I data and Q non zero data. Similary zero Q data and I non zero data. So i am thinking if do timing recovery first i may get proper output as garder looks for zeros crossing between symbols. Is my explanation right

Yes you are correct, differential modulations such as DBPSK, DQPSK, $pi/4$ DQPSK do not need to perform carrier recovery as you can compare one symbol to the next to perform demodulation using a delay and multiply approach. This results in a simple and potentially low power receiver, at the expense of performance compared to coherent modulation schemes.

As far as a static carrier offset this should be no issue- the worst case is if the carrier frequency is changing from one symbol to the next (and at some level it always is, in particular phase noise)- in this regard simply consider how much your symbol rotates from one symbol to the next when you demodulate with this approach and you can directly see the impact of frequency offset (and if the SNR degradation from that rotation based on the resulting symbol distance would be acceptable). For example, a 2 degree rotation would be negligible, while with a 45 degree rotation you be on the threshold of an incorrect symbol decision with no other added noise (using the demodulation approach I describe below). The impact of this can be quantified from the phase noise, and how I approach that with a differential demodulation is to consider with regards to phase how a multiplier for small angles is a "phase adder" thus the delay and multiply demod process is a high pass filter with regards to phase noise. (The filter is given by $1-z^{-1}$). It has a gain of 2 at the high frequency portion (symbol rate) and a null at DC- so it passes (and doubles) the phase noise at those higher offsets). Ultimately you can integrate the high pass filtered phase noise from DC out to the symbol rate and double that to determine the rms phase error, and determine performance issues with consideration to the impact of phase offsets as I described above. What you will find in doing this is the demodulator will likely be very robust against phase noise.

An approach I would consider is doing a complex conjugate multiplication of each new symbol with the previous symbol with the additional +/-45 degree rotation in that process to remove the pi/4 offsets. The output of this would be I and Q with bi-phase signaling on each (similar to QPSK) that you could easily pass into a Gardner (at 2 samples per symbol) for timing recovery followed by matched filters and then data decision. An alternate approach used for low cost / low power non-coherent pi/4 QSPK demodulation consists of hard limiting the waveform and then using a simple FM discriminator followed by a four level detector (since the 4 phase transitions +45°, -45°, +135°, -135°  from each symbol to the next correspond to 4 possible frequencies).


How can I comprehend the complex coefficient of complex Fourier series? Is there a way to understand visually the complex coefficient of Fourier series (shown as $C_n$ in many sources)? I mean something similar to  this which represents $e^{j\omega_0 t}$: .

$\omega_0 \triangleq 2 \pi f_0 .$

But when I read $C_n  e^{jn\omega_0 t}$ my mind shuts down. Here is complex Fourier series: http://www.thefouriertransform.com/series/complexFourier.jpg.
And its one element is $C_n e^{j n \omega_0 t}$. How can I visualize this element? And does $C_n$ have a meaning besides being a coefficient?

if you think of $e^{j n \omega_0 t}$ as a "vector" (complex numbers are not quite the same as 2-dimensional vectors, but they look the same, and are both describable with two real numbers) with length 1 and starting out with an angle of 0 (pointing straight to the right, along the positive real axis, and this vector is rotating counter-clockwise at a rate of $n \omega_0$ radians per unit time.  this vector is what you see looking down the spiral in the graphic referred to in the question, looking along the time axis.

then multiplying by $C_n$ changes the length of the vector to $|C_n|$ and the initial angle of the vector is the angle of $C_n$ or $\arg\{C_n\}$.

that's how i would visualize it.


Compressive Sensing Incoherence Principle As people acquainted with Compressive Sensing would know, incoherence and sparsity are two main principles. I've been reading about compressive sampling and developed an interest into the topic. What I'm about to ask might be quite basic so sorry for that in advance.

As incoherence principle states, sensing basis and sparsifying basis must be incoherent. I can understand why it is so and I also understand why a random i.i.d. matrix would in fact give us a good incoherence. My question is, how would one determine it programmatically?

Let's say my sparse basis is the frequency domain. When one says sensing basis, is it actually the sensing matrix? Assuming that the signal is x, and sensing matrix is A, would I just calculate the coherence between fft of x and my sensing matrix? If yes, how would one do that, say, in MATLAB?

Thanks a lot.

Say you have measurements $y=Ax$, with a basis $\Psi$, and $A\in\mathbb R^{M\times N}$. Then we can write $y=A\Psi\hat x=\Phi\hat x$.

I then believe direct computation of the coherence requires us to solve
\begin{equation}
\mu = \max_{i<j} \frac{| \langle \Phi_i,\Phi_j\rangle |}{\| \Phi_i\|\| \Phi_j\|}
\end{equation}
where $\Phi_i$ is the $i$-th column of $\Phi$, which I think is of order $\mathcal O(N^3)$.

Now, assuming the above is correct you can compute $\Psi_i =\text{ifft}(e_i)$, where $e_i$ is a unit vector, construct $\Phi=A\Psi$, and then with a double for loop find $\mu$.


Regarding the nature of white noise in images I have just started learning the very basics of noise.

I have come to understand what white noise is in sound samples, but I do not understand the meaning of white noise in images.

In sound samples, time advances, and the various sound waves evolve through time. But images don't have time evolution, so I struggle to see how the frequencies come into play. Are the various waves evolving through one axis of the image?

How was the  white noise image generated on the wikipedia page on White Noise?

A signal is just a varying quantity over an axis.

Think of a greyscale image.

Think of a single row of that image. Can you see how the intensity evolves over the distance on the horizontal axis?

In an image, you get two axes (x and y), whereas in sound, you just get one.

In (mono) sound, your instantaneous intensity is the only quantity varying. In an image, that quantity might be brightness alone, or brightness, hue and saturation, or red, blue and green intensities, or whatever channels your color model defines.

Anyway, the point about white noise is it's autocorrelation function; and that's easy to define on the two axes of an image, exactly the same as for sound:

Rows of an image: If you shift the row(s) by any amount (other than 0), then the expected sum-of-products between the original and shifted row is 0.
Columns of an image: If you shift the columns(s) by any amount (other than 0), then the expected sum-of-products between the original and shifted column is 0.

If both properties are given (and in fact, the correlation in each direction, even if you "sliced" your image arbitrarily in diagonals), then your 2D noise is white.


Struggling to implement this simple algorithm for ECG wave delineation I'm trying my best but my maths isn't good enough to implement the algorithm as outlined in this paper in python. It for detecting the onset and offset of a wave on an ECG, and it's using a well validated method.

My data is in a numpy array.

The steps it uses are:

A. Computation of the envelope of the ECG

B. Computation of the auxiliary signal

C. Windowing

I think I've got step A. done:

def calculate_qrs_envelope(self):
    self.hilbert = np.imag(hilbert(self.lead_data['y']))
    self.envelope = np.sqrt(np.add(self.lead_data['y'],self.hilbert))
    self.plot(x=self.lead_data['x'], y=self.envelope, pen='b')

Where my data is stored in lead_data['y'].

If I plot it, it looks correct (I can't show the picture here as I don't have enough karma).

And I'm confident that I won't struggle with C.

However, B is difficult. It says in the paper:

I can calculate AS using a simple derivative, as so:

def calculate_auxiliary_signal(self):
    self.aux_sig = np.append([0],np.multiply(2,pow(np.diff(self.envelope),2)))
    self.plot(x=self.lead_data['x'], y=self.aux_sig, pen='g')

But I don't understand how to do it with regards to what it says about the parabolic fit.

My sampling frequency is 100 Hz.

Would anyone be able to give me a hand with this maths? It would be a great help.

I cannot find any sense on the formula exposed in the paper, but as is, the expression is just a derivative approximation, with $r=r_0$.

$$
x'_k=\frac 1{10}(2(x_{k+2r}-x_{k-2r})+(x_{k+r}-x_{k-r}))\\
y_k=2(x'_k)^2
$$

That is what the paper say.

That is, you have a signal $x$, and you calculate $x'$ as a moving average, taking samples above and below the current index.

For any given signal $x:x_1...x_n$, you will straightforwardly be able to calculate $x'$ and $y$ without trouble.


Equidistant log scale points I have an x-axis which is in logarithmic scale representing frequency, I would like to have NumPoints on this axis that are equidistant on the log scale, I have the following variables:

Fmin = 100
Fmax = 4000 or 8000
NumPoints = 10

What formula should i use to calculate where the NumPoints-2 points between Fmin and Fmax should go? (Fmin and Fmax should each have a point on them)

Additionally, when i switch from Fmax = 4000 -> 8000 i would like to keep the equidistant points already calculated between Fmin->4000 and add extra points between 4000 and 8000.

Summary

Generate NumPoints between FMan -> Fmax(4000), equidistant on a log scale.
Change Fmax -> 8000
Generate more points between 4000->8000.

Many thanks in advance

Andy

Take the log of your extremities, find your points as if it were on a linear scale beween these extremities, take the exponential of the result.

Example:

The log of your extremities are $\log 100$ and $\log 4000$. Your 10 points on a linear scale between these values would be: $\log 100 + k \frac{\log 4000 - \log 100}{9}$, $k \in \{0, 1, \ldots, 9\}$.

Taking the exponential: $e^{\log 100 + k \frac{\log 4000 - \log 100}{9}}$. Notice that the first and last values are 100 and 4000 as expected.

This works well with linspace in programming languages which have this primitive, for example with python/numpy:

numpy.exp(numpy.linspace(numpy.log(100), numpy.log(4000), 10))

The second part of your question is not clear. To generate more points by preserving the same spacing, continue increasing $k$ above 9.


Design a Filter with Resonance I want to design a Filter with a variable Cutoff Frequency and Resonance. However in the source I'm using there is a Bi-Quad Filter with two IIR and two FIR Filters designed like this :

Furthermore I found that I can compute the coefficients for the different Filter types as stated in this table :

Is there a way to alter those coefficients that i can compute the samples with a variant FIlter Resonance?

EDIT Since I'm designing this filter for a software synthesiser, with resonance I mean the typical boost or reduction of the cutoff-frequency in a variable manner.

First, you need to start with in the s domain where a 2nd order transfer function (2 pole low pass filter) which can be written as follows.

H(s) = s^2 + 2zws + w^2

Lets start with omega = 1 and zeta = sqrt(2)/2, then

H(s) = s^2 + 1.414 s + 1

This filter's response and pole locations are shown here.

To peak the response, simply modify zeta. Here zeta = 0.17 and omega = 1.

H(s) = s^2 + 0.34 s  + 1

Then to move the frequncy, change Omega.
Here zeta = 0.17 and omega = 0.39

H(s) = s^2 + 0.1326s + 0.1521

Now all that is left is to determine the z domain coefficients. To do this, use the bilinear transform as shown here.

For the the last example I gave, use  A=1, B=0.136, C=0.1521, D=E=0, F=1, T=2 to calculate the a's and b's.

See this page for a derivation of the equations given here. Here is a rather simple tool that can generate any 2nd order response you want.


Why is $\sin(t)$ a stationary process? I am trying to understand the meaning of the term Stationary Process.
For example, I was told that $\sin(t)$ is a stationary process.

Could someone try to explain, in simple words, why is $\sin(t)$ (for example) is a stationary process?

$\sin(t)$ is no random process, because there's nothing random about it. You could add a random amplitude to get a random process:

$$x(t)=A\sin(t)\tag{1}$$

This is a random process because $A$ is a random variable. However, $x(t)$ is not stationary, but it is cyclostationary, i.e., its statistical properties vary periodically. You can make the process $x(t)$ stationary by adding a random phase:

$$\tilde{x}(t)=A\sin(t+\phi)\tag{2}$$

The phase $\phi\in [0,2\pi]$ is a uniformly distributed random variable that is independent of $A$. It can be shown that the statistical properties of $\tilde{x}(t)$ given by $(2)$ are independent of $t$, and hence, the process is stationary.


Dimensional analysis of integrated white noise process This question is somewhat related to this post.

Let us consider we have a white noise current source $i_n(t)$, with a variance $\sigma_i^2$, and mean, $\mu_n=0$.

Assume that this current is passed through a system, which integrates for a finite time duration ($T_{int}$) on a capacitor $C$. One can then write an equation for the power of the noise voltage (integration of current produces voltage $V=\frac{1}{C}\int_{0}^{T}i(t)dt $  ) produced from this integration process ($\sigma_v^2$):

$$\begin{align}\sigma_{v}^{2} &=E\bigg(\frac{1}{C}\int_{0}^{T_{int}}i_{n}(t_{1})dt_{1}\frac{1}{C}\int_{0}^{T_{int}}i_{n}(t_{2})dt_{2}\bigg)
\\&=\frac{1}{C^{2}}\int_{0}^{T_{int}}\int_{0}^{T_{int}}E(i_{n}(t_{1})i_{n}(t_{2}))dt_{1}dt_{2}
 \\&=\frac{1}{C^{2}}\int_{0}^{T_{int}}\int_{0}^{T_{int}}\sigma_i^2\delta(t_{2}-t_{1})dt_{1}dt_{2}
 \\&=\frac{1}{C^{2}}\int_{0}^{T_{int}}\sigma_i^2dt_{2}
 =\frac{1}{C^{2}}\sigma_i^2\cdot T_{int}\end{align}$$

Here $\delta$ denotes the dirac delta impulse function.

This is a similar result obtained to what was obtained in the post linked above. The problem with the above equation is that it is dimensionally inconsistent (there should be a $T_{int}^2$ factor in the numerator to make it consistent.)

Can someone kindly point out the error in the above argument?

The problem comes from the dimension of $\sigma^2_i$, which you assume to be $[I^2]$. That's, however, not the case if you define

$$E\{i(t_1)i(t_2)\}=\sigma^2_i\delta(t_2-t_1)\tag{1}$$

Note that the dimension of a Dirac delta impulse is the inverse of the  dimension of its argument, which is $[1/T]$ or $[Hz]$. So you get for the dimension of $\sigma^2_i$ $[I^2T]$ or $[I^2/Hz]$. You can also see this by looking at the power spectral density (the Fourier transform of $(1)$), which is simply $\sigma^2_i$, and the dimension of which must be $[I^2/Hz]$. With this dimension of $\sigma^2_i$, the dimension of your result is correct.


Goertzel algorithm: Relationship of magnitude I've written a quick test app that uses the Goertzel algorithm to determine if a given frequency is present in the signal. This is to pick up DTMF tones and various other signals. The app appears to be working. I'm using the 'magnitude squared' approach in my final stage.

http://www.embedded.com/design/configurable-systems/4024443/The-Goertzel-Algorithm

Now, how does this magnitude relate to all the other variables (especially the amplitude of the source signal)? Most of the implementations determine the presence of a given frequency if the magnitude exceeds some threshold. I want to know how this threshold is determined, and how would I go about selecting a suitable threshold.

EDIT:

OK, I think I've found the answer (I hope). Most implementations of the algorithm state that is will return the 'relative power'. I couldn't find any that stated what it was relative to.

This morning, it dawned on me (pun intended) that it's probably relative to the input signal (duh!). So, if we calculated what the power of the input signal is over all the samples taken, and then compare this to the relative power given by the algorithm, we should have a good idea of the contribution of our target frequency?

Is my understanding correct?

You can potentially base the threshold on the overall energy in the frame. If the tone has, say, 10% or more of the total energy, there is a good bet it's there and not just noise.

In general FFT based algorithms may not be the best choice for this. Unless your tone and your local sample clock are phase-locked, you can't line up the tone on one of the FFT bins (or your Goertzel frequency) so the energy of the tone will "smear" out over multiple frequency bins. The FFT is actually a pretty band band pass filter.

The alternative would be a time-domain band pass filter with a well defined bandwidth and steepness. This would also get rid of the framing constraint and you can dial in the time constant based on the properties of the signal rather then internal properties of the algorithm.


complementary IIR filters I would like a pair of complementary IIR filters (lowpass/highpass). By complementary, I mean, when the output from the two filters is summed, the original signal is recovered. I thought I could build such pairs with butterworth filters but using a little math, I discovered only 1st order filters were complementary. I thought I'd done this before, but I'm forgetting how.

Is something wrong with my math? Is there an easy solution I'm forgetting about?

Thanks!

The frequency response of the two complementary filters are $H_2(e^{j\theta}) = 1 - H_1(e^{j\theta})$, or the impulse responses $h_2[n] = \delta[n] - h_1[n]$.

For an IIR filter, $H_1(z)$ can be written as $\frac{b_0 + b_1 z^{-1} + \ldots}{a_0 + a_1 z^{-1} + \ldots}$. Then $H_2(z)$ should be something like $\frac{(a_0 - b_0) + (a_1 - b_1) z^{-1} + \ldots}{a_0 + a_1 z^{-1} + \ldots}$.

So the non-recursive coefficients for $H_2$ are now $(a_0-b_0)$, $(a_1 - b_1)$, etc.

The recursive coefficients are the same for both filters.


How is the symbol error rate for M-QAM, 4QAM,16QAM and 32QAM derived? How do you derive the theoretical symbol error rate as a function of $E_\mathrm{b}/N_0$ for 4QAM? I know that the result should be $Q\left(\sqrt{2E_\mathrm{b}/N_0}\right)$ but I am ĺooking for the derivation. Also, what are the symbol error rates vs $E_\mathrm{b}/N_0$ for 16QAM and 32QAM?

In $2^{2n}$-QAM with a square constellation, there are $4$ "corner" points and
$4(2^n-2)$ "edge" points, and $(2^n-2)^2$ "interior" points.  The conditional
symbol error probabilities given that each type of point is transmitted,
are
$$\begin{align}
P_e(\text{corner}) &= 2Q(x) - Q^2(x)\\
P_e(\text{edge}) &= 3Q(x) - 2Q^2(x)\\
P_e(\text{interior}) &= 4Q(x) - 4Q^2(x)\\
\end{align}$$
where $Q(x)$ is the complementary cumulative probability distribution function
of the standard Gaussian random variable. Combining these using the law of
total probability (with the assumption that all $2^{2n}$ signals are equally
likely) gives
$$P_e\left(2^{2n}\text{-QAM}\right) = 4 \left[1 - 2^{-n}\right]Q(x)
- 4\left[1 - 2^{-n}\right]^2Q^2(x)$$
For $4$-QAM, where $n = 1$, this reduces to
$P_e(4\text{-QAM}) = 2Q(x) - Q^2(x)$. For more details, see, for
example, pp. 147-153 of this ancient lecture note  that I once wrote.