FFT, How to decide if there is a signal among noise? I have sets of data of different deep sky objects. My job is to check for any periodicity.

I use IDL to run an FFT and wavelet methods to check for a signal. To test my code I ran the IDL built in FFT operation on a sine wave, my code works.

I then add noise to the sine wave and increase the noise until I cant distinguish the frequency peak in the FFT output.

My question is what mathematical method can I use to say 'yes there is a signal and this is its frequency' or 'no there is no signal' ? Can or should I get a percent of certainty ?

This post is very close to what I'm asking, but I only have one data set for an object so how could I make a ROC curve?: What statistic is used to determine presence of a signal in noise?

Also I think I dont want a ROC curve.
I will respond to any question you may have about my question. Thank you very much for your help. I am currently do research at UC Santa Cruz as an undergraduate.

As a fellow astrophysicist, I suggest you use the standard 3 sigma above noise for detection, and 5 sigma above noise if you want to do analysis on what you find.

If your noise does not have a flat background, you need to subtract a "Flat", an estimation of that the background might be. If you know what the noise should look like, such as exponential, gaussian, etc., try subtracting such a curve. If not, try a first order polynomial. If not good enough, then go with a second order.

Beyond that, it becomes interpretation... What if that tendency you notice in the noise is actually a signal, spread over many frequencies?

Try computing the FFT of a pure white noise image, to get an idea of what your noise is supposed to look like. (Off course, remember that FFTs assumes your signal is periodic beyond its edges. Don't forget to zero-pad your data, and use a correct apodization, or limit your analysis to something smaller than half the size of your image.)


Short-time Fourier transform Tradeoffs In the Short-time Fourier transform what are the tradeoffs between different:

window functions
window lengths
sample rates

Window functions
Different window functions introduce different amounts of scalloping noise/artifacts. Scalloping noise occurs when the frequency bins of a DFT don't optimally align with the DTFT.

The shape of the DTFT corresponds to the Fourier transform of the window function, and thus window functions whose Fourier transform has shorter side-lobes will have less scalloping noise. However, window functions with less scalloping noise are narrower in the time domain (the Rectangular window has the most scalloping noise) and thus, without window overlapping, will lose the frequencies at the edges of the window.
Terminology
Since the frequency content of overlapping windows are displayed sequentially rather than averaged, I find talking about overlapping windows confusing. I would rather divide the time domain into Frames for which we try to obtain a picture of the frequencies in their range as accurately as possible. This may involve extending the window function past the range of the frames to minimize the scalloping noise/windowing loss tradeoff. However, this will introduce leakage of the frequencies between neighboring frames (reducing time resolution).
Window lengths and sample rates
Edit:
My answer here was bad and I replaced it with this:

Frequency resolution is how similar the window function's magnitude in frequency space is to the Dirac delta function. This is because the product of the window and the signal in the time domain becomes convolution in the frequency domain (and a convolution with the Dirac delta function is a sampling which would give perfect frequency resolution) The fatter the mainlobe (quantified by it's variance), and the higher the sidelobes, the worse the frequency resolution. Additionally, Time resolution can be quantified a the variance of the window function in the time domain.

From wikipedia https://en.wikipedia.org/wiki/Short-time_Fourier_transform#Resolution_issues:

The product of the standard deviation in time and frequency is limited. The boundary of the uncertainty principle (best simultaneous resolution of both) is reached with a Gaussian window function, as the Gaussian minimizes the Fourier uncertainty principle.

The key to note is that Praat's Gaussian Window is twice as long as it's other window functions. In the STFT you can extend you can have an extremely long window function compared to the time interval you want a view of the frequencies for, so long as the window function is close enough to zero for outside of your desired interval. This is optimal. Use the adjustable windows https://en.wikipedia.org/wiki/Window_function#Adjustable_windows
Zero padding
Zero padding doubles the number of bins while keeping the maximum frequency the same, interpolating the DTFT:

From Dan Boschen's answer What happens when N increases in N-point DFT


Sending a real number over a noisy channel with limited energy budget Note: My question is theoretical and my goal is to apply the intuition from the answer to a different field.

Imagine I have a real number $X$ which I want to send it over some noisy channel. (The metric that should be optimized is probably squared error).

I have a limited amount of energy $E$ that I can use to send the signal.

What would be the distribution of the received number $R$ and how would it be affected by E?

Just an example: A kind of formula I'm expecting as an answer ():

$R(E) = X + Noise/E$.
$R(E) = X + Noise/\sqrt E$.
$R(E) = X + Noise/ln(E)$.
something else...

I want to know the correct formula. The most important part for me is the effect of $E$.

P.S. I know that I've omitted some information that might turn out useful. If that's true, I'd like to know how the answer is affected by that information. The noise is... say, some "common" type of noise. The transmission/encoding/correction method is the optimum one given the optimization metric (SE) and the noise. (What would that method be BTW, for some common types of noise?) Does the answer change is the time to send a signal is limited?

Update: Trying to state the question in a more formal way:

We're transmitting a real number.

There is a metric $loss(x_{sent}, x_{received})$ that the transmission should minimize. The $loss$ is a non-negative, symmetric and monotonously increases when  $|x_{sent} - x_{received}|$ increases. An example would be $loss(x_{sent}, x_{received}) = (x_{sent} - x_{received})^2$

There is some transmission algorithm (encoding, transmission, decoding, error correction) that uses energy to send the signal over the noisy channel/medium.
The algorithm must be optimal in that it minimizes mean $loss$ when transmitting over a given noisy channel and using the specified amount of energy $E$: $best\_transmitter:=argmin(mean(loss(x, transmit(x))))$

What would be the relation between the following distributions?: $$(best\_transmitter(x, \alpha*E) - x)$$ and $$(best\_transmitter(x, E) - x) $$

If I spend 10 times more energy on transmission, would the error be 100 times smaller, 10 times smaller, 3.1 times smaller or 2.3 times smaller? (given the optimal transmission algorithm)

The terms of the kind $N(0, 1)*\alpha/E$ are simply representing noise to the system and the way one defines the multiplicative factor $\alpha/E$ simply changes the variance of the gaussian $N(0,1)$ without impacting the mean.

Note: the mean of $R(E)$ is $X$

So the distribution of $R(E) = X + N(0,1)\alpha/E$ is $N(X,\frac{\alpha^2}{E^2})$ and for the received $R(E) = X + N(0,1)\alpha ln(E)$ is $N(X,\alpha^2ln(E)^2)$.

In the first case as the energy of the signal is increasing the effect of noise is decreasing (it's variance). So it makes more sense compared to second one where the variance on noise is increasing with increasing signal energy. However it does not reflect the true SNR

Edit: The correct way with the update in your question is $\sqrt(E)$ term, this is because this truly represnts the actual SNR, and this is what should be used to decide the minimum error detection rule


how were DTMF frequencies determined exactly and how can I extend them? I've been playing with DTMF frequencies and the Goertzel algorithm, and would like to know exactly HOW the DTMF frequencies were determined. I understand that they were chosen such that none of the frequencies would be equal to sum and difference frequencies, or intermodulation products, but I can't find any information on whether the frequencies were simply determined through experimentation, mathematical calculation, etc… does anyone know more about this?

Also, I would like to understand how to extend the two sets of frequencies to two sets of 8 frequencies, or 3 sets of 4 frequencies, such that it's possible to represent 64 keys. I'm not sure whether it would be easier to work with 3 sets of 4 versus 2 sets of 8 each. I guess the analysis would be easier if just using 2 sets of 8? There would be less intermodulation products to keep in mind …

Thanks!
B

EDIT: if it's too difficult to track back how this was achieved, i would like to know what the best way is to find 2 sets of frequencies, giving 64 unique combinations, and such that if you pick 2 frequencies from each of the 2 sets randomly, and mix the sine waves together, and put them through some non linear distortion, they are still easily detectable using goertzel, and the other frequencies are not "hit" (their bins remain at a normal background noise level)

From page 293 of Engineering and Operations in the Bell System (1982-1983) by AT&T Bell Labs:

  The tones have been carefully selected to minimize harmonic interference and the probability that a pair of high and low tones will be simulated by the human voice, thus protecting network control signaling.

P.S. Yes, I know this answer is 6 years late, I just happened to be researching for the answer and I thought this could help.


Can anyone explain me this MDFT polyphase filter bank? I found a very interesting piece of code for MDFT polyphase filter bank here.
Unfortunately, there doesn't seem to be a paper describing the theory. Does anyone know some reference for the code?
I'm especially interested in these 3 topics:

What are the data in the channels exactly? Are they supposed to be real, imaginary or complex?
The code computes only half of the number of bands. Is that due to real-valued signals used?
The result of the synthesis stage is build as a channel-wise difference of the outcome of the two synthesis filter banks. Why is it done like that? I can't find any paper describing this idea.

It's input agnostic, anything will work just as it would with any other real valued prototype filter. I've implemented polyphase filter banks this on radar systems in practice, where we're operating on complex data, both pulse compressed and uncompressed. Filter banks like these have loads of applications due to the i inherent design and theoretical speed.
Polyphase filter banks typically create a phase shifted "copy" of the original version to achieve better reconstruction. From just scanning the code, it looks like this is what he's doing with his X1 and X2 vectors. You'll notice if you step through the code and check out the spectrum, the sub-bands will have no passband overlap. It's even more obvious if you look at the shifted versions of the filter frequency response. Since the "copy" of the original signal was phase shifted, it will also have non-overlapping segments within itself, but tougher the original and copy combine to fully cover from -pi to pic on a normalized frequency scale.

Typically the processing chain is this: get input, create a copy, phase shift the copy (everything for here on out is done twice, one for the original and one for the phase shifted copy; let's call these upper (original) and lower (copy)), apply prototype polyphase filter to the upper and lower signals, apply a DFT to do the polyphase magic, and now you're all channelized. Synthesis is pretty much just this in reverse.

The part people struggle with is typically the polyphase filters and the use of the upper and lower signals, which seems to be the case for this question . The filter bank itself is not complicated but the math behind the polyphase can be if you don't have a multi-rate signal processing background . The prototype filter is typically designed to make say 32 sub-bands, but since we have are upper and lower signals which each have 32 sub-bands, we really have 64. I'll leave the math behind why polyphase filters work to a textbook since that wasn't in your question.

See above, has to do with the phase shift.  Also go a few lines up and you'll see a negation. Check out those indices. The subtraction you're seeing may not be doing exactly what you think it's doing. Try stepping through the code and checking it out.


How did they get this table? I'm trying to work out how this table has been calculated in the attachments under part c. I thought it was $h[0]\star x[n]$ on the 4th row (first row with no label) and then $h[1]\star x[n-1]$ on the row below and then $h[2]\star x[x-2]$ but this doesn't get me the same values.

Would anyone be able to tell me?

Thanks

The point is that

  The output of an LTI system to an input is a weighted sum of shifted versions of the impulse response.

You can see that the input is $$x=(1)\delta[n]+(-1)\delta[n-1]+(1)\delta[n-2]$$
The response to $A\delta[n-n_0],\, \forall n_0$, is $Ah[n-n_0]$. By the principle of superposition, the output to a sum of shifted $\delta$ terms is a sum of shifted $h$ terms. Here:

$$y[n]=(1)h[n]+(-1)h[n-1]+(1)h[n-2]$$
where $h[n]$ is $$ \cdots \quad 0 \quad 1 \quad -2 \quad 3 \quad -4 \quad 2 \quad0\quad\cdots$$


SVD vs matched filter Reading about singular value decomposition (SVD) in the context of signal processing applications, one can separate the signal from the noise into orthogonal subspaces. On the surface this sounds like using the decomposed output of the SVD is more optimal than matched filter (or correlation) for detection, demodulation, timing estimation, etc. Is this true or is there more nuance to this problem?

That's not true, it's not better.
The thing is: the matched filter just implements the projection in the signal vector space, onto the signal vector itself (or a multiple thereof). (You'll find correlation is just an inner product in that space.)
The line through that vector is the signal subspace, the plane to which that vector is normal is the noise space. Matched filtering is a method of linear algebra, and not different than using the SVD to find the subspaces. It's just that for matched filtering, the signal space is known.


Chebyshev Filter Transfer Function I'm trying to derive the transfer function for Chebyshev filter.
$$|H(\Omega)|^2=\frac{1}{\sqrt{(1+\epsilon^2T_n^2(\frac{\Omega}{\Omega_c})}}$$
where   $$T_n(x)=\cos(N\cos^{-1}(x)) \forall x \le 1$$
$$T_n(x)=\cosh(N\cosh^{-1}(x)) \forall x \ge 1$$
$$H(s)=\frac{1}{\sqrt{(1+\epsilon^2T_n^2(\frac{s}{j\Omega_c}))}}$$
To calculate the poles I made the following substituion which looks like this:
$$\cos^{-1}(\frac{s}{j\Omega_c})=\alpha + j \beta$$
So
$$s=j\Omega_c(\cos\alpha \cosh\beta - j \sin\alpha \sinh\beta)$$
$$s=\Omega_c \sin\alpha \sinh\beta + j \Omega_c \cos\alpha \cosh\beta$$
where $\alpha=\frac{(2k-1)\pi}{2N}$ and $\beta = \frac{1}{N} \sinh^{-1}{\frac{1}{\epsilon}}$
But with this substitution my poles yield to be stability or produce an ambiguity on the stability end.
But if I make the substitution as $$\cos^{-1}(\frac{s}{j\Omega_c})=\alpha - j \beta$$
poles seem stable.
Can someone help me with this. is my substitution wrong?

Note that the squared magnitude of the frequency response is given by
$$\big|H(j\omega)\big|^2=\frac{1}{1+\epsilon^2T^2_N\left(\frac{\omega}{\omega_c}\right)}\tag{1}$$
In the $s$-domain we have
$$H(s)H(-s)=\frac{1}{1+\epsilon^2T^2_N\left(\frac{s}{j\omega_c}\right)}\tag{2}$$
Computing the zeros of $(2)$ does not only result in the zeros of the filter's transfer function $H(s)$, but also in the zeros of $H(-s)$. The zeros of $(2)$ lie on an ellipse, and for each zero $s_0$ in the left half-plane, there's also a zero $-s_0$ in the right half-plane. Since we want a stable filter, we assign all zeros in the left half-plane to $H(s)$.


Using Convolution as Feature Extraction I  am now studying image processing in my spare time. My understanding of convolution is about 'response to a specific filter':

When we have a raw image, or raw signal; and a filter, aka kernel; we apply a 'moving dot product' between the image and the filter.

The filter/kernel represents our domain-knowledge about the raw signal, which is usually some kind of features or patterns. For image processing, such features are often edge, corner, SIFT, etc.

When we apply the "moving dot product", we practically divide the image into "patches" and ask each patch of the image: "Do you look like the feature I am looking for?", high response value means "yes, I am the feature (edge) you are looking for"; while low response value means "no, I am not"

My first question is : is my understanding correct?

My second question is: If my understanding is correct, then we should be able to design much more fancy kernels which represents some very specific features: such as eye (dark in the middle, surrounded by white areas), or nose (triangular shape). But I haven't found a lot of literatures in the area of defining your own kernel.

Can anyone provide any insight here?

Thanks!

Your understanding of the convolution process is correct. However, note that in the convolution the kernel is first mirrored before the dot product. Though, for a symmetric kernel this does not matter.

Understanding the convolution in your way, of finding specific features, is a similar interpretation as that of a Matched filter. Furthermore, using the convolution to find a specific part of the image (e.g. eye), is sometimes called Template matching (but here, usually correlation is used, which is equal to convolution except the kernel is not mirrored).

Another interpretation of the convolution is done via the convolution theorem, which says that in the frequency domain, the frequency response of the filter and the image are multiplied. Essentially, you can then understand the filter to enhance / reduce higher or lower frequencies, (or do more fancy stuff in the frequency domain). With this interpretation, it sometimes becomes more understandable why a certain kernel has a certain effect (e.g. the Gaussian filter is easier to understand as a lowpass than a filter to look for Gaussian shapes in the image).


Check whether a given equation is linear $ a = (x, y) \in \mathbb{Z}^2 $ is given as a pixel. My equation in which $g$ is image(matrix)  is defined as,

$f(x, y) = 56g(x,y)+93g(x−1,y)+92g(x+1, y)−57g(x, y−1)+555g(x, y+1) $

How can we know whether $f(x,y)$ is linear, proving? By saying linear, it is like a filter which is linear. $g(x,y)$ means a pixel at $x$ and $y$. I'm really newbie for image processing and the question will help me a lot what's going on underneath.


Implementing an L-System conversion from R2 to R1 for symmetric, square matrices I have been reading for some time about L-Systems, and specifically the Hilbert Space filling curve.  I am interested in writing a function to convert upper-triangular matrix coordinates into an ordered, 1-dimensional space filling curve.

The usual Hilbert space filling curve is derived with a single, first-order curve:

Then

I have successfully implemented the mapping algorithm listed on Wikipedia,

//rotate/flip a quadrant appropriately
void rot(int n, int *x, int *y, int rx, int ry) {
    if (ry == 0) {
        if (rx == 1) {
            *x = n-1 - *x;
            *y = n-1 - *y;
        }

        //Swap x and y
        int t  = *x;
        *x = *y;
        *y = t;
    }
}

//convert (x,y) to d
int xy2d (int n, int x, int y) {
    int rx, ry, s, d=0;
    for (s=n/2; s>0; s/=2) {
        rx = (x & s) > 0;
        ry = (y & s) > 0;
        d += s * s * ((3 * rx) ^ ry);
        rot(n, &x, &y, rx, ry);
    }
    return d;
}

//convert d to (x,y)
void d2xy(int n, int d, int *x, int *y) {
    int rx, ry, s, t=d;
    *x = *y = 0;
    for (s=1; s<n; s*=2) {
        rx = 1 & (t/2);
        ry = 1 & (t ^ rx);
        rot(s, x, y, rx, ry);
        *x += s * rx;
        *y += s * ry;
        t /= 4;
    }
}

However, I need to add two more first order curves to this process: one that is a right angle and the other which is a wedge.  These two new curves do not need to be rotated or flipped, and in combination with the usual Hilbert 1st order curve, one can fill an upper triangular matrix (see my sketch below):

That is, the upper left of the red right-angle again becomes a right angle, the upper right becomes a Hilbert curve, and the lower right part of the right-angle becomes a wedge, and so on.

I would like to implement the above proposed L-System conversion, but every single tutorial on L-Systems is about STARTING AT ZERO and then DRAWING A LINE that fills the space without crossing itself (i.e. writing a string of steps to trace the whole space, instead of converting a subset of coordinates in the space to there 1-D counterpart).

Can anyone offer any intuition on how I might incorporate these new restriction into the existing code?  Are there any resources that explain L-System conversions (as opposed to drawing systems), and how to construct them in code??

I've defined the following function in R

State 0 indicates that one is located within a right-angle region

State 1 indicates that one is located within a Hilbert curve

State 2 indicates that one is located within a wedge region

rot <- function(n, x, y, rx, ry) {
  if (ry == 1) {
    if (rx == 0) {
      x <<- n-1 - x
      y <<- n-1 - y
    }

    temp <- x
    x <<- y
    y <<- temp
  }
}

xy2d <- function(n, x, y) {
  # rx = x region (right = 1, left equals 0)
  # ry = y region (top = 1, bottom = 0)
  # s = half the dimension
  # d = linear coordinate value
  d=0L
  state = 0L
  for (s in 2L^(log2(n/2L):0)) {
    rx = bitwAnd(x, s) > 0
    ry = bitwAnd(y, s) > 0
    if (state == 0) {
      # *--2  | 0  ->  *--2   | 1 -> 0--3   | 2 -> *--2   |
      # 0--1  |        0--1   |      1--2   |      1--0   |

      # if x on right region, add lower left triangle (region 0)
      #      (s * (s + 1) / 2) * rx
      # if y on upper region, add lower right square as well (both regions 0 and 1)
      #     (s * s) * ry
      d = d + (s * (s + 1) / 2) * rx + (s * s) * ry
      d
      # change the state
      #  rx  ry | (rx + ry)
      #  0   0  |  0
      #  1   0  |  1
      #  1   1  |  2
      state = rx + ry
      state
    }else if (state == 1) {
      # 0--3  | 0  ->  0--1   | 1,2 -> 0--3   | 3 -> 2--3   |
      # 1--2  |        3--2   |        1--2   |      1--0   |

      # Determine how many squares to add
      # equals number in s-by-s square times it's order in the curve
      #  rx  ry | (3 * rx) ^ (1 - ry)
      #  0   1  |  0 ^ 0 = [00] ^ [00] = 00 = 0
      #  0   0  |  0 ^ 1 = [00] ^ [01] = 01 = 1
      #  1   0  |  3 ^ 1 = [11] ^ [01] = 10 = 2
      #  1   1  |  3 ^ 0 = [11] ^ [00] = 11 = 3
      d = d + s * s * bitwXor(3 * rx, 1 - ry)
      d
      # ROTATE
      rot(n, x, y, rx, ry)
      x
      y

    }else if (state == 2) {
      # *--2  | 0  ->  2--1   | 1 -> *--2   | 3 -> *--2   |
      # 1--0  |        3--0   |      1--0   |      0--1   |

      d = d + (s * s) * (rx == ry) + (s * (s + 1) / 2) * (rx == 1 && ry == 1)

      #  rx  ry |  Target | state - (rx + ry)
      #_________|_________|__________________
      #  1   0  |  1      | 2 - (1 + 0) = 1
      #  0   0  |  2      | 2 - (0 + 0) = 2
      #  1   1  |  0      | 2 - (1 + 1) = 0

      # new_state = current_state - (rx + ry)
      state = 2 - (rx + ry)
      if (state == 1) {
        # If the state goes from state-2 to state-1, then rotate the hilbert
        # curve 180 degrees (i.e. reflect horizonally and vertically)
        x <- n-1 - x
        y <- n-1 - y
      }
    }
  }
  return(d)
}

xy2d(8,5,0)

[1] 16
```


Integer and a half periods of sine wave in DFT window Why does a pure sine wave starting and ending at 0 in the sample window of a DFT but containing 1.5 cycles (or periods) produce the same spectral leakage (i.e. not more) as a pure sine wave with integer number of periods in the sample window of a DFT.

That's not true. Take a simple example:

N = 8; n = 0:N-1;
w1 = pi/4; w2 = 3*pi/8;
s1 = sin(n*w1); s2 = sin(n*w2);
S1 = fft(s1); S2 = fft(s2);
subplot(2,2,1),stem(n,s1)
subplot(2,2,2), stem(n,s2)
subplot(2,2,3), stem(n,abs(S1))
subplot(2,2,4), stem(n,abs(S2))


Oversampling in quantization Q:  It is said that "to maintain the same quality in the two cases, we require that the power spectral densities remain the same". Why is this a measure of the same quality? Why is not
the integral of  each power spectral density over the corresponding Nyquist interval not used, instead?

Let $x(t)$ be your continuous-time bandlimited signal, modeled as a WSS random process with a PSD of $S_{xx}(\Omega)$ in the band $\Omega \in[-W,W]$, $\Omega$ in radians per second.
Sampling $x(t)$ at the critical rate $T_s = \frac{\pi}{W}$, yields the unquantized, discrete-time sequence $x[n] = x(nT_s)$. Let the associted PSD with $x[n]$ is
$$\Phi_{xx}(e^{j \omega}) = \frac{1}{T_s} S_{xx}(\frac{\omega}{T_s})  ~~~,~~~\omega \in [-\pi, \pi] \tag{1}$$
Assume some $N$ bits of quantization is used with a step-size $\Delta$. The process of quantization converts the unquantized sequence $x[n]$ into quantized (digital) $\hat{x}[n]$ with the quantization error denoted as $e[n]$ such that: $$\hat{x}[n] = x[n] + e[n] \tag{2}$$
For simplification reasons, the quantization error $e[n]$ is assumed to be an i.i.d., white-noise of independent (from $x[n]$ and itself at nonzero lags) WSS random process uniformly distributed in $[-\frac{\Delta}{2}  ~,~ \frac{\Delta}{2}]$. Which, therefore has the PSD given by :
$$ S_{ee}(e^{j\omega}) = \frac{\Delta ^2}{12} ~~~,~~~\pi < \omega < \pi \tag{3}$$
Notice, how I modeled the quantization noise process, directly in discrete-time, with a power spectral density independent of sampling period $T_s$. This means that whatever sampling rate is used to obtain the quantized samples $\hat{x}[n]$, the assoiced noise power at the discrete-time will be the same as $\sigma_e^2 = \Delta^2/12$ in the range $-\pi < \omega < \pi$. It only depends on the number of bits; hence the step-size $\Delta$.
However, the same is not true for the PSD of the signal $x[n]$. As the sampling rate changes, its PSD will scale both in amplitude and in the frequency band as given by the relation. Assume we employ $M$ times oversampling so that $T' = T_s/M$, then:
$$ x[n] = x(n T_s) \longleftrightarrow \Phi_{xx}(e^{j\omega})  \implies \\ x_o[n]=x(n \frac{T_s}{M}) \longleftrightarrow  \begin{cases}{ M \cdot \Phi_{xx}(e^{j\omega/M}) ~~~,~~~|\omega| < \frac{\pi}{M} \\ 0 ~~~,~~~ \frac{\pi}{M} < |\omega| < \pi }\end{cases} \tag{4}$$
Now, if you lowpass filter the quantized oversampled sequence $x_o[n]$ and decimate the result back into its critical rate, both PSDs of $x_o[n]$ and $e[n]$ will be reduced by $M$, however, this will only revert PSD of $x[n]$ into its original, while reducing the equivalent noise power of quantization into $1/M$ of its original value, thus enabling one to use more bits to represent the original signal samples $x(nT_s)$.


Introduction to speech processing with code and examples? Can anyone suggest an introduction to speech processing,
class notes or books
accompanied by programs that can be run to see what the algorithms do ?
For example, in Burrus et al., Computer-Based Exercises for Signal Processing Using Matlab,
20 of the 400 pages are on speech processing. But it's from 1994;
there must be other hands-on books ?
Books with reviews like the following may have their uses, but are not what I'm looking for:

  I found this book comprehensive but incomprehensible, primarily because of the
  lack of real-world examples. It would benefit from a series of programming exercises
  with training sets

(This question may be a tall order — not many people ar good at theory and practice and explaining.)

(Should this be a community wiki ?)

You can try Prof. Rabiner's slides at http://www.ece.ucsb.edu/Faculty/Rabiner/ece259/ . They also have some matlab code examples which you can try.


How is the DTFT of a periodic, sampled signal linked to the DFT? I am trying to understand the connection between FT, DTFT and ultimately the DFT. But I am failing to link the DTFT to the DFT.
This is how far I am getting: Say I have a function $f(t)$, and its Fourier Transform $F(t)=\mathcal{F}\{f(t)\}(\nu)$ defined as $F(t)= \int_{-\infty}^{\infty} dt \, f(t)\exp(-i \,2\pi\nu t)$. Lets then say we have ideally sampled the signal at time intervals $\Delta T$ like so:
$$\bar{f}(t)=f(t)\sum_{n=-\infty}^{\infty} \delta(t-n\Delta T).$$
I understand that the FT of this signal leads me to the DTFT:
$$\bar{F}(\nu) = \mathcal{F}\{\bar{f}(t)\}(\nu) = \sum_{n=-\infty}^{\infty} f_n \exp(-i\,2\pi \nu n\,\Delta T), \, f_n = f(n \, \Delta T)$$
I also understand that $\bar{F}(\nu)$ is linked to the (continous time) Fourier Transform via:
$$\bar{F}(\nu) = \frac{1}{\Delta T}\sum_{n=-\infty}^{\infty}F\left(\nu - \frac{n}{\Delta T}\right),$$
which is a periodic function with period $1/\Delta T$. If $f(t)$ is bandlimited then the values of the DTFT at discrete frequencies $\nu_m = m/M \cdot \Delta T, \, m=0,1,...,M-1$ will give me a value proportional to the FT of $f$ at these frequencies. I also understand that I only need $M$ samples in Fourier space due to the periodicity of $\bar{F}$
Where I am struggling is to understand how to get from here to the DFT. I'll plug the discrete frequencies into the DTFT above and I further assume that the sequence $f_n$ is $M$-perdiodic such that $f_{n-kM}=f_{n}, \, k\in \mathbb{Z}$. Following this wiki article I write:
$$\begin{eqnarray}
\bar{F}(\nu_m)=\bar{F}_m &=& \sum_{n=-\infty}^{\infty} f_n \exp\left(-i \,2\pi \frac{m}{M}n\right) \\
&=& \sum_{k=-\infty}^{\infty} \left(\sum_{n=0}^{M-1}\exp\left(-i \,2\pi \frac{m}{M}n\right) f_{n-kM}\right) \\
&=& \sum_{n=0}^{M-1} \exp\left(-i \,2\pi \frac{m}{M}n\right)  \sum_{k=-\infty}^{\infty}f_{n-kM}
\end{eqnarray}$$
The wikipedia article then calls the inner sum $\tilde{f}[n]$, which leads us to the DFT. But I fail to see how this infinite sum of a periodic series relates to the cyclic extension of a finite series. I just don't see it.
I know that this topic is a perpetual source of confusion for people and we have excellent answers e.g. here, here, and here. But I don't quite get it yet.

I think that your misunderstanding comes from the fact that you assume that $f[n]$ is periodic. If $f[n]$ were periodic then its DTFT would be a Dirac comb. What is happening here is that you want to find the $N$-periodic sequence whose DFT coefficients equal equidistant samples of the DTFT of $f[n]$. And it turns out that that sequence is given by
$$\tilde{f}[n]=\sum_{k=-\infty}^{\infty}f[n-kN]\tag{1}$$
which is - as expected - an aliased version of the original sequence. Remember that sampling in one domain corresponds to aliasing in the other domain.
A derivation of $(1)$ is given in this answer.


Filtering out a specific Gaussian from a signal using Fourier concepts I have been wondering about a problem and need some help.
Let's say I have a 1-d time signal (discrete), I also have the knowledge that the signal consists of repeated Gaussians with known fixed standard deviation and for starters let's say of known a fixed amplitude. So the signal consists of series of similar Gaussians and then obviously there is the random noise.
Is there any way to use that knowledge in frequency filtering to get rid of the noise?

You can try performing a matched filtering with the gaussian of interest. See the below python code for an example. It is not exactly a work in the Foruer domain, (but can be seen as such) and remove thes noise just leaves the Gaussians. You need to adapt the sigma of your gaussians (i just did some rough trial-and-error here)

signal = np.loadtxt("C:/local/mmatthe/work/sheet1.csv")
t = np.arange(-signal.shape[0]/2, signal.shape[0]/2)
t2 = np.arange(-20, 20)
sigma2 = 10
gauss = lambda t: 1/np.sqrt(2*np.pi)*np.exp(-(t)**2/(2*sigma2))
gauss_samples = gauss(t2)
gauss_samples /= np.sqrt(sum(gauss_samples**2))

plt.figure(figsize=(10, 5))
plt.subplot(1,2,1)
plt.plot(t, signal)
plt.plot(t, gauss(t))

plt.subplot(1,2,2)
plt.plot(np.convolve(gauss_samples, signal));


ARMA (Auto Regressive Mean Average) Process Representation as AR (Auto Regressive) Let's say we have an ARMA (Auto Regressive Moving Average Model) process where the transfer function is a minimum phase system (Namely Invertible).

By Wold's Decomposition it is guaranteed to have MA representation of that system (Infinite order).

My question is could it be also represented by AR (Auto Regressive) system (Finite or infinite order)? If so, how could the coefficients be calculated?

Ok,
The answer is as following, any ARMA process which is Minimum Phase system could be represented as MA or AR system.

The MA part is easy by Wold's Decomposition.
Since the filter in Wold's Decomposition is Invertible (Minimum Phase by itself), Shall be called $ H \left( z \right) $, It is easy to build the optimal linear predictor from it, $ P \left( z \right) = 1 - \frac{1}{H \left( z \right)} $, which a linear combination of the past of the process -> AR Model.

It is all coming from the property of the Minimum Phase system being analytic in the unit circle.


Why remove energy at Nyquist before ifft? I've heard that you should remove the energy at Nyquist before performing the ifft.  Why is it that you can't leave the Nyquist energy there, is there danger of aliasing somehow?

you don't have to set $X\left(\frac{N}{2} \right)=0$ if you don't want to.  it will correspond to this component:

$$ X(k)\frac{1}{N}e^{j 2 \pi \frac{nk}{N}}\bigg|_{k=\frac{N}{2}} = X\left(\frac{N}{2} \right)\frac{1}{N}(-1)^n $$

but when you sample some $x[n]$, FFT it and find that $X\left(\frac{N}{2} \right) \ne 0$, you do not know the phase of that Nyquist component before it was sampled.  the Nyquist component actually aliases with itself in such a way that the quadrature part of the negative frequency cancels the quadrature part of the positive Nyquist frequency component.  such would look like:

$$ = \frac{X\left(\frac{N}{2} \right)}{N}(-1)^n = \frac{X\left(\frac{N}{2} \right)}{N}\left( \cos(\pi n) + A \sin(\pi n) \right) $$

$A$ could be anything and you wouldn't know it because $\sin(\pi n)$ is always zero for any integer $n$.

so if you're doing audio, i would likely set both the Nyquist and the DC components to zero and also insure the complex-conjugate symmetry with the other bins before iFFT.  just to make sure that, within rounding error, the imaginary part of your result is zero (the result is purely real) and that you have an unambiguous magnitude and phase for every non-zero component.


PWM samplerate to equal PCM accuracy How fast samplerate 1 bit PWM needs to equal 16 bit 48khz PCM?

1 second signal = 48000 samples,each sample can be anywhere from -32767 to +32768,therefore 1 second of 16bit 48khz PCM contains 65536^48000 possible combinations,thats very big number!

I was thinking about this,correct me if I am wrong,for 1 bit PWM signal to equal acuracy of 16bit 48khz PCM,it would need to run at 3145728000hz samplerate,that is 3 GHZ.

Now I would like to mention music format known as DSD,its 1 bit,running at 64x sampling rate of RedBook,that is 44100 x 64 = 2822400hz... thats 3 MHZ not GHZ,that is orders of magnitude bellow what I calculated to be minimum for 1bit signal to equal acuracy of PCM.

DSD is using noiseshaping,it helps to push the noisefloor down in the audible range,but that doesnt make up for fact that the samplerate is far too low to be able to reproduce all the possible combinations of sound in given lenght of signal equal to PCM ( right? )

PWM in audio is widely used in Class D amplifiers,doing a research I found out that they run even slower,some of them bellow the MHZ range,250khz for example,I am thinking,PWM freqency that slow must have much lower accuracy than PCM

It must be throwing away information because the 1bit PWM even when running at high speeds just doesn't equal the data bandwidth of PCM if it's not in the gigahertz range. Is what I have written true or false?

Short answer: For a direct implementation of PWM in discrete-time, yes, for 48 kHz, 16-bit audio, a clock of 3.2 GHz is required.

Direct Stream Digital (DSD) format uses 1-bit $\Delta$-$\Sigma$ modulation (DSM), which is somewhat different than pulse-width modulation (PWM), although both generate an on-off switched sequence as the output.

Regular PWM uses a comparator (the sign function), where the input is a modulating (the desired signal) and a high-frequency carrier signal. The harmonic and intermodulation distortion between the desired signal and the carrier signal are what is going to determine the unwanted signal output limiting the signal-to-noise-and-distortion (SINAD) ratio. With a triangle-wave carrier with frequency $\omega_c$, the output for a single sinusoid with frequency $\omega_v$ and amplitude $M$ between 0 and 1 is:
$$
PWM(\omega) = \frac{1}{2} M \cos{\omega_v t} \\
+ 2 \sum_{m=1}^{\infty} \frac{J_{0}\left(\frac{1}{2} m \pi M \right)}{\pi m} \sin\left(\frac{1}{2} \pi m \right) \cos\left(m \omega_c t \right) \\
+ 2 \sum_{m=1}^{\infty} \sum_{n=\pm 1}^{\pm \infty} \frac{J_{n}\left(\frac{1}{2} \pi m M \right)}{\pi m} \sin\left(\frac{1}{2} \pi (m+n) \right) \cos\left((m \omega_c + n \omega_{v} ) t\right)
$$
$J_n$ denotes the n-th order Bessel function.
The first line has your desired sinusoidal signal, the second and third line are harmonic and intermodulation distortion. This is entirely deterministic and similar spectra can be computed for any input signals that can be described by Fourier-series. Sample rate does not play a role in the ideal, continuous-time case. The distortion products can typically be moved to higher frequencies and filtered out by a low-pass filter by choosing a higher carrier frequency $\omega_c$, but there is no guarantee that intermodulation products will not appear in the baseband. The third line contains what is called foldback distortion, and as you can see, there will always be distortion in the baseband, and there is no way to guarantee that some input signal will not cause some combination of superpositions that cause severe foldback distortion. The best you can hope for is that the foldback distortion is small for most signals (which it usually is, and it is dominated by other noise sources).

Discretization (quantization and sampling) effects cause additional distortion, and several PWM methods try to take this into account. As you have discovered, you need very high time-resolution in order to represent large word-widths where there is discretization in time. This does not mean that the (average) switching frequency has to be high, only that the time-resolution is high. That is, you need to represent your output sequence as a non-uniformly sampled signal. Hence you don't need a high sampling rate, only a mechanism on the output that can cause non-uniform switching with that accuracy. This does makes it difficult to accurately simulate PWM. There are several PWM-like methods that seek to solve these distortion problems.

DSM (noise-shaping using 1-bit) also uses the comparator, but there is no carrier signal. Instead there is feedback from the error between the input and output of the comparator. By virtue of Bode's sensitivity integral, it is then possible to spectrally shape the error, dependent on the feedback filter used. As the comparator is a discontinuous non-linearity (no Lipshitz constant), there does currently not exist any theory with regards to the stability of DSM, and the successful design is entirely dependent on simulations. It works most of the time because in most cases you have an actuating signal to the comparator that causes some sort of averaging, causing the discontinuity to become effectively continuous. The DSM does go unstable and produce weird results in some pathological cases. So in general, it does not work, it just works most of the time. In DSM the distortion is also entirely deterministic, although it can be chaotic. It is impossible to de-correlate the distortion in DSM from the input.

I do not know of any methodology other than (Monte-Carlo) simulation for reliably determining an approximate SINAD in either case.

I've added a couple of figures comparing PWM and DSM with 1-bit truncation.

Time-series after reconstruction:

Power Spectral Densities (up to 2 kHz):

For DSM a 5th order feedback-filter was used, oversampling rate of 128, and sample rate of 1 MHz. With this it is possible to achieve an ENOB in excess of 17 bit with a 2-kHz baseband. For audio this would have to be 20 kHz, hence a sampling rate of 10 MHz. The PWM carrier was at 4 kHz. Most of the distortion (the noise-floor) stem from time discretization, not having introduced proper pre-distortion needed for discrete-time PWM, and will not represent ideal (natural) PWM very well. I will see if I can update this later. Using proper pre-distortion you should see a SINAD in excess of 100 dBc (ENOB > 16 bit). A 1-bit truncated signal can have an ENOB less than one since the ENOB is based on the uniform quantization error distribution assumption with is only approximately valid for 7-bit quantization or higher.

Why not just use a 1-bit DSM re-quantization (e.g. DSD) which is much more linear than a PCM-PWM conversion? The answer to this is switching losses. With high OSR and high and irregular switching activity, the DSM bit-stream in its basic form is not very suitable to drive high power Class-D amplifiers. Likewise, PWM due to its two-level representation and high jitter susceptibility (the need for good switching time resolution) is not very suitable for high resolution small signal DACs.

It might also be noted, that it is quite common for Class-D amplifiers to re-quantize to e.g. 8-bit using DSM before producing the PWM signal, instead of using complicated pre-distortion and non-uniform switching techniques, as a low word length reduces the required timing-accuracy significantly.

Some literature to check out:

Direct Digital Pulse Width Modulation for Class D Amplifiers

Signal processing for high resolution pulse width modulation based digital-to-analogue conversion

Digital-to-Analog Conversion in High Resolution Audio

Understanding Delta-Sigma Data Converters

Oversampling Delta-Sigma Data Converters

Lee's rule extended (on DSM stability and design)


please give the reason why every notation for DFT is valid? $X$ represents sample in frequency domain and $x$ represents samples in time domain.

NOTATION 1

$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \ e^{-j \frac{2\pi}{N} n k} $

$ x[n] = \frac{1}{N} \sum\limits_{k=0}^{N-1} X[k] \ e^{j \frac{2\pi}{N} n k} $

NOTATION 2

$ X[k] = \frac{1}{N} \sum\limits_{n=0}^{N-1} x[n] \ e^{-j \frac{2\pi}{N} n k} $

$ x[n] = \sum\limits_{k=0}^{N-1} X[k] \ e^{j \frac{2\pi}{N} n k} $

NOTATION 3

$ X[k] = \frac{1}{\sqrt{N}} \sum\limits_{n=0}^{N-1} x[n] \ e^{-j \frac{2\pi}{N} n k} $

$ x[n] = \frac{1}{\sqrt{N}} \sum\limits_{k=0}^{N-1} X[k] \ e^{j \frac{2\pi}{N} n k} $

NOTATION 4

$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \ e^{j \frac{2\pi}{N} n k} $

$ x[n] = \frac{1}{N}\sum\limits_{k=0}^{N-1} X[k] \ e^{-j \frac{2\pi}{N} n k} $

Please observe the scaling factor $\frac{1}{N}$ and change of a negative sign $-$ over the exponent term ${e^{\pm j\frac{{2\pi }}{N}nk} }$. Why is every notation of DFT valid?

It simply does not matter and it is just a matter of convention which version is used. I always use the first definition, which is also the one you will find in (almost) all textbooks on DSP. If you take always the same time domain signal $x[n]$ then it is obvious that $X_2[k]$ (according to notation 2) and $X_3[k]$ (according to notation 3) are simply scaled versions of $X[k]$ (according to notation 1). The one with scaling $1/\sqrt{N}$ is the unitary version of the DFT (i.e. the DFT becomes a unitary transformation). It's just notation number 4 that is a bit weird. In this case, $X_4[k]=X[-k]$, but I have never seen the sign convention like this. In this case, $x[n]$ can be interpreted as the Fourier series coefficients of $X[k]$. But again, all four versions are OK and you're free to use any of them as long as everybody knows which version of the DFT you mean.


Are there equations for placing sounds in the stereo field using delays? I'm thinking of the typical stereo widener (haas, binaural) effects here, those that extend from the typical pan knob to create "out of the headphones" -experience.

These processors can be automated to move sounds around, but I was wondering whether there was mathematics that described how to set the delays in order to move the sound into a certain spot in an angle -like fashion.

Such as in this plug-in:

What about modelling the "distance" from the listener?

Spherical head, distant source

For a spherical head with the ears at opposite points on its surface, a planar wave from a distant sound source reaches an ear on a straight line through air or if blocked by the head then partially through the shortest path along the surface of the sphere. If the sound source is at the same vertical coordinate $z$ as the ears, we can think of the head as a circle on the $x, y$ plane:

For simplicity, the head is a unit circle centered at $x = 0, y = 1$, turned so that the nose moves by an angle $\alpha$ counterclockwise from the origin. We are looking for the distances the sound still has to travel after the planar wavefront first reaches the head. There is a shortest straight line path (red dashed) of length $1-\sin(\alpha)$ to the right (red) ear. The shortest path (white dashed) to the left (white) ear consists of a straight line segment of length $1$ plus a curved path of length $\alpha$. The difference between the lengths to the two ears is $\alpha + 1 - (1 - \sin(\alpha)) = \sin(\alpha) + \alpha$. This holds for $0 \le \alpha \le \frac{\pi}{2}$, but because the function is antisymmetrical, it applies also to $-\frac{\pi}{2} \le \alpha \le 0$ giving a negative value for that range.

Modifying the equation to include the speed of sound $c$ and the head radius $r$, the inter-aural time difference (ITD) turns out as:

$$\text{ITD}(\alpha) = \frac{r}{c} (\sin(\alpha) + \alpha)$$

The formula is applicable for angles $-\frac{\pi}{2} \le \alpha \le \frac{\pi}{2}$, where $\alpha = 0$ means the source is straight ahead and $\alpha = \frac{\pi}{2}$ means it is 90° to the right, if the left channel is delayed for positive ITD. The plotted $\sin(\alpha) + \alpha$ part swings from $-\frac{\pi}{2} - 1$ to $\frac{\pi}{2} + 1$, so for a realistic maximum ITD of 0.6 ms (for distant sources) one can use simply:

$$\text{ITD}(\alpha) = \frac{0.6\ \text{ms}}{\frac{\pi}{2} + 1} (\sin(\alpha) + \alpha) \approx 0.23\ \text{ms}\ (\sin(\alpha) + \alpha)$$

For sources that are behind, before using the formula, add or subtract $\pi$ to/from $\alpha$ to bring it to $-\frac{\pi}{2}\ldots\frac{\pi}{2}$ and flip its sign.

Spherical head, nearby source

If one wants to take into account the distance to the source, the wavefront becomes spherical, or circular on the $x, y$ plane. Let's look at a head with just one ear considered:

The head is again a unit circle, but this time centered at the origin. The source (blue) is on the negative $y$ axis at a distance $d$ from the center of the head. The ear (red) is rotated by an angle $\alpha$ counterclockwise from the positive $x$ axis, or clockwise if $\alpha$ is negative. We shall constrain $-\frac{\pi}{2} \le \alpha \le \frac{\pi}{2}$ so that the shortest air path to the ear is always via the right side of the head as pictured. At the time the sound first reaches the head, the wavefront (blue) is cicular with radius $d-1$. We don't need that for the calculations, because we can use the time the sound leaves the source as the starting point. There is some negative angle $\beta$ such that if $\alpha \le \beta$, there is a direct air path of length $\sqrt{\cos(\alpha)^2 + \left(d + \sin(\alpha)\right)^2}$ to the ear located at $x = \cos(\alpha)$, $y = \sin(\alpha)$. Otherwise the sound travels a distance $D = \sqrt{\cos(\beta)^2 + \left(d + \sin(\beta)\right)^2}$ to the point $x = \cos(\beta)$, $y = \sin(\beta)$ where the straight path and the circle meet tangentially, and transitions to a further curved path of length $\alpha-\beta$. There is a right triangle (pink) that enables to solve $\beta$ using the pythagorean theorem:

$$d^2 = D^2 + 1^2\\
\Rightarrow \beta = - \text{asin}(\frac{1}{d}),\ D = \sqrt{d^2-1}$$

To recap, modifying the equations to include head radius $r$ and speed of sound $c$:

$$-\frac{\pi}{2} \le \alpha \le \frac{\pi}{2},\\
T_l = \left\{\begin{array}{ll}\frac{r}{c}\sqrt{\cos(\alpha)^2 + \left(\frac{d}{r} + \sin(\alpha)\right)^2}&\text{if } \alpha \le - \text{asin}(\frac{r}{d}),\\
\frac{r}{c}\left(\sqrt{(\frac{d}{r})^2-1} + \alpha + \text{asin}(\frac{r}{d})\right)&\text{otherwise,}\end{array}\right.\\
T_r = \left\{\begin{array}{ll}\frac{r}{c}\sqrt{\cos(\alpha)^2 + \left(\frac{d}{r} - \sin(\alpha)\right)^2}&\text{if } \alpha \geq \text{asin}(\frac{r}{d}),\\
\frac{r}{c}\left(\sqrt{(\frac{d}{r})^2-1} - \alpha + \text{asin}(\frac{r}{d})\right)&\text{otherwise,}\end{array}\right.\\
\text{ITD}(\alpha) = T_l(\alpha) - T_r(\alpha)$$

Again, for sources that are behind, before using the formula, add or subtract $\pi$ to/from $\alpha$ to bring it to $-\frac{\pi}{2}\ldots\frac{\pi}{2}$ and flip its sign.

For increasing distances, this definition of ITD approaches the previous one:


Is there a need of Point interpolation before proceeding for gaussian smoothing of an incomplete distribution? Suppose there is a distribution that has values sampled on the interval 1-25 with corresponding sample values that have to be smoothed. For example:

n    f[n]

3    123
5    234
8    246
10   567
14   253
18   251
24   132

The first column has sample points between 1-25 but does not have samples for all points on the interval. Should interpolation of other points in the distribution be done so that we have approximate values for all points between 1-25 before applying Gaussian smoothing?

Irregularly sampled data need not necessarily be interpolated prior to smoothing. You can try to smooth your data in the following way:

$$y(t_n)=\frac{\sum_{k}g(t_n-t_k)x(t_k)}{\sum_{k}g(t_n-t_k)}\tag{1}$$

where $t_n$ are your sample times, $x(t_n)$ are your data, and $g(n)$ is your smoothing kernel. Note that $t_n$ need to be integers (as in your example), which is no restriction, just a matter of formulation. Since $t_n$ are irregular, you do not always use all values of the chosen kernel $g(n)$, but this is compensated for by the normalization in (1). The kernel width has to be chosen such that the kernel always covers at least several of your irregularly sampled data. The smoothing in (1) is also used in the Nadaraya-Watson estimator.


Phase versus time graph of QPSK In wikipedia, there was an example on QPSK modulation where the main signal was 11000110 and I contained 1001 and Q contained 1010, now their output was shown by the Signal in the figure.

Now according to constellation diagram Signal can be represented as

11 -> pi/4
00 -> -(3*pi)/4
01 -> (3*pi)/4
10 -> -pi/4

But my confusion was in the phase vs time diagram shown in another figure given below

For 11 it starts from phase pi/4, but for 00 shouldn't the phase start from -(3*pi)/4 = -135 ?

You're right that the phase plot doesn't match the signal plot. Note that it also shows more symbols than the signal plot. Those two plots do not refer to the same bit sequence, and the phase plot is only meant to illustrate the difference between QPSK and OQPSK (the original figure consists of two plots, one for QPSK and one for OQPSK).


FFT, interpolation and phases I am trying to understand how to interpolate discrete sinusoidal data that I (for now) generate myself, and I get problems that I don't understand. I'm a newbie in this field so do not hesitate to point out even the most basic errors.

I have tried two approaches.

The first is to take (Python's/numpy's) fft(), get amplitudes with abs() and phases with
$$
\frac{\text{angle} - \pi/2}{2 \pi} * \text{signal_length}
$$

Then, for each sampling point $p$ I do a sine series,
$$
\sum_i \text{amplitude}_i * \sin(2 \pi i * \frac{p - \text{phase}_i}{\text{signal_length}})
$$

For artificial data where the (only) sine lambda is signal_length the fit is perfect regardless of the artificial data phase. However, if I generate the data with sine that has lambda of half the signal length, the fit is far from good for most phase offsets. Obviously, the fit is also not good when I generate a signal from several added sines of different lambdas.

Edit 1: Note that I can artificially correct the phases of different lambdas by dividing them by n, where n is signal_length / lambda of the generating function, but this doesn't help when I don't know the lambda of the generating function and/or when the signal is a mix of several sines.

The second approach was to get separate sine and cosine amplitudes from fft() output by real() and -1 * imag() and then do a summation similar to above for both sine and cosine terms. This works well whenever the phase offset of the generated data is n*2*pi but unwell when it's anything else.

Eventually I'd like to be able to interpolate data generated by a superposition of several sinusoidal functions of various lambdas and phases, but so far it's not looking good. I think there is something that I don't understand about Fourier/DFT.

Note that I can't/don't want to use any built-in ifft functions but I need to do some kind of a sine/cosine summation manually. I tried to write the definition of inverse DFT, but the built-in functions don't understand correctly things such as e^i. Ideally I would want to sum just the sine or the cosine terms, but not both.

Edit 2 LaTeXed some math

I resolved this problem by myself. Using approach 1, I replaced $$phase_i$$ with $$\frac{phase_i}{i+1}$$ in the summation formula, and suddenly everything worked perfectly.


Low pass vs band pass filter I've been reading some EEG papers where the authors use a bandpass filter at, say, 1Hz and 30Hz. I've also seen some accelerometer papers where they bandpass at 0.5Hz and 10Hz

What is the difference between a bandpass (1Hz, 30Hz) and just a regular low pass filter at 30Hz?

It seems like they're trying to remove really low frequencies, but does it really make that much of a difference? Is there a particular reason why this is done?

A bandpass filter allows signals between two specific frequencies to pass, but that discriminates against signals at certain frequencies.

A lopass filter is a filter that passes signals with a frequency lower than a certain cutoff frequency and minimizes signals with frequencies higher than the cutoff.

There's a slight distinction. What's useful depends on the situation.


Why do we need deterministic measurement matrices in compressed sensing? I recently introduced myself into the field of CS, but I do not understand why some people try to find deterministic measurement matrices?
If I am correct, gaussian random matrices are very powerful with their properties. They are easy to create, their reconstruction power is very good (having a low RIP constant) and so on.

As far as I know there are two reasons:

In sensing part: For practical implementation, usage of random matrices
is hard, so people try to come up with simpler matrices that are fixed,
this is thought to be helpful for designing efficient hardware.
In reconstruction part: To find deterministic matrices which their
structure can be exploited for faster signal recovery. Using a
structured matrix, might pave the way for efficient recovery. The
goal is to achieving something like FFT for compressive recovery.


Are MATLAB function zp2tf() and tf2zp() are complementary or not? I was under impression that given, pole, zero and gain the transfer function (filter coefficients b and a) is fixed. Therefore, if we designing a filter using cheby2() for some coefficients to get the pole zero and gain by using

[z, p,k ] = cheby2( N, Rs, Ws, 'stop');

I assumed that we can get back the [z,p] by

[hB,hA] = zp2tf(z,p,k);

[z1,p1,Hek] = tf2zp(hB,hA);

In other words z1=z and p1=p. But this is not true as shown by my MATLAB code.
Can anyone please tell me what am I missing?

This happens frequently if your poles are reasonably close to the unit circle. Consider the following example

%% TF2ZP is problematic
fs = 44100;
% 6th order lowpass, fc = 50Hz, sampled at 44.1kHz
[z,p,k] = cheby2(6,80,50*2/fs);
% to transfer function
[b,a] = zp2tf(z,p,k);
% back to zpk
[z1,p1,k1] = tf2zp(b,a);
display([p p1]);

Displaying the poles side by side, you get

         Correct Poles                  TF2ZPK
      0.99729 + 0.00078062i       1.0007 +          0i
      0.99729 - 0.00078062i      0.99983 +  0.0033549i
      0.99814 +  0.0019985i      0.99983 -  0.0033549i
      0.99814 -  0.0019985i      0.99691 +  0.0030633i
      0.99936 +  0.0025696i      0.99691 -  0.0030633i
      0.99936 -  0.0025696i      0.99535 +          0i

The poles are very different and one pole has actually moved out of the unit circle and the resulting filter is unstable.

This is caused by limited numerical precision. $tf2zp()$ requires calculating the roots of a polynomial which is a numerically tricky problem especially if the roots are very close together. $zp2tf()$ works perfectly fine but  $tf2zp()$ doesn't.

Morale of the story: If you have poles close to the unit circle, keep your filters in ZPK or SOS representation and avoid transfer functions. Avoid $tf2zp()$ if possible.

  Are MATLAB function zp2tf() and tf2zp() are complementary or not?

Sort of. They are intended to be complimentary but like any complimentary pair this is subject to numerical noise of the implementation. In this particular instance there are reasonable use cases, where the numerical noise can be substantial to the point where the function becomes useless.


Dealing with echoes / dead-spots in digital transmissions I'm looking at communicating digital data through the air as soundwaves.

To illustrate the problem, I generate a 10-cycle 2kHz wave, play it out through my MacBook's speakers and pick it up on my iPod Touch:

The echoes are about half the amplitude of the probing signal.

HotPaw pointed out to me in a comment on a recent question that there is a potential problem with dead-spots: regions of destructive interference.

If I use a particular carrier frequency, transmitting from one point in the room to another, there will be certain locations where acoustic cancellation / destructive interference occurs, resulting in a dead spot.

Surely other technologies such as Wi-Fi must encounter the same difficulties. How do engineers deal with this phenomenon?

The only solution I can think of is to simultaneously use say three different carriers, and have the listener select the strongest.

PS severely hampered on the tags: data transmission, communication, radio, information transfer, echo, Radio, rf, signal, bleh

If the channel is static and you are in a deep fade -- then, yes, your only option is to increase the bandwidth or go multicarrier like you suggest.

In fact, modern WiFi systems use that exact technique: OFDM ("Orthogonal frequency division multiplexing" -- sometimes known as multi-tone or multicarrier transmission).

Alternatively, you could go with two receive antennas (microphones in your case) spaced far enough. If you use 2000Hz carrier, the wavelength is v/f = 300[m/s] / 2000 Hz = 0.15 m. So a distance of around 10cm should give you enough diversity. This can also be implemented at the transmitter.


I want to find the first period $T_0$ in a signal $x[n]$ and save it to a vector? I want to find the first period $T_0$ in a signal $x[n]$ and save it to a vector $\hat x$?  $x[n]$ is read in by the command:

[x, fs] = audioread('Data/sample1.wav');

First of all, I have calculated the main frequency of the signal x[n]. I did it by

[x, fs] = audioread('Data/sample1.wav');
N = length(x);
f = (0:(N-1))*fs/N;
xdft = fft(x);
[~,index] = max(abs(xdft(1:length(x)/2+1)));
freq = 0:(fs/length(x)):fs/2;

Ok. freq is now my main frequency. As a next step, I want to calculate T_0.

T_0 = 1/fs;

And now I have to find this sequence, but how to do this? If I guess the signal in the wave is sines or cosines, but I do not know how to check this.

The time period should be inverse of the signal's frequency

  T_0 = 1/freq

not sampling frequency


Bandwidth computation for NRZ baseband transmission scheme In Andrew Tanenbaum's book on Computer Networks, while explaining about bandwidth of a NRZ (Non Return to Zero) baseband transmission scheme, he says the following:

  With NRZ, the signal may cycle between the positive and negative levels up
  to every 2 bits (in the case of alternating 1s and 0s). This means that we need a bandwidth of at least B/2 Hz when the bit rate is B bits/sec. This relation comes from the Nyquist rate [Eq. (2-2)].

I have 2 questions here.

What does the first statement mean exactly? Does the author imply that for a signal that switches between 1s and 0s alternately, if a particular bit is a 1, then the next time we get a 1 again is after 2 bits?
How does the first statement suggest the second statement?

The following excerpt from the book explains Nyquist's theorem and is followed by Eqn. (2-2).

  Nyquist proved that if an arbitrary signal has been run through a low-pass filter of bandwidth B, the filtered signal can be completely reconstructed by making
  only 2B (exact) samples per second. Sampling the line faster than 2B times per
  second is pointless because the higher-frequency components that such sampling
  could recover have already been filtered out. If the signal consists of V discrete
  levels, Nyquist’s theorem states:

$$maximum\; data\; rate = 2B \log_2 V \;bits/sec$$

The first statement says that for the given binary sequence of $$a_k = \{ 1,0,1,0,1,0,...\}$$ the transmitted NRZ line-signal would oscillate at the maximum rate :
$$V = \{+A,-A,+A,-A,...\}.$$ which is alternating at a rate half of the binary data transmission bit-rate. This maximum frequency of line voltage oscillation should be capable of passing through the channel filter. Hence, the channel bandwidth should be at least 500 Hz in that case.

For example, if you want to transmit binary data at 1000 bits per second, by using NRZ, then the line pulse voltage alternates maximum at 500 Hz fundamental frequency.

And the second statement then follows the first: for such an alternating signal, your channel bandwidth should be at least 500 Hz (half the binary rate in Hz.) to a first approximation ideally. You should consider noise, pulse shaping and ISI issues for determining the required bandwidth more realistically.

Nyquist sampling theorem states that for a baseband signal of bandwidth, say, 4000 Hz, you need to take at least 8000 samples per second to exactly reconstruct the continuous signal back. And if you use 8 bits per sample, then that requires 8 x 8000 = 64 kbps binary data transmission rate. Note that such a binary rate would require a channel bandwidth of at least 32 kHz if NRZ is used.


autocorrelation of multiple signals Problem: I am looking at an adaptive filtering application where the eigenvaluespread of the autocorrelation matrix $R$ is important for the convergence of the algorithm. For a single channel system the autocorrelation matrix $R$ for iterationstep $n$ can be calculated by $R=E\{ x(n)  x^H(n)\}$ where $x(n)$ is the input signal of the adaptive filter at iterationstep $n$ consisting of a number of samples $N$ recorded over a timespan. The calculation of the eigenvalues is straight forward.
Question: What is the "multichannel equivalent" for $R$ in the case of e.g. an adaptive filtering multichannel application? Do I need to calculate some sort of autocorrelation tensor?

You apply the same formula, but instead of using a scalar $x(n)$, you will have $x(n) \in \mathbb{C}^{M \times 1}$, where $M$ is the number of channels, and $x^{H}(n)$ is the $1 \times M$ matrix whose entries are the complex conjugate of the entries in $x(n)$. As a result $R \in \mathbb{C}^{M \times M}$ is an Hermitian matrix, and its eigenvalues are all real.


Term for integral multiples of fundamental period? As we know that a signal has fundamental frequency and integral multiples of fundamental frequency are known as harmonics

I am curious about the term that will be used for integral multiples of fundamental time period?What is that term?

If $T$ is a positive real number such that
$$x(t) = x(t+T) ~\text{for all}~t \in \mathbb R\tag{1}$$
for continuous-time signals, or if $T$ is a positive integer and
$$x[n] = x[n+T]~\text{for all}~n \in \mathbb Z \tag{2}$$ for discrete-time signals, then $x$ is said to be periodic with period $T$ or to have period $T$.
A signal $x$ that has period $T$ also has period $2T$, period $3T$, $\ldots$, but might (or might not) have smaller period such as $T/2$ or $T/3$ etc. The smallest value of $T$ for which $(1)$ or $(2)$ is called the fundamental period of the signal, and all other periods necessarily are integer multiples of the fundamental period.

So, to answer your question about the name for integer multiples of the fundamental period, they can be called a period of $x$ which is perfectly valid, or, if it is necessary to emphasize the fact that $T$ is not the fundamental period, one can say more prolixly, "$T$ is a period but not the fundamental period of $x$." The more one want to nitpick (or strive for clarity and completeness), the more wordy one must be.

Note that if $T$ is the fundamental period of a signal $x$ and if $T^\prime = mT$ (where $m$ is an integer greater than $1$) is a period, then we can express the continuous-time signal $x(t)$ of period  $T^\prime$ as a Fourier series with fundamental frequency $$f_0^\prime = \frac{1}{T^\prime} = \frac{1}{mT}= \frac{1}{m}f_0$$ where $f_0 = \frac 1T$ is the fundamental frequency of $x(t)$ of fundamental period $T$. $f_0^\prime$ is called a subharmonic of $f_0$.  Notice, however that the Fourier series that we get by regarding $x(t)$ as a signal of period $T^\prime$ is exactly the same as the Fourier series that we get by regarding $x(t)$ as a signal of period $T$. In
$$x(t) = \sum_{n=-\infty}^\infty c_n \exp(j2\pi nf_0^\prime t)
= \sum_{n=-\infty}^\infty c_n \exp(j2\pi \frac nm f_0t)\tag{3}$$
the $c_n$'s have value $0$ except when $n$ is an integer multiple $km$ of $m$ and so the right side of $(3)$ can be expressed as
$$x(t) = \sum_{k=-\infty}^\infty c_{km} \exp(j2\pi k f_0t)$$
which is precisely the Fourier series for $x(t)$ regarded as a signal of period $T$. In short, while we have introduced the concept of subharmonics of $f_0$, $x(t)$ has no frequency content at these subharmonic frequencies.

A similar result holds for discrete-time signals $x[n]$ but I am too lazy to write it all out. Perhaps @OlliNiemitalo will resuscitate his answer and spell it all out for those the benefit of those who are interested.


removing vertical lines in grayscale images by thickness Suppose we have a grayscale image that contains vertical lines. Now suppose that not all vertical lines are the same, some of them have different thickness.

Question is, is there a way, in MATLAB or other programming language, to filter the vertical lines by thickness? In other words, based on the user selected thickness, only thin vertical lines are removed, while all other thick vertical lines are kept intact? or all thick vertical lines removed, all other thin vertical lines intact?

I tried experimenting with MATLAB with the morphological operators, but nothing seems to really work, perhaps I am using the wrong method.

edit:
I am thinking about one possible way to solve this, but not in code, just as an idea. Every vertical line has 2 vertical edges. For every pair of 2 edges (1 line), calculate distance between two edges. Then, based on user threshold, remove or keep corresponding lines based on these calculated distances.

edit2: here is sample image:

How to remove only 3 thin lines and keep 2 thick lines OR how to remove 2 thick lines and keep 3 thin lines? That is what I have to accomplish.

edit3: what I mean by remove is in this context interpolate them with the background (average value of left and right close pixels).

Here is some Python code that does mostly what you want.  You need to add the thick/thin line selection logic.  You'll probably want to tweak it some more too.  There seem to be some artifacts when it is turned back into a .jpg.

from PIL    import Image

import numpy as np

#============================================================
def main():

#---- Load the Image

        theImage = Image.open( "test.jpg")

        thePixels = theImage.load()

        theRangeMin = 0
        theRangeMax = 7

#---- Get the Size

        theSize = theImage.size

        theWidth  = theSize[0]
        theHeight = theSize[1]

        print theWidth, theHeight

#---- Scan some Rows

        theSum = np.zeros( theWidth )
        theScanCount = 0

        for n in range( 100, theHeight, 100 ):
            print n
            theScanCount += 1

            for m in range( theWidth ):
                theShade = thePixels[m,n]
                theSum[m] += theShade[0] + theShade[1] + theShade[2]

#---- Calculate Derivatives

        theFirst  = np.zeros( theWidth )
        theSecond = np.zeros( theWidth )
        theSum /= ( 3.0 * theScanCount )

        for m in range( 1, theWidth-1 ):
            theFirst[m]  = ( theSum[m+1] - theSum[m-1] ) * 0.5
            theSecond[m] = theSum[m+1] - 2.0 * theSum[m] + theSum[m-1]
            print m, theSum[m], theSecond[m], theFirst[m]

#---- Find Vertical Lines

        theThreshold = 10

        theLastFirst = 0

        theLines = []

        m = 1
        while m < theWidth-1:
            if theSum[m] < theThreshold:
                while m < theWidth-1:
                    if theSum[m] > theThreshold and theFirst[m] < theThreshold:
                        break
                    m += 1

                print theLastFirst, m
                theRange = m - theLastFirst
                theRate = 1.0 / float( theRange )
                if theRangeMin <= theRange and theRange <= theRangeMax :
                    theLines.append( [ theLastFirst, m, theRange, theRate ] )

            elif abs( theFirst[m] ) < theThreshold : theLastFirst = m

            m += 1

#---- Hide the Lines

        theLineCount = len( theLines )

        for n in range( theHeight ):
            for p in range( theLineCount ):
                theLine = theLines[p]

                theSpot = theLine[0]
                theFirstShade = thePixels[theSpot,n]
                theLastShade  = thePixels[theLine[1],n]
                theRange = theLine[2]
                theRate  = theLine[3]
                theFirstVal = ( theFirstShade[0] + theFirstShade[1] + theFirstShade[2] ) / 3.0
                theLastVal  = ( theLastShade[0]  + theLastShade[1]  + theLastShade[2] ) / 3.0
                theSlope = ( theLastVal - theFirstVal ) * theRate

                if n == 100: print theSpot, theFirstVal, theLastVal

                for r in range( 1, theRange ):
                    theVal = theFirstVal + r * theSlope
                    theGrey = int( theVal ) & 0xFF
                    theShade = ( theGrey << 16 ) + ( theGrey << 8 ) + theGrey
                    thePixels[theSpot+r,n] = theShade
                    if n == 100: print "     ", theVal

#---- Write the Image to a File

        theImage.save( "out.jpg" )

#============================================================
main()

Followup:

Here is the debug output at the first thin line:

194 119.857142857 8.85714285714 -1.14285714286
195 123.142857143 -3.71428571429 1.42857142857
196 122.714285714 -64.7142857143 -32.7857142857
197 57.5714285714 13.8095238095 -58.2380952381
198 6.2380952381 88.0 -7.33333333333
199 42.9047619048 42.7142857143 58.0238095238
200 122.285714286 -79.2380952381 39.7619047619
201 122.428571429 0.428571428571 0.357142857143
202 123.0 -0.142857142857 0.5
203 123.428571429 0.0 0.428571428571
204 123.857142857 0.142857142857 0.5
205 124.428571429 -0.285714285714 0.428571428571
206 124.714285714 1.42108547152e-14 0.285714285714


How to convert floating point code to fixed point? I am working on converting a floating point MATLAB code into fixed point C code. The code is to be run on a DSP platform, the specifics of which are yet to be decided.

I'm new to floating point to fixed point conversion and am facing quite a few challenges.

My questions are:

What is the best way for me to start and where should I start, because the topic is quite confusing.
My MATLAB code is basic filtering code; there's FFT, windowing followed by IFFT. Windowing seems simple enough as its multiplication, but I'm totally lost when it comes to FFT and IFFT.

Thanks!

Having a standard fixed point library available as a reference implementation would be helpful

There are more than a few fixed point dsp libraries, some with source code and documentation. Google returns many links.

TI and Analog Devices have fixed point libraries that support their products. The documentation and support is very helpful.

Matlab also offers fixed point add ons. Student pricing is very affordable.


Minimum Sampling Rate of Bandpass signal Consider a DSBSC signal s(t) = m(t).cos(2.pi.20k.t) where say m(t) is multi-tone signal with frequencies less than 5 kHz.

The spectrum gives a typical bandpass say triangle centred around 20k and -20k.
Now, when I calculate the minimum sampling rate;

Since the highest frequency is not an integer multiple of bandwidth of the bandpass signal; I end up getting the sampling rate as 25 KHz.

But when I draw the spectrum for the sampling rate of 10 kHz, the message signal centred at origin was doubled in amplitude and there was no aliasing as well.

Does that mean the minimum sampling rate is indeed 10 kHz?

I have used the typical bandpass sampling rate procedure many times, this was the first time it failed in giving me the minimum sampling rate!

edit: I doubt that the 10 kHz would be the minimum sampling rate to recover the message signal while the bandpass sampling rate of 25 kHz is the one which recovers the bandpass spectrum!

Also, if the above conclusion is correct, is it independent of the carrier frequency (Fc>>Fm)?

The problem is that the general criterion for bandpass sampling (which results in $25\text{ kHz}\leq f_s\leq 30\text{ kHz}$) applies for general signals, whose spectrum may not be symmetric around $f_c$. In your example, sampling at 10 kHz would fail because spectrums of different shapes overlap.

The most complete resource I know on bandpass sampling is "The theory of bandpass sampling", by Vaughan et al, IEEE Trans. on Signal Processing, 1991, Vol. 39, Issue 9 (http://dx.doi.org/10.1109%2F78.134430). At the start of page 4, it specifies that in the special case of DSB signals, the minimum sampling interval is $1/B$, where $B$ is the bandwidth (for both sidebands). This corresponds to your case, where the minimum sampling rate is 10 kHz, equal to the signal's bandwidth.


Confusion between samples per second and samples per symbol My understanding of sampling frequency is when we have an analog signal x(t) we convert it to a digital signal. Assume the sampling frequency is denoted by $f_s$ then the sampling time $T_s$. In my understading this is denoted the number of samples per second.

In MATLAB programs, I have come across generation of base band signal. The authors have upsampled the base band signal and the sampling rate was stated to be 3 samples per symbol.

My question is what does samples per symbol mean.Assuming we have QAM symbols are we repeating it three times?

Thanks

A digital communications signal such as QAM can be generally written as $$s(t) = \sum_{n=0}^{N}a_n p(t-nT),$$ where $a_n$ are (possibly complex) amplitudes taken from a finite alphabet, and the symbol rate is $1/T$. The pulses $p(t-nT)$ are often called "symbols". Matlab commands that ask you to specify "samples per symbol" are asking how many samples you want to use to define $p(t)$.

Of course, since you know the pulse duration, specifying the samples per symbol is an indirect way of specifying the sampling rate.


Finding Harmonics Analytically I've always obtained harmonics using simulation and never by hand.

My question is how can someone find the harmonics for a signal (say a square wave) by hand? What equations come in handy for finding the magnitudes at which each harmonic and sub-harmonic occur?

I want equations or techniques that are useful to find a closed form equation that gives us the magnitude of the harmonic and sub-harmonic for a given harmonic n or sub-harmonic $n\pm \alpha  \ \ \ \  \alpha \in  \mathbb{N}$

What you're asking for is:

  given a periodic signal (e.g. square wave), where are the discrete frequency components (i.e. harmonics)?

You'll find that the harmonics you're looking for are the coefficients of the Fourier series expansion of your signal.

So, if you have a signal that has a known Fourier series expansion, or can be constructed from such signals, it's easy.

But how to find the coefficients in general?

This is answered exactly by the Fourier transform.

So, the relevant formulas are, given a signal $s(t)$ with period $P$

the coefficients $c_n$ of the Fourier series expansion $s(t) = \sum\limits_n c_n e^{j2\pi n\frac tP}$: $$c_n=\frac{1}{P}\int_{t_0}\limits^{t_0+P} s(t)\cdot e^{-j \tfrac{2\pi nt}{P}}\ dt$$ describing the $n$th harmonic. You just plug in your square wave for $s(t)$ and get out the coefficients of your harmonics.

and in the general, you don't get discrete harmonics, but a continuous spectrum, and thus, instead of a discrete Fourier series the Fourier transform.

Not posting info on the Fourier Transform here; it's pretty much common knowledge, and you wouldn't benefit – your profile indicates that you're studying physics, and thus, you've already met the Fourier transform, for example when describing electrons in periodic crystals.


Hanning and Hamming Window I am new in the field of signal processing. What's the different between Hamming and Hanning window? When we use the former and latter?

Max is right that the difference between Hamming and Hann windows are small. (BTW, I am a proponent of the movement to totally do away with them term "Hanning".  There is no Dr. Hanning nor Mr./Ms. Hanning that the window is named after.)

The Hamming window is 92% Hann window and 8% rectangular window.  Hamming found out that he was able to reduce the height of the maximum side lobe by doing that.

Another issue that I just thought of:  The Hann window (and some others) is a complementary window; that is that the latter half of the window of one frame adds to the first half of the following frame to 1.  For analysis, the property of being complementary is usually not important, but for synthesis (or reconstruction), this property is salient.

So if you like symmetry about 0, the Hann window having even length $N$, is

$$ w\left( \tfrac{n}{N} \right) = \begin{cases}
\tfrac12 + \tfrac12 \cos\left( 2 \pi \tfrac{n}{N} \right) \qquad & |n| \le \frac{N}{2} \\
0 \qquad & |n| \ge \frac{N}{2} \\
\end{cases} $$

The property of being complementary is:

$$ \sum\limits_{m=-\infty}^{+\infty} w\left( \tfrac{n-mH}{N} \right) = 1 $$

then

$$\begin{align}
x[n] &= x[n] \times 1 \\
     &= x[n] \times \sum\limits_{m=-\infty}^{+\infty} w\left( \tfrac{n-mH}{N} \right) \\
     &= \sum\limits_{m=-\infty}^{+\infty} x[n] w\left( \tfrac{n-mH}{N} \right) \\
     &= \sum\limits_{m=-\infty}^{+\infty} x_m[n-mH] \\
\end{align}$$

where

$$\begin{align}
x_m[n-mH] & \triangleq x[n] w\left( \tfrac{n-mH}{N} \right) \\
x_m[n] & = x[n+mH] w\left( \tfrac{n}{N} \right) \\
\end{align}$$

So $x_m[n]$ is a finite-length frame of signal of length $N$ and represents the audio in the neighborhood of $x[n+mH]$ or $m$ hops from the beginning.  $H = \frac{N}{2}$ is the frame hop displacement.

The Short-Time Fourier Transform (STFT) is:
$$\begin{align}
 X_m[k] &= \sum\limits_{n=-\tfrac{N}{2}}^{\tfrac{N}{2}-1} x_m[n] e^{-j 2 \pi nk/N}
\\
\\
   &= \sum\limits_{n=-\tfrac{N}{2}}^{\tfrac{N}{2}-1} x[n+mH] w\left( \tfrac{n}{L} \right) e^{-j 2 \pi nk/N} \\
\\
&= \sum\limits_{n=0}^{N-1} \hat{x}_m[n] e^{-j 2 \pi nk/N} \\\end{align}$$

where

$$ \hat{x}_m[n] = \begin{cases}
x_m[n] \qquad  & 0 \le n < \tfrac{N}{2} \\
\\
x_m[n-N] \qquad & \tfrac{N}{2} \le n < N \\
\end{cases}$$


Why these cross-connections in the allpass filter? Why these cross-connections in the allpass filter?

and what do negative and positive $k_1$ really refer to?

Is it when the signal is going downwards or is it when the signal is below the zero line?

This is called a lattice structure implementation of all-pass filters.

$k_1$ and $-k_1$ are scalar gains. The input to the delay element is $$a[n]=x[n]-k_1a[n-1]$$ taking the $z$-transform, the transfer function is
$$A(z)(1+k_1z^{-1})=X(z)\Rightarrow A(z)=\frac{X(z)}{1+k_1z^{-1}}$$ and the output is
$$\begin{align}
y[n]&=k_1a[n]+a[n-1]\\
\end{align}$$
Thus
$$\begin{align}
Y(z)&=A(z)(k_1+z^{-1})\\
&=X(z)\frac{k_1+z^{-1}}{1+k_1z^{-1}}
\end{align}$$
$$H(z)=\frac{Y(z)}{X(z)}=\frac{k_1+z^{-1}}{1+k_1z^{-1}}$$
So it is an all-pass filter ($|H(e^{j\omega})|=1,\ \forall \omega$) with mirrored coefficients in numerator and denomerator.


Speech discrimination from animal sound i am looking for discrimination of speech from animal sound if you take a 150ms frame of speech and bird as shown below, what feature can i extract from them so that i can effectively state the difference

zero crossings ?
cross correlation of successive frames seems high in bird sound ? isn't it
any other feature from your side ?

one more example in above fig black-speeches,yellow siren,red bird,safron heavy wind the spikes are foot steps

above figure shows the difference in formant positions on spectrogram
in bird sounds the formants are differently spaced and less formants

will anything could be extracted from this analysis

the reason is i am looking for unsupervised(with out any training or SVM) speech discrimination form animal sounds because both are voice but computers cant classify the difference unless they are trained by any nueral network which need huge memroy and training database

The essence of your problem is one of Voice Activity Detection (VAD). This is a well researched topic and you should investigate the available algorithms - there are many papers available online, here for example

You should apply a VAD algorithm to the signal and replace non-speech areas with silence, you can then look at the specific edge cases which are concerning you, such as animal sounds have being incorrectly considered speech.

For the case of animal sounds you may want to investigate Formant Detection (for example this paper) since as you pointed out there is a distinction between them. You may also want to consider developing a heuristic based on your own empirical data about the differences. You can take many signal features such as the fundamental frequency, the signal power/amplitude, peakedness etc... and weight these appropriately to give you a confidence rating about the likelyhood of it being animal rather than human speech. This is a fuzzy approach, you would need to iterate a lot on the weightings and metrics.


Calculation of cosine in frequency domain instead of calculatin in time-domain followed by a FFT I got an $N$, in my case 512, point FFT of a real-valued signal. Based on some calculation in my application I determine the parameters $k \in [1, N-1]$, the number of oscillations per period, $\phi \in [0, 2\pi[$, phase. I use these parameters and the function
$\cos\left(2\pi\cdot \frac{k}{N}\cdot t + \phi\right)$
to describe a signal in time-domain. Instead of calculating all the points of the signal and then doing an FFT I would like to directly calculate the result in frequency domain.

I calculated the FFT coefficients for the function as following:

$$c_n = \frac{N}{2}\left(\mathrm{e}^{i\phi} \cdot \mathrm{sinc}\left(n-k\right)+\mathrm{e}^{-i\phi} \cdot \mathrm{sinc}\left(n+k\right)\right), \quad\text{ with }\mathrm{sinc}(x) = \frac{\sin(\pi x)}{\pi x}$$

Note: $N$ is in here because the output of the FFT algorithm I use is not scaled and scaling is done in IFFT:

I verified my coefficients using Mathematica with the FourierCoefficient function.

I calculated a 512 elements long vector of complex doubles in Matlab and set my first vector to zero as there is no (dc) offset.

Is this valid or do I have to calculate bin 0 too if my $k$ is not a whole number?
Can I use $c_0$ for it?

The next 255 values I filled $c_n$ with $n$ from -1 to -255. The next value I filled with zero as I was unsure which coefficient should go there. And then the next 255 values I filled with $c_n$ with $n$ from 255 to 1.

To verify the correctness of my generated signal I plotted the real value of the IFFT of my vector. The output was as expected as long as $\phi$ was 0 and $k$ was a whole number.

Plot (blue is what was expected and red the result from IFFT) with $\phi=0$ and $k = 8.3$:

Plot with $\phi=0.3$ and $k = 8$:

Plot with $\phi=0.3$ and $k = 8.3$:

For reference my matlab code is:

k = 8.3; %number of oscillations per period
phi = 0.3; %phase in radians

t = linspace(0,511,512);
n = linspace(-1,-255,255);
fft_coeff0 = 256*(exp(i*phi)*sinc(-k)+exp(-i*phi)*sinc(k));
fft_coeffs(1:255) = 256*(exp(i*phi)*sinc(n-k)+exp(-i*phi)*sinc(n+k));

reconst_fft = [fft_coeff0 fft_coeffs 0 conj(fft_coeffs(end:-1:1))];

y = cos(k/512*2*pi*t+phi);
plot(t, ifft(reconst_fft), 'r', t, y, 'b')

What could cause the difference between my expected output and the IFFT?
How can I fix it?

EDIT: Changed code to only calculate one side and use complex conjugating for the other side and added $c_0$ coefficient to vector.

Plot with $\phi=0.0$ and $k = 8.5$:

Adding $c_0$ fixed the wrong amplitude.

If you have a pure tone signal defined by:

$$ S[n] = M \cos( \alpha n + \phi ) $$

Where $f$ is the frequency in cycles per frame and $\alpha$ is in radians per sample:

$$ \alpha = f \cdot \frac{ 2\pi }{ N } $$

And use a $ 1/N $ normalized DFT:

$$ Z[k] = \frac{ 1 }{ N } \sum_{n=0}^{N-1} { S[n] e^{ -i\beta_k n } } $$

Where:

$$ \beta_k = k \cdot \frac{ 2\pi }{ N } $$

The most concise formula for calculating the DFT bin values is:

$$  Z[k] = \frac{ M }{ 2N } \left[ \frac{ Ue^{ i\beta_k }  - V }{ \cos( \alpha ) - \cos( \beta_k ) } \right] $$

Where:

$$ U = \cos( \alpha N + \phi ) - \cos(\phi) $$

$$ V = \cos( \alpha N -\alpha + \phi ) - \cos( -\alpha + \phi ) $$

Notice that there is a pluggable discontinuity at the bin value of a whole integer frequency.  When the frequency is very close to an integer value there is a more computationally robust form of the formula:

$$ Z[k] = \frac{ M }{ 2 }
\cdot \frac{ \sin \left( \frac{ \delta N }{2} \right) }{ N \sin \left( \frac{ \delta }{2} \right) }
\cdot \left[ \frac{ \sin \left( \frac{ \delta N }{2} + \phi \right) e^{ i\beta_k }  - \sin \left( \frac{ \delta N }{2} - \delta - \beta_k + \phi \right) }{ \sin \left( \frac{ \delta  }{2} + \beta_k \right) } \right] $$

Where:

$$ \delta = \alpha - \beta_k = ( f - k ) \cdot \frac{ 2\pi }{ N } $$

The derivations of these formulas can be found in my blog articles "DFT Bin Value Formulas for Pure Real Tones" and "An Alternative Form of the Pure Real Tone DFT Bin Value Formula".

I realize this is a reply to a very old post, but Community dredged it up and finding these formulas, and the subsequent inverses which give exact frequency formulas is what got me into DSP.

Ced


Decomposition of **3D** structuring elements for morphological operations I am struggling to implement a mathematical morphology toolset in an image processing package. I base my implementation on what I saw in MATLAB, and on several papers on the subject.

There seems to be abundant literature on morphological operation optimization through structuring element (strel) decomposition. For instance, one can get a tremendous speed bonus by using two orthogonal lines instead of a square as structuring element for dilation. Several papers give methods for optimization through strel decomposition:

Rolf Adams, "Radial Decomposition of Discs and Spheres," CVGIP:
         Graphical Models and Image Processing, vol. 55, no. 5, September
         1993, pp. 325-332.
Rein van den Boomgard and Richard van Balen, Methods for
       Fast Morphological Image Transforms Using Bitmapped Binary
       Images, CVGIP: Models and Image Processing, vol. 54, no. 3, May
       1992, pp. 252-254

etc...

However, all these publications are about 2D structuring elements. I could not find much on 3D decompositions.

Do you have any clues on how to decompose:

a 3D sphere. Not a ball that is used in 2D grayscale morphology, but an actual flat 3D sphere;
a 3D diamond?

I guess this depends on the digital distance transform that one is approximating on the 3d grid and there are various local connectivities possible. There is an implementation in ImageJ here.

It would also be good to verify if you are using a non-flat structuring element or a correct 3d structuring element. Read Matlab reference here. In the place of euclidean distance sqrt(x.^2 + y.^2 + z.^2) one could add the Manhattan distance using this. Quick 3d distance here.

For 3d structuring element decomposition one can see. But the 3d decompositions of a convex shapes like sphere into separable 2d lines is non-trivial, certain shapes are easier than others like the cube. One can refer here for efficient algorithms for spheres.


Modulation and Noise I am fairly new to studying communication systems.

When it comes to different digital modulation techniques, does the effect of an additive complex noise (a) introduced by the system, and (b) by the channel vary depending on the type of modulation used?

Moreover, is the noise introduced by the non-ideal system can be treated exactly the same as with the noise introduced by the channel (that is, the approach to work around and solve them are the same)?

This is up to how you model noise and channel.

In the classical, most-commonly-used models, the channel is actually noise-free and the noise is only added at the receiver.

If you have a different channel, that's kind of specific to how your special problem is being modeled, and we don't know that.

But since you say, you're new to communications systems:
The cases you're most likely to encounter are actually that the noise only happens at the receiver.

If that's not the case, look for the context in which you're working. For example, if you're in the context of diversity combining, and are considering a channel that has interferers, e.g. because  you're trying to implement optimum combining, then your overall noise is usually described by a covariance matrix – in which you don't really care about where the noise energy actually comes from (there's the assumption that receiver noise is uncorrelated and hence manifests as diagonal entries in the matrix).


Power spectrum of uniform white noise Given a white noise image $W_{i,j} \sim U[a,b]$ where each pixel is distributed uniformly in $[a,b]$, how would I go about computing its power spectral density? That is, I want to find $E[|\hat{W}_{i,j}|^2]$. Am I supposed to write it out through the autocorrelation function and then use the fact that pixels are independent of each other to split the expectation? Would this just yield a non-zero term at the DC? That is, the squared mean grey value $E[W_{i,j}] = \frac{a+b}{2}$.

The fourier transform of the two dimensional autocorrelation function should do it.


Wireless Channel Multiplicative Model Let's assume a SISO channel in the equivalent baseband, the standard channel model (discrete) is:
\begin{align}
y(k) = h(k)x(k) + z(k) \tag1
\end{align}
Where $y(k)$ is the detected signal, $h(k)$ denotes the channel, $x(k)$ is the sent signal and $z(k)$ is additive noise. If we consider reflections on buildings, cars, etc., why isn't there a convolutional model like in acoustics:
\begin{align}
y(k) = h(k)*x(k) + z(k)\text?\tag2
\end{align}

What's labeled $(1)$ in your question is a special case, the flat channel. It can be represented as a single coefficient.

In general, channels aren't flat, and we then need to apply $(2)$ instead. That's no different from acoustics.

So, your claim that $(1)$ generally applies to SISO channels is plain wrong.

However, when naming something "SISO", one usually means "in contrast to MIMO"; and in the context of MIMO, we basically need our channels to be flat, so that the MIMO channels can be represented as a $m\times n$ matrix of single coefficients between $m$ and $n$ transmitting/receiving antennas, respectively.¹

But, you say,

  if the individual channel is like $(2)$, how can we model it like $(1)$?

We can't. But we can apply multi-carrier techniques (mostly, OFDM), to first divide a wideband channel with a convolutive channel impulse response into many narrow channels, which individually are flat.

¹ In theory, one can do MIMO on multipath / frequency-selective channels instead of first dividing them in many flat subchannels, but as far as I know that's not done commercially. OFDM implementations are mature and cheap, and it's much easier to come up with working math that allows us to decompose the superpositioning MIMO channel when it's flat. In fact, there's no standard textbook I'm aware of that even does MIMO on frequency-selective channels.


CIELAB Color Clustering I am trying to understand if there is an inherent advantage to using CIELAB for color-based segmentation.

I segmented the cytoplasm from the image above, where the left image below is from RGB color clustering, and the right one is from CIELAB clustering, where a and b matrices were used.

You can see that there are very slight differences between the methods, where in this case, CIELAB clustering seems to work a bit better, but there are other examples where it doesn't.

So my question is, can you give me the reasons for using CIELAB?

If it's too much of a bother to explain, a good article reference would be nice, it would be better if it includes some mathematical explanation as well.

Generally speaking, a change of color space is essentially a nonlinear transform that stretches and distorts the coordinates. The transform is also continuous, except maybe on a singularity line or half-plane. It does not preserve the distances, but preserves the neighborhoods.

In the case of RGB > Lab, you can see it as a linear transform (change of basis, from RGB to XYZ) followed by a nonlinear transform of every coordinate independently (rescaling of XYZ to Lab). The topology is essentially preserved.

On the other hand, image segmentation algorithms do their best to delimit clusters of similar colors and separate the clusters from other clusters. So if you just deform the clusters, a good segmentation algorithm will adjust the boundaries and somehow undo the transform. If you repeat the experiment with XYZ or YUV coordinates, you should see no big difference either.

A stronger qualitative difference can be observed if the topology of the space is modified, for example by moving to the HSV space (Cartesian to cylindrical coordinates), of by performing some projection to lower the dimensionality such as XYZ > xyz. But beware of the physical interpretation of what you are doing.


How do I convert a real baseband signal to a complex baseband signal? I have radio telescope observations that have resulted in two real-valued signals (corresponding to the right- and left-handed circular polarizations).

The signals are sampled at rate $2B$, and provide a bandwidth of $B$.

I wish to get the full Stokes parameters from these signals, which requires me to convert these real signals into complex signals sampled at $B$ (half that of the real signals).

However, I don't know the procedure to do this. I just know that it is possible.

In general, how do I take a real-valued baseband signal and convert it to a complex-valued signal?

To convert a real signal sampled at rate $2B$ to its complex baseband representation (sampled at rate $B$), you want to map the frequency content in the range $[0, B)$ in the real signal to the range $[-\frac{B}{2}, \frac{B}{2})$ in the resulting complex signal. This can be done in a couple different ways:

Design a linear filter to approximate a Hilbert transform. Run your real signal $r[n]$ through the filter to yield the transformed signal $\tilde{r}[n]$. Use this to form the analytic signal:
$$
r_a[n] = r[n] +j\tilde{r}[n]
$$
$r_a[n]$ will contain only the positive frequency components of the original real signal $r[n]$; all of the negative frequencies will be zero (assuming a perfect Hilbert transformer; in practice, the effect will not be perfect). Thus, you have isolated the desired frequency band $[0, B)$.

Multiply $r_a[n]$ by $e^{-j\frac{\pi}{2}n}$ to effect a frequency shift of $-\frac{B}{2}$, shifting the desired frequency content to the range $[-\frac{B}{2}, \frac{B}{2})$. Then, decimate the signal by 2 by discarding every other sample. The result is a complex baseband signal sampled at rate $B$.
My preferred approach is a more straightforward implementation of the shift that you're looking for:

Multiply $r[n]$ by $e^{-j\frac{\pi}{2}n}$ to effect a frequency shift of $-\frac{B}{2}$, shifting the desired frequency content to the range $[-\frac{B}{2}, \frac{B}{2})$. The result is a complex signal that is centered in the appropriate place (the center of the band is at zero frequency).
Apply a lowpass filter to pass only the content in the range $[-\frac{B}{2}, \frac{B}{2})$ (to the extent required while meeting your application's antialiasing requirements).
Decimate the signal by 2 by discarding every other sample. The result is a complex signal sampled at rate $B$.

In practice, I always use some variant of strategy #2. You can make it even more computationally efficient by implementing the decimation as part of the lowpass filtering process, for instance with a polyphase filter.