Confusion about ensembles and averages in autocorrelation matrices I just noticed that until now I often don't cared about the scaling of the autocorrelation Matrix in Matlab. I then checked with the book Statistical Digital Signal Processing from M. H. Hayes. There it is defined like this:

First of all, why does the author note it like $x(0)x^*(0)$? Written like this wouldm't this mean first element of the vector x times first element of the vector x? And the expected value is then what $\frac{x(0)x^*(0)}{2}$?

$r_x(\tau)$ can be calculated as $E[X(t)X(t+\tau)]$ This is more what I would expect. A pointwhise vector multiplication of a $\tau$ shifted vector.

Here comes my actual question. For e.g. $x = [1 ,2, 3]^T$, $r_x(0)$ should become $\frac{1^2+2^2+3^2}{3}=\frac{14}{3}$ in my oppinion. Matlab shows for xcorr(x)

3.0000    8.0000   14.0000    8.0000    3.0000

$r_x(0)$ should be the value in the middle. But why don't they take the expected value? On this page they suggest that the autocorrelation matrix can be calculated with

z=autocorr(x)

and

Rxx=toeplitz(z)

Instead of autocorr, xcorr(x) should work too stated other posts. As expected this results in

Rx =

    3.0000    8.0000   14.0000    8.0000    3.0000
    8.0000    3.0000    8.0000   14.0000    8.0000
   14.0000    8.0000    3.0000    8.0000   14.0000
    8.0000   14.0000    8.0000    3.0000    8.0000
    3.0000    8.0000   14.0000    8.0000    3.0000

Again, no expected values. And this solution is also wrong because it is shifted. $r_x(0)$ is not in the upper left corner. I could shift it by hand, but is that how autocorrelation matrixes are calculated in Matlab? I checked with some other sources in the internet. Often the algorithm to calculate an entry of an autocorrelation contains just the sum.

EDIT: sorry that I didn't fully get your question. Now I got your points of conufion and try to provide a shortest possible answer:

First, you should better use xcorr function when estimating ACS $r_x[k]$ from data sequence $x[n]$. The output of xcorr will return the lag-0 at the middle sample as you recognize. Look at the input argument scaleopt which provides normalization to the computed output. The default is none ! and causes a confusion for you. You should select biased so that equation 7 below and xcorr output will be identical. Like this:

rx = xcorr([1,2,3],'biased'); % identical to (1/3)*xcorr([1,2,3])

Second,  do not use the raw output of xcorr as the input to the toeplitz function. As the output of xcorr will include lags for rx[-2],rx[-1],rx[0],rx[1],rx[2]. But you only need rx[0],rx[1],rx[2]. So you better call it like:

L = 3;
x = [1,2,3];
rxL = xcorr( x, 'biased');
Rx = toeplitz( rxL(L:end) );  % now yo uare ok ;-)

Now it could be argued that the Matlab's default xcorr scale option none is misleading. And yes, I would say so. In its default form it computes the convolution and leaves it unscaled...

Finally, please distingusih a theoretical expectation from practical estimation. When algorithms work on data, they estimate and they don't compute expectations, which is a theoretical tool. Hence, all expectations are replaced with estimations when algorithms are run for the practical applications...
It seems (to me) that your confusion is about three things.

Theoretical definition of the autocorrelation sequence $r_{XX}[k]$ of a random process $X[n,s]$
Practical estimation of it from a sample function $x[n]$ via ergodicity.
Interpretation of a random vector (including samples of a random process)

Consider the interpretation of a random process $X[n,s]$ as an ordered collection (set) of indexed random variables $X_n(s)$ for each time index $n$. Each random variable $X_n$ may, in principle, has its own unique probability distribution/density function $F_{X_n}(x)$. Consider the enormous implication of this on the number and type of all possible joint density functions.

To make things extremely simplified, we make the assumptions that all those random variables have the same PDF $F_{X_n}(x)$ for all $n$, and independent of each other, also further that first and second moments are independent of partciular time instants, which yields an i.i.d. , WSS random process $X[n,s]$ at hand.

The theoretical auto-correlation sequence (ACS) of this WSS random process is defined as:

$$ r_{X_n X_m}(n,m) = E\{X_n X^*_m\} = E\{X[n]X^*[n+k] \} = r_{X X}[k]   \tag{1}$$

The computation of this theoretical ACS $r_x[k]$ requires you to evaluate the integrals with PDF functions. This is possible with questions on the paper, but with applications in working systems, where you won't solve (mathematical) questions but perform statistical calculations, you need to estimate it from data, based on ergodicty assumption, which states that theoretical ensemble averages can bu estimates from practical sample function time-averages. Then you have

$$ \hat{r}_X[k] = \lim_{N \to \infty} \frac{1}{2N+1} \sum_{n=-N}^{N} x[n]x^*[n+k] \tag{2}$$

Note that $x[n]$ is a particular sample function (a realization) of the hrandom process $X[n,s]$. Note also that the infinite limits are replaced by data observation boundaries in practice.

Now consider the random vector $X(n) = [X(n),X(n-1),..,X(n-N+1)]^T$ which includes $N$ consecutive random variables from the random process $X[n,s]$, beginning at index $n$.

From this vector build up the outer product $X \cdot X^H$ , which includes a matrix of pairwise multiplications of random variables $X_i X_j$ as in

$$
RM =
  \begin{bmatrix}
X_1 X^*_1 & X_1 X^*_2 & ... & X_1 X^*_N \\
X_2 X^*_1 & X_2 X^*_2 & ... & X_2 X^*_N \\
... & ... & ... & ... \\
X_N X^*_1 & X_N X^*_2 & ... & X_N X^*_N \\
  \end{bmatrix} \tag{3}
$$

Now taking the expectation of this $RM$ matrix yields the theoretical auto-correlation matrix of the WSS random process as:

$$ R_X = E\{ X \cdot X^H \} = E\{RM\} = E\{\begin{bmatrix}
X_1 X^*_1 & X_1 X^*_2 & ... & X_1 X^*_N \\
X_2 X^*_1 & X_2 X^*_2 & ... & X_2 X^*_N \\
... & ... & ... & ... \\
X_N X^*_1 & X_N X^*_2 & ... & X_N X^*_N \\
  \end{bmatrix} \} \tag{4}$$

and this results in:
$$ R_X = E\{ X \cdot X^H \} = \begin{bmatrix}
r_x[0] & r_x[-1] & ... & r_x[-N+1] \\
r_x[1] & r_x[0] & ... & r_x[-N+2] \\
... & ... & ... & ... \\
r_x[N-1] & r_x[N-2] & ... & r_x[0] \\
  \end{bmatrix} \tag{5}
 $$

for real data it becomes:
$$ R _x =
 \begin{bmatrix}
r_x[0] & r_x[1]^* & ... & r_x[N-1]^* \\
r_x[1] & r_x[0] & ... & r_x[N-2]^* \\
... & ... & ... & ... \\
r_x[N-1] & r_x[N-2] & ... & r_x[0] \\
  \end{bmatrix} \tag{6}
$$

Again this theoretical ACM $R_X$ cannot be computed but is estimated from sample function data...

Estimation of the auto-correlation matrix $R_x$ can be done in a number of ways (both related) but at the basic level you see that what you should actually estimate is the auto-correlations sequence $r_x[k]$ for the lags $k=0,1,...,N-1$. Then from the Toeplitz and Hermitian symmetric property of $R_x$, you can construct it.

One way to practically estimate ACS $r_x[k]$ from the data $x[n]$ of length $L$ is :

$$ \hat{r}[k] = \frac{1}{L} \left( x[k] \star x^*[-k] \right) \tag{7}$$

This yields an estimate of ACS for the lags $-L < k  < L$, (you only need $N$ of those for $k=0,1,...,N-1$)

Alternatively, ACM can be estimated as a whole similarly to ACS from:

$$ \hat{R}_x = \frac{1}{K} \sum_{k=1}^{K} X(k) \cdot X(k)^H \tag{8}$$  where data vector $X(k)$ includes $N$ consecutive samples from $x[n]$ (specially processed at the edges)


Comparing FFT for increasing time intervals in a signal Assume a dataset where the rows represent signals: each signal $s$ is sampled at a sampling rate $f_s$ and available as an array of common length $T$. For an increasing sequence of integer times $t_i < T$, I would like to find out the information gain reached by considering the rest time $T-t_i$. Intuitively, I would like to find out a $t_{\dagger}$ for which the amplitude of the different frequencies in the signal has already converged. The motivation behind this is that I would like to potentially shorten the measurement time in the future, assuming the data distribution remains the same.

My approach has been to analyze the convergence in the FFT spectra of the increasing time slices, with the FFT computed from time $t=0$ to the corresponding $t_i$.

Are there any more adequate ways of doing this?

I'm afraid the premise of your approach is faulty as can be shown with a single counterexample.  Suppose your signal is a simple real pure tone sinusoidal.  When the interval length reaches 3 cycles, bin index 3 (zero based) of the DFT will be non-zero and all the rest of the DFT bins will be zero.  (Considering only the lower half of the DFT values because a real valued signal as the upper half is a conjugate mirror of the lower half.)  Then as your interval increases the DFT starts to smear with leakage.  When you reach 3 and 1/2 cycles, the smearing will be the greatest and the magnitudes of bins 3 and 4 will be nearly equal and their phases nearly 180 degrees apart.  When you get up to 4 cycles in your sampling interval, once again all the bin values are zero except for bin 4.  That behavior will repeat for each whole number of cycles your interval reaches and therefore there will never be any kind of convergence happening.

It is not clear to me from your question, but are you assuming your signals are mixes of steady pure tones?


Method for modifying frequency-domain signal to produce predetermined time-domain signal We are performing a convolution of audio signal, with 511th order FIR
filters, using fast convolution methods: overlap-add or overlap-save.
Filters are being designed in frequency domain several thousand times
per second. Filters need to be zeros padded to be used in order to be
used in fast convolution algorithms, so to acquire appropriately
padded filters we do:

frequency domain filter design, filters of 511th order
IFFT into time domain
zeros padding to 2048 length (1536 zeros)
FFT back into frequency domain (for multiplication with spectral
representation of input audio signal)

Computational complexity of such operation occurred be to high for our
needs. We're trying to figure out if there's a method to avoid
IFFF->FFT steps and reduce the overall complexity from O(NlogN) to
O(N). Any ideas if it is possible, and, if so, how can it be achieved?

Zero-padding is one domain is equivalent to Sinc (or Dirichlet) convolution in the other, which is O(n^2), so that won't save you.  For less computation, you might want to try different filter design methods (compositing existing FIR filters, IIRs, etc.) or finding a way to recompute your filter less often.


How to derive the results that averaging $N$ signals yields a $\sqrt{N}$-fold increase in signal-to-noise ratio? If I have understood correctly, averaging $N$ noisy independent segments or signals increases the signal-to-noise $\sqrt{N}$-fold.

How does one derive this result?

I will show how to calculate the SNR for the case of $N=2$ measurements; it is easy to extend the result to general $N$. Assume a signal $s(t)$ has power $S$, and the noise $n(t)$ has variance $\sigma^2$ and zero mean. Then, the signal $s(t)+n(t)$ has SNR equal to $S/\sigma^2$.

Now assume you observe $s(t)$ twice, each time with different, uncorrelated noise, $n_1(t)$ and $n_2(t)$, each with zero mean and variance $\sigma^2$. You average the two observations to get $$\frac{2s(t)+n_1(t)+n_2(t)}2 = s(t)+\frac{n_1(t)}2+\frac{n_2(t)}2.$$

The variance of both $n_1(t)/2$ and $n_2(t)/2$ is $\sigma^2/4$, so the total noise variance is $\sigma^2/2$ and the SNR is $2S/\sigma^2$, for an improvement of 2.

For some reason I don't understand, the wikipedia page on signal averaging defines the SNR as $S/\sigma$. Under that definition, the improvement will indeed be $\sqrt{N}$.


Fourier Transform of infinite sum I'm looking at last year's exams, and I found an exercise I can't solve:
(Roughly translated)

Consider $x(t)=\sum\limits_{k=1}^{+\infty}\left(\frac{1}{2}\right)^k \cos(k2\pi t)$ the input to a LTI system, with the transfer function shown in the picture below

a) Determine (and sketch) the Fourier transform of x(t) and y(t).

b) Determine the Fourier coefficients of y(t). Justify.

Edit:

The only thing that comes to mind is using the Linearity of the Fourier Transform and calculate the Transform of $\left(\frac{1}{2}\right)^k \cos(k2\pi t)$, k dependant, and $X(jw)=\sum\limits_{k=1}^{+\infty}X_k(jw)$, so it would be something like
$$
X(jw)=\sum\limits_{k=1}^{+\infty}\left[2^{-k}\pi\left[\delta(\omega-2k\pi) + \delta(\omega+2k\pi)\right]\right]
$$

However, I have never seen a Fourier Transform that looked like this, so I'm not quite sure that's correct, or, if it is, if it can be used (efficiently) for calculations, in this form.

After I get $X(jw)$, I can get $Y(jw)=X(jw)H(jw)$, so right now I'm pretty much just stuck in calculating $X(jw)$. Any help is much appreciated, and feel free to ask for clarification, if something I wrote is not clear enough.

i think it's pretty obvious that only the component with angular frequency of $4 \pi $ gets through the LTI system to the output.  so the output should be $\frac{1}{4} \cos(4 \pi t) $.


find power of the signal I have this signal:
$$x(t)=\sum \mathrm{rect}\left(\frac {t-kT_0} {20}\right)+\mathrm{rect}\left(\frac {t-kT_0} {40}\right)+\mathrm{rect}\left(\frac {t-kT_0} {80}\right)$$
with $T_0=140$, pass it through the filter $H(f)={\mathrm{rect}}\left(\frac f F\right)$.
What is the value of $F$ such that power is $P=1$?

Fourier Trasform is: $\sum C_n\delta(f-nf_0)$
$C_n=\frac 1 {T_0}[sinc(20nf_0)+sinc(n40f_0)+sinc(n80f_0)]$

if call y(t) the signal output, I have:
$Y(f)=X(f)H(f)$.

the power spectral density is:  $P(f)= \sum (|C_n|^2\delta(f-nf_0))|H(f)|^2$
 with  $P_x(f)=\sum (|C_n|^2\delta(f-nf_0) $
don't know find F

if don't error for n=0 $C_n=1$
 so for $0<F<\frac 1 {T_0}$ the power should be 1. but I'don't sure

Since this is clearly homework, and the limits on the $\sum$ are crucial to figuring out whether there's an answer:

Hint: Parseval's Theorem states that power in time domain is equivalent to power in frequency domain. What is the power in $x(t)$? Is $x(t)$ periodic? If so, what's the period?

As usual, it often helps to make a drawing (that's, by the way, the minimum I'd have expected from your question!) to get a feeling for the signal you're dealing with.


Upsampling/Interpolation and Downsampling/Decimation I am having trouble figuring out exactly what is happening the process of downsampling/upsampling and the output in the example given below:

So the 100Hz sine wave is sampled at 2000Hz resulting in a fft output mirror about 1000Hz, the nyquist frequency.
It is then downsampled, effectively sampling this signal at 500Hz (factor of 4)/ or is it more a appropriate to say the signal now has characteristics like the 100Hz wave was sampled at 500Hz? This is achieved by discarding samples.
This wave is then upsampled by four, however I do not see what this does to benefit the signal? I understand afer this it is low-pass filtered to remove the new alias images due to downsampling, but surely you you could low-pass filter after the downsampling and still end up with the same frequency output, i.e. 100Hz? Or does it allow a 100Hz output with increased data points, which is better?

Thanks for taking the time to read my question and any help would be greatly appreciated!

So the 100Hz sine wave is sampled at 2000Hz resulting in a fft output mirror about 1000Hz, the nyquist frequency.

no. The FFT doesn't give you that. It's the spectrum as visualized here. It's an interpretation of the FFT (which is just an implementation of the DFT). For the DFT, it's just a single periodic function with a period of 5. That simple. No repetitions.

The repetitions in the visualization stem from the fact that mathematically, discrete spectra are periodic, and that the spectrum graphs were chosen to have frequency axis longer than a single Nyquist bandwidth. Again, the DFT doesn't care about what it means – it will just "show" that there is a periodicy of frequency $\frac 15$.

  It is then downsampled, effectively sampling this signal at 500Hz (factor of 4)/ or is it more a appropriate to say the signal now has characteristics like the 100Hz wave was sampled at 500Hz? This is achieved by discarding samples.

You're starting to get into the right mindset. Yes, this is like a 100 Hz signal sampled at 500 Hz – count the samples! There's a repetition every five samples!

Generally, I recommend just "dropping" knowledge of the original frequency at the ADC. Really, it makes things easier. You're sampling a 100 Hz tone with 2000 Hz, so your period is 20 samples, and the relative frequency is $\frac 1{20}$. Then you decimate by 4 — leading to a relative frequency of $4\frac1{20}=\frac15$.

The DSP algorithms don't have any of that "original frequency" knowledge; as said, the DFT really doesn't care whether it's a 100 Hz signal sampled at 2 kHz, or it's a 4 MHz signal sampled at 80 MHz – all that counts is how long a period is, in samples, because all the DSP only sees your signal as sequence of numbers. One after the other – no matter how fast the orignal samples came in.

  This wave is then upsampled by four, however I do not see what this does to benefit the signal?

Wrong question to ask. It's done to illustrative interpolation.

  I understand afer this it is low-pass filtered to remove the new alias images due to downsampling, but surely you you could low-pass filter after the downsampling and still end up with the same frequency output, i.e. 100Hz?

Nononononono!

You didn't learn your lesson; yes, your signal after decimation could be filtered, but that would have no effect, because it still fulfills the bandwidth criteria set by Nyquist. That's why there's no anti-aliasing filter here – the signal is already sufficiently band-limited as is.

However, after interpolation, you will need an anti-imaging filter.

Because the inserting of zeros did exactly that, adding images where there should be none.

  Or does it allow a 100Hz output with increased data points, which is better?

Again, not the question. This is an illustration, and the fact that you wanted to omit the anti-image filter shows that it's a good illustration: compare $|X_i|$ and $|X_u|$. You only want to see what's left in $X_i$, not all the spectrum that is in $X_u$!


Golay complementary sequences in 802.11ad The Golay complementary sequences consist of two sequences of the same length $N$ such that there acyclic autocorrelation sequences have sidelobes equal in magnitude but opposite in sign. So, when the two autocorrelations are summed up, the result has a peak of $2N$ and zero sidelobe level. As an example, $Ga_{128}$ and $Gb_{128}$ are two complementary Golay sequences of length 128. If I sum up the autocorrelations of $Ga_{128}$ and $Gb_{128}$, the sidelobes are exactly zero.

The IEEE 802.11ad protocol uses Golay sequences for channel estimation. However, it doesn't simply use the basic $Ga_{128}$ and $Gb_{128}$ sequences. Instead it concatenates these sequences to generate following longer sequences:

$Gu_{512} = [-Gb_{128}\phantom{1} -Ga_{128}\phantom{1}  Gb_{128}\phantom{1} -Ga_{128}]$
$Gv_{512} = [-Gb_{128}\phantom{1}  Ga_{128}\phantom{1} -Gb_{128}\phantom{1} -Ga_{128}]$

(a) However, when I add the autocorrelations of $Gu_{512}$ and $Gv_{512}$, I don't get zero sidelobes. So, what is the value of such longer sequences?

(b) Regardless of the sidelobe performance, I would assume that only $Gu_{512}$ and $Gv_{512}$ sequences should be transmitted. However, the IEEE 802.11ad single-carrier (SC) modulated signal transmits a third sequence consisting of $-Gb_{128}$. Any explanation why this different length signal should be transmitted for channel estimation apart from $Gu_{512}$ and $Gv_{512}$?

Any explanation would be greatly appreciated.
~ryan

For part (a), I discovered that the IEEE 802.11 protocol doesn't recommend adding auto-correlations of $Gu_{512}$ and $Gv_{512}$. Instead the correct method is to autocorrelate the first 256 bits of each of these sequences with the next 256 bits. So, $Gu_{512}$ and $Gv_{512}$ are not complementary Golay sequences themselves. But, essentially, both $Gu_{512}$ and $Gv_{512}$ consist of two 256-bit Golay complementary pairs.

~ryan


what's the baseband FFT size of LTE? In standard of 3GPP LTE physical layer(3GPP TS 36.211), the FFT size for baseband signal generation for both Dl and UL is 2048 and CP size is also determined. However there is a table in several resources said we have different FFT and CP size. E.g.

http://lteuniversity.com/ask_the_expert/f/59/t/405.aspx

How is it possible?

Read the book 4G Lte Advanced by Stefan Parkvall and Eric Dahlman. According to the book,
FFT size of 2048, with a corresponding sampling rate of 30.72 MHz, is suitable for the wider LTE carrier bandwidths, such as bandwidths of the order of 15 MHz and above. However, for smaller carrier bandwidths, a smaller FFT size and a correspondingly lower sampling rate can very well be used.(Lowering the sample rate simplifies the implementation in dual-mode HSPA/LTE terminals)

Also,nothing forbids the implementation of an OFDM modulator as a set of parallel modulators,
Also, nothing prevents the use of a larger IFFT size, for example a size-2048 IFFT size, even in the case of a smaller number of OFDM subcarriers.

Tcp= 160·Ts and Tcpe=512.Ts. Change in CP induces change in the total number of REs(6X12, instead of 7X12). Note that #sub carriers remain 12 in both cases and FFT size is dependent on BW allocated which is expressed in #sub carriers.


Derivation of downsampling in the frequency domain I'm having some trouble with this derivation

$$s_d = \begin{cases}1&\text{for}\quad m\quad\text{multiples of}\quad D\\ 0 &\text{otherwise}
\end{cases}$$

then $s_d$ is somehow rewritten to:

$$\frac{1}{D}\sum_{k=0}^{D-1}e^{jk2\pi \frac{m}{D}} $$

Can someone explain this?

EDIT:

Basically the same derivation as done here
Frequency Representation of Downsampled Signal

The expression

$$s_d[m] = \frac{1}{D}\sum_{k=0}^{D-1}\exp(j2\pi \frac{km}{D})$$

is the Fourier series expansion of the periodic function $s_d[m]$ in the discrete case (i.e. only up to a certain amount of frequencies).

$s_d[m]$ is a also called comb function, and I have written more about it (including the relation you ask for) in one of my articles.

Furthermore, here's some example Python code that verifies the correctness:

m = np.arange(64)
D = 16
sdm = 1/D*sum(np.exp(2j*np.pi*k*m/D) for k in range(D))
plt.plot(m, sdm.real)

The picture shows the Comb function, where every 16 samples the values becomes one, else zero.


Demodulating a LoRa data symbol LoRa integrate data into symbols by selecting the starting frequency of an up chirp. The resulting up chirp traverses across BW and arrives at the same place it started. At the demodulator, this chirp is multiplied by a down chirp followed by the fft to extract the data bits.
I tried to do this myself on paper but I can't help but find two different fft bins. I understand if both -ve & +ve bins are same, the demodulation is successful however, they don't aren't. I illustrate this at my best below.
As per DSP principles, the multiplication of received and locally generated base band chirps is the addition of the two instantaneous frequencies of the chirps in red and purple. Then, I have drawn a line which represents their summation in yellow. As we can see, the multiplication produces two products (in yellow), one positive and one negative frequency component that are not equal.

But according to literature , multiplication with a down chirp should produce only a single frequency bin, meaning, the positive and negative frequency products resulting from the multiplication must be of same magnitude. But, I see two different products here. How is that possible?
Would be so much helpful to me if someone can point out where I am wrong.

The resulting up chirp traverses across BW and arrives at the same place it started.

not quite exactly right, but very close: it ends up in the DFT bin before the bin it started from. This makes the chirps (frequency domain) cyclically shifted versions of the straight-up "prototype chirp". (my name invention)
Maybe this way of looking at it helps intuition:
Let's assume the system is synchronized in time¹ and let's look what happens when you multiply the unshifted downchirp with the upchirp.
You're multiplying a complex sinusoid with a frequency slope that is the exact inverse of the slope of the other sinusoid. In other words, you're multiplying two sinusoids which only differ in the sign of their exponents. Now, multiplying two exponential functions with the same base ($e^\cdot$) will lead to a new exponential function with the sum of the exponents of the factors as exponent. Yay! That means up- and down-chirp cancel out, you get a constant $e^{j0}=1$.
Now, if the chirp is shifted, the sum of the exponents doesn't cancel - but it becomes a linear function – as argument of the complex exponential: a tone!
A tone that on top of that falls in the raster of the DFT. So, a single peak after FFT.

As per DSP principles, the multiplication of received and locally generated base band chirps is the addition of the two instantaneous frequencies of the chirps in red and purple

yep.

Then, I have drawn a line which represents their summation in yellow.

And here you simply forgot to realize that your positive frequency is outside Nyquist bandwidth. You need to subtract the sampling rate, and you'll see that it aliases to the same frequency as the rest of the yellow line.
Don't forget that your in discrete time and frequency domain: the spectrum must be imagined to repeat infinitely every multiple of the sampling rate.

¹ Unsolicited words of advice: if you find literature that proposes cool signalling schemes and uses the sentence Assume the system is synchronized, don't invest too much time into that scheme until you find information on how that can be synchronized... synchronization is one of the hard parts, and it has broken many a system's promises.


Is there a way to protect your signal? I'm looking for a way to shield or protect produced signal or wave. In mechanics suppose you run the wave in tunnel. Is there a way to do it in transmitting the signals or waves created by mobile phone or wifi or in the internet?

Signal integrity is a full area of electrical engineering and signal engineering.
A few examples :

To protect a propagating signal in free-space you can emit it louder (=more transmitting power, hence increasing signal to noise ratio (SNR)).
Numerical transmissions often include error correcting codes.
To protect propagating signals in cables, most of them are electromagnetically shielded.

I may detail my answer according to what you want to know exactly.


How can I generate audio wave files of harmonic series/partials for a specific frequency in Matlab this is the only code I know, but it is for one specific frequency and not for its harmonic series, I am trying to get all the partials as individual files, eg. If I want to generate 40hz sine wavwrite file and its harmonic series which would be 40-80-120-140hz etc, about 551 individual audio sine wav files all at harmonic series frequency values in one go generated with an Matlab Code...

fs = 44100; % sampling frequency
nbits= 16;
file = 'output.wav'
d = 40; % duration
f = 60; % wave frequency
phi = pi/2; % phase shift - will produce cosinusoid for pi/2
t = linspace(1/fs, d, d*fs); % time vector
y = sin(2*pi*f*t + phi);
plot(t,y);
wavwrite( y, fs, nbits, file)

use a for loop to save each frequency into a different named wav file, use string manipulation command num2str() of Matlab. The following will help:

clc; clear all; close all;

Fs = 44100;         % Sampling frequency
F0 = 200;            % base frequency
t = [0:1/Fs:0.1];   % time index
k = 0;              % harmonic index

for fsine=F0:F0:Fs/2
    x = sin(2*pi*fsine*t);
    filename = ['wav' num2str(k)];
    wavwrite(x, Fs,16, filename);
    k = k+1;
 end


Reconstruction of audio signal from Spectrogram I have a set of songs for which I extracted the magnitude spectrogram using a Hamming Window with 50% overlap. After extracting the spectrogram, I did some dimensionality reduction using Principal Components Analysis (PCA). After reducing it to lower dimensionality, I reconstructed the spectrograms from lower dimensions. So now, there would be some error between the original spectrogram and the reconstructed spectrogram. I would like to convert this spectrogram back to the audio signal and play it, so that I would be able to know when reconstructed from lower dimensions, how does the audio sound.

Is there any function available in say Matlab. to convert a magnitude Spectrogram to an audio signal ??

If the spectrogram was computed as the magnitude of short time fourrier transforms from overlapping windows, then the spectrogram contains implicitly some phase information.

The following iterations do the job :

$$x_{n+1} = \text{istft}(S\cdot\exp(i\cdot\text{angle}(\text{stft}(x_n))))$$

$S$ is the spectrogram, $\text{stft}$ is the forward-short time Fourier transform, $\text{isft}$ is the inverse-short time Fourier transform.


Isolating frequency-specific oscillation from transient I am performing time-frequency analysis on electrophysiological data with complex morlet wavelets. Each decomposition in the attached figure shows total power following the onset of a stimulus (top row) or preceding an eye movement (bottom row). There is a broadband, transient increase in power following stimulus onset.

I am interested in determining whether there are frequency-specific changes in power (i.e., 'oscillations') during the 0-200 ms epoch as well; therefore, I have parsed the total power into phase-locked and non-phase locked components. Because variability in response latencies may bleed into the estimate of non-phase locked power (David et al. 2006),

I want to compute adjusted power as recommended by these authors (pg 1589). I am unsure of how to implement this orthogonalization in Matlab-- it is slightly beyond my mathematical savvy.

I attached a snapshot of the related paragraph in the paper as well.

Is anybody familiar with this computation?

The notation and terminology in this and the closely related paper could have been more simplified I think.

To understand what's going on there, we need to start from page 1583. Specifically:

$y(t)$ is basically your measured signal (the output from each neuronal population). A scaled version of that is what you would sense via some modality (in an ideal world).
$s(\omega, t)$ is the frequency-time representation of a $y(t)$. This is complex.
Given that this is about Event Related Responses (ERR):

$g(\omega, t)_T$ is the average absolute value of $s(\omega,t)$ after a number of trials. This is real valued.
$g(\omega, t)_e$ is the average of $s$ times the average of the conjugate of $s$ (both of which will be complex, the result of which however will be real).
These look similar but their magnitude is not.

Then the explanation of the breakdown of that "power" follows. What this means is that, based on the assumptions in the paper, we think that the neuronal activity is due to three things: a) Background activity, b) Activity because of connections between neuronal populations and c) Induced activity by the stimulus. Furthermore, these three components are assumed orthogonal. This means that if one "happens", the rest do not. This is important for later on.
WHERE/WHEN this holds then we could estimate the background activity by looking at the spectrum BEFORE the stimulus onset, we could estimate the activity due to the stimulus (the part that is phase locked to the stimulus), because we know the "shape" of the stimulus and therefore what is left should be due to structural connections.

These are captured in $g(\omega, t))_i$ and $g(\omega,t)_T$.
Notice here that the baseline is only $g(\omega)$

So, when you arrive at equation 23, what you see is the explanation of that "adjustment" where:

$g(\omega,t)_T$ is the frequency-time representation of the averaged stimulus response after $T$ trials as before.
$g(\omega,t)_e$ is the frequency-time representation "power of average" as before
$\hat{\eta}$ now denotes the proportion of the frequency-time representation of $g(\omega,t)$ that is due to $g(\omega,t)_e$. The "trick" here is this, you have already assumed that the three components are orthogonal. WHERE/WHEN this holds, the "pseudoinverse" ($g(\omega, t)^+$) of a given $g(\omega, t)$ multiplied by the "shape" of the stimulus response $g(\omega,t)$ over the trials would be zero and you would be left with the "shape" of the stimulus response. So, that $\hat{\eta}$ business there (and the subtraction) is what allows you to do that "separation".

Perhaps it is now clear how to produce that $g(\omega,t)^+$. Provided that you have a given $g(\omega,t)_e$, invert it with pinv (or anything else equivalent to a pinv in another platform).

But, I wonder how effective would that be for what you are trying to do. Irrespectively of this specific technique, a "pulse" or a step or any abrupt change of a signal in the time domain usually has a broad representation in the frequency domain (whether that is Fourier or Wavelet, the point is that you need more components to synthesize the step). You are looking for "oscillations" in a 0-200ms interval. What is the frequency of the oscillations you are looking for?

Hope this helps


Two-dimensional wavelet analysis Could you tell me why wavelets without scaling functions can't be used for two-dimensional analysis? What is the role of scaling function?

,,A wavelet'' actually is two functions that are orthogonal to each other. That way, a signal can be split into two parts. One part is the wavelet proper, while the other is usually called ,,scaling function'' (sometimes ,,mother wavelet'' as well). When you have just one of both, you can not perform a meaningful wavelet transform.


Dynamic convolution vs Volterra series I'm trying to understand how a dynamic convolution model relates to something like a Volterra series, and what kinds of effects the latter can capture that the former can't (and vice versa). By dynamic convolution, I mean a system where the impulse response varies with the amplitude of an impulse.

I care only about dynamic convolutions that are well-behaved and "smooth" in some sense, i.e. where the impulse response varies continuously or smoothly as a function of the amplitude.

I have this intuition that a Volterra series with nonzero coefficients only on the diagonals are somehow related to dynamic convolutions. In other words, a series that reduces to the following:

$y = h_0 + h_1 \ast x + h_2 \ast x^2 + h_3 \ast x^3 + ...$

where there are never any terms like $x[t] \cdot x[t-1]$ and so on.

I'm curious if I'm on the right track with that, and in general, whether any dynamic convolution can be transformed into a Volterra series in that form or vice versa.

A flavor of dynamic convolution (is that a trademark by the way?) has a different impulse response $g_i$ associated with each range of instantaneous input. A number of ranges can be defined by fuzzy membership functions $f_i(x)$ (Fig. 1).

Figure 1. Amplitude ranges that each use a different convolution kernel.

Omitting time indices, the input $x$ and output $y$ are related by:

$$y = \sum_{i=1}^Ng_i\ast \left(f_i\small(x\small)\times x\right).$$

In the notation used, convolution has higher precedence than summation. Reading the signal flow from the above, $f_i(x)\times x$ is a memoryless waveshaper. It is followed by regular convolution with $g_i$, and finally, the outputs of the convolutions are summed. For a practical, limited absolute input amplitude, well-behaved waveshapers, such as the ones needed here, can be approximated to arbitrary precision by polynomials. With the appropriate polynomial series with coefficients $a_{i,j}$:

$$ = \sum_{i=1}^N g_i\ast\left(\sum_{j=0}^\infty a_{i,j}\times x^j\right),$$

and reorganizing:

$$= \sum_{j=0}^\infty\left(\sum_{i=1}^N a_{i,j}\times g_i\right) \ast x^j.$$

Recognizing $h_j = \sum_{i=1}^N a_{i,j}\times g_i$ we have your diagonal Volterra series:

$$= \sum_{j=0}^\infty h_j \ast x^j.$$

Unfortunately, polynomials are not very good at approximating functions of compact support, like those needed here (at least before combining the ranges), so for a reasonable level of error a diagonal Volterra series approximation of dynamic convolution may require taking quite large powers of the input signal, which is computationally expensive and can lead to numerical problems.


concurrent multiple sound sources seperation for TDoA estimation consider a room with some sound resources and several mics, these sound resources almost concurrently emit sound signals. As for one mic, it may pick up sound signals from different sound sources which are overlapped.

I am just asking for a brief idea of the separation of these overlapped sound signals for the further process such as TDoA estimation.
I have no idea which is the current mainstream technology to deal with this kind of problem?
I was thinking about the acoustic signal's encoding or modulation so mics can identify which signal from which sound source.

Thank you in advance

Following your idea, I tried 2 chirp signals with frequency division in Matlab, then use band filter+cross correlation algorithm, the result is quiet promising.

So instead of the chirp signal, did you just use sine wave signal(tone)?

Shuo,

Here is a couple of insights for your system:

1.- The different sound sources must play different signals (beacons), that may be resolved at the receiver device. A simple example would be to transmit a set of tones (one by speaker) with different frequencies. A slightly more complex solution would be to transmit other kind of signals with low cross-correlation (e.g. Gold, Kasami or Zadoff-Chu sequences).

A different approach at all would be to transmit actual modulated messages that somehow don't interfere each other at the receiver. This is usually even more complex.

2.- Regarding the method for distinguishing the different signals at the receiver, the typical solution is to use a set of matched filters (one filter for each beacon signal). Then, determine the TDoA by measuring the differences between the time indexes in which you find peaks at the matched filters' outputs.

On the other hand, if you really need to use modulation/encoding, a simple set of matched filter won't be enough, you will probably need to implement actual demodulation schemes, which may include synchronization algorithms and other practical restrictions.

According to my personal experience, and considering your introduction lacking much detail, you may successfully use a set of (frequency separated) tones at the transmitter, and a set of matched filters at the receiver. At least as a first approach, and if the performance is not good enough for you, you may use more advanced techniques. I honestly believe that in most cases no further cumbersome techniques are required.

I hope this helps!


How to Estimate Carrier Frequency & Amplitude for Several Overlapped Template Signal As is shown in the picture,
a template signal is delayed and overlapped in time domain,
the shape of  template signal is already known
The positions of each template are already known.

Each template has its own carrier frequency, and amplitude.
Or the carrier has a Gaussian distribution $\mathcal{N}$($\mu$,$\sigma^2$).

How to estimate them？ Please provide some recommendations, thank you.

If the auto correlation of the signal is sharp enough, you can do Matched Filter and search for local extreme points.

Yet it seems the figure you'e displaying is in the frequency domain.
But we can treat the frequency domain as it was time domain and work with the same "Trick".

Namely searching for the frequency in the frequency is like searching for the delay in the time domain.
Since the signal is known, its multiplication by a carrier is known whihc will be used as the Matched Filter.
Applying it will yield 3 peaks each with its own "Delay" which corresponds to its frequency.
Going back to that frequency bin (Or calculating the DTFT at that point) will yield the amplitude.


DPSK vs PSK Error Probability BACKGROUND

The bit error probability for BPSK under AWGN is easily derived from tail probabilities of Gaussian distributions and results in

$$P_e = Q\biggr(\sqrt{\frac{2E_b}{N_o}}\biggl)$$

The equivalent bit error probability for DBPSK is given as follows but much more complicated to derive:

$$P_e = \frac{1}{2}e^{-E_b/N_o}$$

A complete derivation for the DBPSK case is here:

http://staff.ustc.edu.cn/~jingxi/Lecture%209_10.pdf

With the same formula and plotted in comparison to BPSK on Wikipedia (https://en.wikipedia.org/wiki/Phase-shift_keying#/media/File:DPSK_BER_curves.svg):

MY QUESTION

I incorrectly thought I could simplify this derivation by extending the simpler BPSK $P_e$ result through understanding what occurs when you multiply two signals with independent noise  (Matt L has provided that here: SNR After Multiplying Two Noisy Signals), since such a product results when performing non-coherent demodulation for DBPSK.

I show this in the block diagram below:

This is the non-coherent structure for DBPSK demodulation. The transmitter is also differentially encoded to minimize error propagation (so that errors always occur in pairs rather than propagate until the next transition).

Here we can see that given an input DBPSK signal with $SNR = SNR_1$, the signal after being delayed one bit period $T$ will also have $SNR = SNR_1$, but the noise component will be independent (assuming AWGN, the noise is one symbol period is independent of the noise in the next symbol period). With reference to Matt L's result linked above, the predicted SNR at the output of the multiplier would be:

$$SNR_2 = \frac{SNR_1 SNR_1}{SNR_1+SNR_1+1}$$

For real signals, the frequency at the output of the multiplier is the sum and the difference of the input frequencies, so in this case the difference is baseband signal of interest while the sum is the double of the carrier that we filter out with the low pass filter (LPF). The power of both the signal and noise components would be effected the same way in this process, so the SNR at the output of the LPF would still be $SNR_2$.

Note for SNR>>1, $SNR_2$ approaches $SNR_1/2$ or 3 dB worst.

Given this, combined with the double error property that a single bit error always results in 2 errors assuming we use differential encoding in the transmitter- I can convince myself that the predicted bit error rate for DBPSK would be as follows (reducing the SNR by 2 when SNR >> 1 and doubling the resulting P_e) but from the detailed derivation this is clearly incorrect. I understand the detailed derivation- my question is not with that but what is the flaw with this alternate approach?

$$P_e = 2Q\biggr(\sqrt{\frac{x}{2x+1}}\biggl)$$

where $x = \frac{2E_b}{N_o}$

It is interesting to note that for higher order M-PSK, this 3 dB result does match (note the difference between QPSK and DQPSK in the plot above). Perhaps this is a clue that real versus complex is a factor?

The flaw in the problem posed starts right from the block diagram  that you have posted: it is not the case in DPSK receiver that the RF signal is multiplied with a delayed version of itself and then passed through a low-pass filter followed by slicing etc.

Instead, the RF signal is "demodulated" through a variant of a standard coherent QPSK receiver whose local oscillator is assumed to be synchronized in frequency to the incoming carrier signal but not necessarily synchronized in phase. Thus, the incoming $$r(t) = g(t)\cos(2\pi f_c t+\theta) + n(t), 0 \leq t < T,$$ where $g(t)$ is the baseband pulse, $T$ the bit duration, $\theta$ the unknown RF phase, and $n(t)$ the AWGN (or band limited AWGN if you like to think that the incoming signal has been bandpass-filtered, e.g. as in a tuned amplifier in an IF stage of a superheterodyne system), is multiplied by $2\cos(2\pi f_c t)$ in the I branch and separately by $-2\sin(2\pi f_c t)$ in the Q branch. Since
\begin{align}
2\cos A \cos B &= \cos(A+B) + \cos(A-B)\\
2\sin A \cos B &= \sin(A+B) + \sin(A-B)
\end{align}
we get that the signal parts of the two mixer outputs are
\begin{align}
2\cos(2\pi f_c t) g(t)\cos(2\pi f_c t+\theta) &= g(t)\cos(4\pi f_c t+\theta) + g(t)\cos(-\theta)\\
-2\sin(2\pi f_c t) g(t) \cos(2\pi f_c t+\theta)&= -g(t)\sin(4\pi f_c t+\theta) - g(t)\sin(-\theta).
\end{align}
The double-frequency terms in the mixer outputs can be eliminated explicitly by low-pass filtering the results before doing anything else, or they get eliminated implicitly when we do matched filtering on the two mixer outputs. Now, matched filtering is typically done by via correlation
(multiply by $g(t)$ and integrate from $0$ to $T$) and thus we effectively get a complex-valued decision variable $Xe^{j\theta} + \mathcal CN(0,\sigma^2)$ at the end of the integration period. We can't make a decision just on this decision variable because we don't know the value of $\theta$; depending on what $\theta$ is, the decision might be wholly bass ackwards from what it should be! But what we can do is save the decision variable and use it $T$ seconds later to answer the question

  Did the incoming signal change phase from the unknown $\theta$ during $[0,T)$ to (also unknown) $\theta+\pi$ during $[T,2T)$, or did the phase stay the same unknown $\theta$ during $[T,2T)$?

The final expression for the error probability of DBPSK is quite simple:
$$P_{\text{e, DBPSK}} = \frac 12 e^{-E_b/N_0}$$
but, as you say, is harder to derive than the $Q
\left(\sqrt{\frac{2E_b}{N_0}}\right)$ error probability for coherent BPSK.

Turning to the question of why the trick of using the coherent BPSK error formula as a guide and simply replacing the SNR for BPSK with the SNR for DBPSK in the error probability formula doesn't give the right answer for the DBPSK error probability, the issue is that the two systems make decisions differently and there is no obvious reason why the error probabilities of the two systems should be given by the same formula in terms of SNR. Note, for example, that $P_e = Q(\sqrt{2E_b/N_0})$ for coherent BPSK but only $Q(\sqrt{E_b/N_0})$ for coherent FSK.

Finally, note that for $x > 0$,
$Q(x)  < \frac 12 \exp\left(-\frac{x^2}{2}\right)$
and so, setting $x=\sqrt{2E_b/N_0}$, we get that
$$P_{\text{e, BPSK}} = Q(\sqrt{2E_b/N_0}) < P_{\text{e, DBPSK}} = \frac 12 e^{-E_b/N_0}$$


Significance of an impulse in the frequency domain I know that $X(f)$ gives the amplitude associated with the frequency component $f$ of a signal $x(t)$.

Now, a sinusoidal signal in time $x(t) = A \cos (2 \pi f_0 t)$, has a Fourier transform $X(f) = \frac{A}{2}[\delta(f-f_0) + \delta(f+f_0) ]$.

My question is that the Dirac Delta funcion tends to $\infty$ at $0$. Then, multiplying it by $A/2$ should also result $\infty$. If this was the case, then what is the significance of $X(f)$ for $x(t) = A \cos(2 \pi f_0 t)$.

More specifically what is the significance of the impulses (having infinite amplitudes) in the frequency domain?

Thank you. :)

The Dirac delta is not strictly a function but a distribution. The Dirac delta is such that $\delta(x)=0 \ \forall x\neq0$ and it has to meet the following restriction:

$$\int_{-\infty}^{\infty} \delta(x) \ \mathrm{d}x = 1$$

This means that the unit impulse must integrate $1$ over all the real numbers.

Let's define this other "function" $$\tilde{\delta}(x)=2\delta(x)$$

As we can easily see, this new "function" must integrate $2$:

$$\int_{-\infty}^{\infty} \tilde{\delta}(x) \ \mathrm{d}x = 2\int_{-\infty}^{\infty} \delta(x) \ \mathrm{d}x = 2$$

And this is what we get when multiplying an impulse by a constat number. The area under the "function" changes and, due to the Dirac delta being non-zero only at the origin, this means that the impulse must change its height (if we can call it so) in order to change the area it integrates.

In the frequency domain, these constants act as weights that determine something like the relevance of a frequency in a given time-domain signal. If we have a small (i.e. it integrates a small area) impulse in a determined frequency, then there is a pure sinusoid of that frequency present in the signal whose amplitude is rather small. On the other hand, a large impulse would correspond to a pure sinusoid with a large amplitude.


Distortion in sound after multiplying frequency spectrum by constant I make a simple sound equalizer that operates in frequency domain and lets user to adjust frequencies in sound by using 4 sliders. The first one responsible for 0 - 5kHz, the fourth one for 15-20kHz.

Steps are as follows:

I read wav file and store it in float array
I perform complex fft on that array (separately for left and right channel)
I multiply real and imaginary parts of bins representing 0-5kHz frequencies (both positive and negative) by 1.1 3.981 to increase these low frequencies by 10% 12dB in the final sound.
I perform ifft on array
I alternate real parts of left and right channels (returned by ifft) to create the final audio

The problem is that after this process the sound is distorted. It sounds like the speakers were not plugged in correctly. I found that if I divide values returned by ifft by arbitrary constant then the final sound is right, but is much quieter. I make the division in time domain, on the results from ifft.

The problem doesn't occur if I multiply frequencies by a number less than 1. So if frequencies are attenuated no further division in time domain is needed.

I suppose there is a mistake in the whole process. But if all steps are fine, how should I deal with distorted sound? Is dividing in time domain a proper solution? What number should I use to divide the results then so the sound is not distorted?

EDIT

This is the code I use to perform presented steps. I use Apache Commons math implementation of FFT and SimpleAudioConversion class taken from there http://stackoverflow.com/a/26824664/2891664

Unfortunately, code highlighting doesn't work.

// read file and store playable content in byte array
File file = new File("/home/kamil/Downloads/Glory.wav");
AudioInputStream in = AudioSystem.getAudioInputStream(file);
AudioFormat fmt = in.getFormat();
byte[] bytes = new byte[in.available()];
int result = in.read(bytes);

// convert bytes to float array
float[] samples = new float[bytes.length * 8 / fmt.getSampleSizeInBits()];
int validSamples = SimpleAudioConversion.decode(bytes, samples, result, fmt);

// find nearest power of 2 to zero-pad array in order to use fft
int power = 0;
while (Math.pow(2, power) < samples.length / 2)
    power++;

// divide data into left and right channels
double[][] left = new double[2][(int) Math.pow(2, power)];
double[][] right = new double[2][(int) Math.pow(2, power)];

for (int i = 0; i < samples.length / 2; i++) {
    left[0][i] = samples[2 * i];
    right[0][i] = samples[2 * i + 1];
}

//fft
FastFourierTransformer.transformInPlace(left, DftNormalization.STANDARD, TransformType.FORWARD);
FastFourierTransformer.transformInPlace(right, DftNormalization.STANDARD, TransformType.FORWARD);

// here I amplify the 0-4kHz frequencies by 12dB
// 0-4kHz is 1/5 of whole spectrum, and since there are negative frequencies in the array
// I iterate over 1/10 and multiply frequencies on both sides of the array
for (int i = 1; i < left[0].length / 10; i++) {
    double factor = 3.981d; // ratio = 10^(12dB/20)
    //positive frequencies 0-4kHz
    left[0][i] *= factor;
    right[0][i] *= factor;
    left[1][i] *= factor;
    right[1][i] *= factor;

    // negative frequencies 0-4kHz
    left[0][left[0].length - i] *= factor;
    right[0][left[0].length - i] *= factor;
    left[1][left[0].length - i] *= factor;
    right[1][left[0].length - i] *= factor;
}

//ifft
FastFourierTransformer.transformInPlace(left, DftNormalization.STANDARD, TransformType.INVERSE);
FastFourierTransformer.transformInPlace(right, DftNormalization.STANDARD, TransformType.INVERSE);

// put left and right channel into array
float[] samples2 = new float[(left[0].length) * 2];
for (int i = 0; i < samples2.length / 2; i++) {
    samples2[2 * i] = (float) left[0][i];
    samples2[2 * i + 1] = (float) right[0][i];
}

// convert back to byte array which can be played
byte[] bytes2 = new byte[bytes.length];
int validBytes = SimpleAudioConversion.encode(samples2, bytes2, validSamples, fmt);

You may listen to the sound here
https://vocaroo.com/i/s095uOJZiewf

Chance are your output is clipping. Most commercial wave files are highly compressed (in terms of dynamics) and any type of processing is likely to cause clipping, unless you remove a lot of gain or apply a suitable limiter.
Processing directly in the frequency domain is difficult and likely to result in time domain aliasing. That seems like a poor choice for an audio equalizer
Your band selection is unusual. Most audio equalizer use log-spaced frequency bands to better match human auditory perception. You'll find, that there is almost no energy at all in the 15k-20k band.


Phase difference measurement of a 100 kHz tone I would like to measure the phase difference of a 100 kHz tone. The tone is transmitted and then received.

The received tone is sampled using an 1 MHz 12 bit ADC(IQ) and the SNR is at least 35 dB.

Our aim is to measure the phase difference in at least 360/500= 0.72 degree steps.

Is this possible?
Should I use a PLL or a FFT based approach?

The achievable rms accuracy is approximately the reciprocal of the SNR as an rms quantity, specifically  $10^{-35/20}$ in your case assuming the small angle criteria, if your noise is all phase noise. It is likely your SNR is equally phase and amplitude noise, so the phase result as limited by your noise would be 3 dB better after hard limiting the signal to eliminate the AM component.

I linked other posts that include multiple variants for phase detector implementations, depending on if you have a complex or real representation of the signal. The simple approach is to hard limit the signals (as most phase detectors are sensitive to AM) and then use a multiplier. The output of a linear multiplier is proportional to the cosine of the phase angle  between the two inputs.

Phase Detectors: Phase synchronization in BPSK

Phase Detection for complex signals: How to correct the phase offset for QPSK I-Q data

Phase Detection for Carrier Recovery: High modulation index PSK - carrier recovery


compute a difference( dB) between two BER curves For example, given these 4 curves and I would like to compute a loss/difference between BER of BPSK and BER 8PSk in dB.
How to do that? let say for BER 10^-2

For the segment you show, there's no single number that describes what you want – your lines need to be parallel for that to become a single number.
So, your question can't be answered in general from the diagram, as you did, you'll have to take a fixed BER; and then, I don't really understand the question – that's a horizontal line in your plot. You can read the difference directly from where the curves hits that line.
You can, however, talk about an asymptotic gain if you are able to figure out how these curves continue for large SNRs (you'll have to run with sufficiently many points at say around 23 dB Eb/N0, in my experience, for 16-PSK); at some point, these curves become sufficiently parrallel. Then you can measure the difference in Eb/N0 necessary to achieve the same BER. (that's a horizontal line through the graph; you measure the distance of intersection of curves with the horizontal line.)
Note that this does not work the same way for BER achieved with a given Eb/N0. (i.e. you can't just measure the vertical distance, that doesn't converge as well.)


$N$ point moving average filters in state space I am new to this filter, I did read about them, but could find out a state space notation of these:

$$y(n)=\frac{1}{N}\sum_{m=0}^{N-1}x(n-m)$$

Are moving average filters an LTI systems?
And how do I represent them in state space representation? i.e with $A,B,C,D$ matrices.

Yes sure they are LTI. Let $A$ be the $(L-1)\times (L-1)$ shift matrix
$$
A := \begin{pmatrix}0 & 1 & 0 && \dots & 0\\0 & 0 & 1 & 0 &\dots &0\\\vdots &&& \ddots && 0\\0&&\dots&&1&0\\0 &&\dots&&0&1\\0&0&\dots&0&0&0\end{pmatrix}
$$
and
$$
B = \begin{pmatrix}0\\0\\\vdots\\0\\1\end{pmatrix},\qquad C=\begin{pmatrix}1/L&1/L&\dots&1/L\end{pmatrix},\quad D=1/L
$$
Call $z$ the state variables, then a state space model is
\begin{align*}z(k+1) &= Az(k) + Bx(k)\\y(k) &= Cz(k)+Dx(k)\end{align*}
Infact one has
\begin{align*}
z_1(k) &= z_2(k-1) = z_3(k-2) = ... = z_{L-1}(k-L) = x(k-(L-1))\\
z_2(k) &= x(k-(L-2))\\
&\vdots\\
z_i(k) &= x(k-(L-i))\\
&\vdots\\
z_{L-1} &= x(k-1)
\end{align*}
thus
\begin{align*}
y(k) &= \dfrac{1}{L}\left(\sum_{i=1}^{L-1} x(k-L+i) + x(k)\right)= \dfrac{1}{L}\left(\sum_{m=1}^{L-1} x(k-m) + x(k)\right) \\&= \dfrac{1}{L}\sum_{m=0}^{L-1}x(k-m)
\end{align*}


Single/ Double sided spectrum plot After performing the fft transform there are 2 options to plot the spectrum

the complex (Real and imaginair) spectrum and the magnitude with the phase.
(I neglect the magnitude only for this question)
My question is why is the spectrum of magnitude with phase always single sided.
and why is the complex spectrum always double sided.

second questsion is
make it sense to create a double sided magnitude phase plot or an single sided complex plot?

For real-valued sequences, the DFT is always conjugate symmetrical, i.e.

$$X[k]=X^*[N-k],\qquad k=0,1,\ldots N-1\tag{1}$$

where $N$ is the DFT length. From $(1)$ with $X[k]=X_R[k]+jX_I[k]$ it follows that the real part of $X[k]$ is symmetrical, and the imaginary part is anti-symmetrical:

$$X_R[k]=X_R[N-k]\\
X_I[k]=-X_I[N-k]$$

Furthermore, with $X[k]=M[k]e^{j\phi[k]}$ we get for magnitude and phase

$$M[k]=M[N-k]\\
\phi[k]=-\phi[N-k]$$

So for real-valued sequences, you get the above symmetries in the real and imaginary parts, and also in the magnitude and phase. So no matter if you represent the DFT in terms of real and imaginary parts, or in terms of magnitude and phase, a one-sided spectrum (i.e. for indices $0,\ldots, \lfloor N/2\rfloor$) is always sufficient, as long as the sequence is real-valued.


Are there any quick way to tell if a kernel is a high-pass filter? Im new in this image processing field so please be gentle with terminology, i have a kernel of 3x3 and i want some way in order to test if its a high pass filter or not

Thanks in advance

Proper way: do a Fourier transform and look at the magnitude of the result
Quick and dirty: take the sum of all all elements. If it's zero or small compared to the absolute sum, it has some high pass-ness to it.


Is this an infinite response or a finite response? If you're given $h[n] = 1$ for every even number of $n$, otherwise it'll be 0. Is that considered an infinite response or a finite response and how do you know?

A finite impulse response (FIR) filter has an impulse response of finite length. It can be implemented by a non-recursive structure, but there are also cases where the designer would choose a recursive structure for implementing an FIR filter (e.g., a recursively implemented moving average filter, as explained in this answer).

An infinite impulse response (IIR) filter has an impulse response of infinite extension. Consequently, it cannot be implemented by a non-recursive structure, but it requires a recursive structure.

With the above explanation in mind, I think it should be straightforward to decide which class the system in your question belongs to.


WFM SDR and Audio resolution How does SDR resolution translate to audio resolution when using FM demodulation?

For example: using a 8-bit/2.048 Mhz IQ stream, what is the maximum audio sample rate and bit depth that can be generated from a +/-75kHz mono WFM signal with 20kHz audio?

My understanding is that the SDR bit depth is irrelevant but for link quality purposes. But how does the SDR sample rate translate to sample rate and bit depth of the resulting audio signal?

how does the SDR sample rate translate to sample rate and bit depth of the resulting audio signal?

Extremely indirectly. I haven't got precise mathematical formulas for you, but here's an overview of the topic which should let you figure out where you want to dig deeper.

First of all, FM is an analog modulation. This means that there is no inherent sample rate or bit depth of the signal recovered by demodulating, but only an audio bandwidth and noise floor. So let's look at what influences that:

The bits per sample (bit depth) of the digitized RF determines the amount of error introduced by the quantization (conversion to discrete samples), which can be described as quantization noise.

Note that quantization noise is at a constant level in the digital signal coming from the ADC. This means that — unlike radio interference — the signal-to-quantization-noise ratio improves the more you turn up your receiver's gain. But you can only do that so far before you get clipping (signal level exceeding the representable sample values).

Therefore, the best results are obtained when you have a narrow analog filter in your receiver to remove all but the desired signal before digitization — the wider bandwidth you receive, the more total power is arriving, and so the less gain you can use before clipping.

Of course, using a narrow bandwidth prevents you from having those pretty waterfall displays of ten stations at once. It's a tradeoff.

For more words on the subject, see my previous answer over on Amateur Radio Stack Exchange, to the question Why is dynamic range relevant for an SDR?
A simple analysis would suggest that the sample rate of the digitized RF has no effect, because it's just bringing in a wider bandwidth, full of unrelated signals and noise, than we actually need. However, this is not true, if the ADC has fewer bits per sample than the actual DSP processing (for example, 8-bit integer samples from the ADC going into 32-bit floating-point DSP operations). The extra samples contain some of the information that was lost in quantization; resampling the unnecessarily-high-rate signal to a lower rate will automatically use that information, so those 32 bits aren't actually wasted.

So to combine this with the previous point, the best case given fixed hardware is to configure it (if possible) for the narrowest analog filter bandwidth and the highest sample rate (creating a signal whose spectrum has an obvious filter hump in the middle) and then resampling it before demodulation.

To sum up, the sample rate and bit depth of the incoming RF determines how much noise is added to the signal by the process of digitization. The next question is how that noise, the RF noise entering the software demodulator, affects the noise level of the resulting audio, and I'm afraid I don't know that and haven't found an obvious reference.

Now that we know that our input can be considered as a signal with more or less noise depending on its parameters, let's look at what to do with the output.

As Nyquist tells us, the sample rate of the output audio signal does not need to be greater than twice the highest frequency present in the signal. That's the audio bandwidth, which for broadcast FM is, let's say, 15 kHz. Therefore, the required sample rate is 30 kHz, ideally. However, the signal must be filtered to remove noise that is outside the expected bandwidth, and using exactly 30 kHz sample rate means that the noise which passes the necessarily-imperfect filter will be aliased down into lower frequencies. Thus, it is better to use a higher sample rate to allow some room for the filter transition band — perhaps 44.1 or 48 kHz, because those are common audio sample rates.
Quantization introduces quantization noise — but we already have noise in our signal. Therefore, the bit depth of the audio should be chosen so that it introduces an insignificant amount of quantization noise compared to the noise already present. (Note: You can't avoid introducing any noise at all. Your finite-precision arithmetic is doing that at every step.)


SNR for complex valued signal and when to use SNR and when to use Eb/N0? This question is in continuation to the one asked earlier Help in understanding the formula of Signal-to-Noise-Ratio (SNR) - Part 1

Consider a model
$$y[n] = x[n] + w[n]$$

Signal means the desired information which in this case is x and y is the output noisy observation.

Question A

Case 1: When noise is AWGN $ \sim C\mathcal N(0,2\sigma^2_w)$; $y$ and $x$ is also complex valued i.e.,
$w \sim CN(0,2\sigma^2_w), x \sim CN(0,2\sigma^2_x)$

Then would SNR = $\frac{E[x^2]}{2\sigma^2_w} = \frac{2\sigma^2_x}{2\sigma^2_w} $ ?

Case 2: AWGN $ \sim \mathcal N(0,\sigma^2_w)$; $y$ and $x$ are real valued. Then would SNR = $\frac{E[x^2]}{\sigma^2_w}$?

These topics are not clearly mentioned in the Answers given for the other Question. The Answers explain what is SNR but my confusion remains that if the signal and the noise are complex valued, then would there be any change in the SNR due to the real and imaginary component? In the denominator of the formula would there be a 2? That is why I have asked this here.

Question B: https://en.wikipedia.org/wiki/Signal-to-noise_ratio shows 2 formulae for SNR. In many research articles, the formula of SNR is not used, instead EbN0 is used. I am using Least squares estimator to estimate channel coefficients and then trying to plot the mean square error (in db) of the channel estimates with varying SNR. NExt, I want to plot the bit-error rate (using simulation) for bpsk and say QAM modulation. I used SNR for the plot of mean square error. So, should I use SNR for Bit error rate? OR is there no such rule that both the plots should be using SNR and not EbNo formula?

Question A:
For the considered system model
$$ y[n] = x[n] + w[n], $$
as already pointed out in your previous question, it holds:
$$ E\{ y[n]y[n]^* \} = E\{ (x[n]+w[n])(x[n]+w[n])^* \} \\= E\{ x[n]x[n]^* \} + E\{ x[n] \}\underbrace{E\{ w[n]^* \}}_{= 0} + \underbrace{E\{ w[n] \}}_{= 0}E\{ x[n]^* \} + E\{ w[n]w[n]^* \} \\= E\{ x[n]x[n]^* \} + E\{ w[n]w[n]^* \}.$$
As a result, your $SNR$ will always be$$ SNR = \frac{E\{ x[n]x[n]^* \}}{E\{ w[n]w[n]^*\}}. $$
So for both of your cases the $SNR$ you specify is correct. It just depends on your definition of $x[n]$ and $w[n]$.

Question B:
For your application I would recommend to use $\frac{E_b}{N_0}$, because the energy per bit depends on the modulation scheme (BPSK, QAM etc.). Therefore comparing the $BER$ in terms of $\frac{E_b}{N_0}$ will be more fair.

You can find further information:

https://en.wikipedia.org/wiki/Eb/N0

https://de.mathworks.com/help/comm/ug/awgn-channel.html

http://www.gaussianwaves.com/2008/11/relation-between-ebn0-and-snr-2/


FFT on non-uniformally sampled signals Most fft tutorials use sinusoidal signals for demonstration, which makes the user already know what the fft is supposed to give as an output, but what about real time signals where we don't have any previous knowledge except the sampled data with the following problems:

No previous knowledge
Variable sampling time
Unknown Noise

I have a signal which has the previous characteristics, so I used the following Matlab code to interpolate and select a proper period (visually):

clear all
load('ExperimentData.mat')
close all
% Because the sampling period is not constant, I have first to interpolate the data
Ts=min(diff(t));Fs=1/Ts;
t_new=t(1):Ts:t(end);
x_new=interp1(t,x,t_new);
%% Extracting one period and Applying ZERO PADDING
tStart=0.7086;[~,idxStart]=min(abs(t_new-tStart));
tEnd=0.7686;[~,idxEnd]=min(abs(t_new-tEnd));
x1p=x_new(idxStart:idxEnd);
t1p=t_new(idxStart:idxEnd);
F=figure;plot(t,x,'b');hold on
plot(t_new,x_new,'r:')
plot(t1p,x1p,'.-g');
ylabel('x');xlabel('time[sec]');legend('Original signal','Interpolated Signal','One Period')

The outcome is the following:

I then applied the FFT once on the complet signal an once on one period with ZERO padding using the following code:

X = fft(x_new-mean(x_new));
N = size(x_new,2);
f=Fs/2 * linspace(0,1,N/2+1);
figure;
AX(1)=subplot(211);semilogx(f,abs(X(1:N/2+1)));ylabel('Magnitude');xlabel('Frequency(Hz)');title('Entire signal Without Zero Padding')
AX(2)=subplot(212);semilogx(f,angle(X(1:N/2+1)));ylabel('Phase');xlabel('Frequency(Hz)')
linkaxes(AX,'x');axis  tight;
%Zero Padding
x_1PZP=[x1p-mean(x1p) zeros(1,100000)];
X_1PZP = fft(x_1PZP);
figure;
AX(1)=subplot(211);semilogx(f,abs(X_1PZP(1:N/2+1)));ylabel('Magnitude');xlabel('Frequency(Hz)');title('One Period with Zero Padding')
AX(2)=subplot(212);semilogx(f,angle(X_1PZP(1:N/2+1)));ylabel('Phase');xlabel('Frequency(Hz)')
linkaxes(AX,'x');axis  tight

The results are:

My wonderings are:

Does the interpolation of a signal alter the FFT ?
From the Time Signal, the FFT on the complet signal is closer to reality. but it still doesn't show all frequencies. Why?
I almost can't extract any information fro the phase plot. What am I missing?

The data can be obtained from this link.

"Does the interpolation of a signal alter the FFT?"
Yes, it does.
The key point is that in a system that uses sampling, you simply do not know what does the signal "do" between two sampling instances. For a high enough sampling frequency Fs, it is just assumed that the corresponding sampling period is high enough to be approximated by a series of sinc pulses. For more information please see: http://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem

By using interpolation, you only assume that the signal behaves in a certain way during the time that it is not observed. Therefore, the result of the FFT will be different, even for different interpolation methods.

Furthermore, interpolation has a "natural" smoothing action because all it can do is "degrade" the existing data, it cannot add anything more to it. For more information about this notion, please see the Data Processing Theorem (https://books.google.co.uk/books?id=F4j8EHBKS1IC&lpg=PA313&ots=tHO2C1ZJs1&dq=the%20data%20processing%20theorem&pg=PA313#v=onepage&q&f=false)

This is why, interpolation will tend to represent lower frequencies (for which you have more data anyway, due to the sampling frequency) better.

By the way, 'interp1' by default uses linear interpolation, that is, intermediate values between sampling instances are assumed to lie on a straight line that is defined by the two known points. An alternative to this would be spline interpolation which takes into account the slope by which the interpolated curve should enter and depart from the interpolating area. For more information about this, please see http://en.wikipedia.org/wiki/Cubic_Hermite_spline

(Please also note Trigonometric Interpolation, a technique very closely associated with the Fourier Transform. For more information please see: http://en.wikipedia.org/wiki/Trigonometric_interpolation)

From the Time Signal, the FFT on the complet signal is closer to reality. but it still doesn't show all frequencies. Why?

Because it is not exactly the same signal, as far as the FT is concerned.

A key property of the Discrete Fourier Transform is "Periodicity". According to this, the signal at the input of the Fourier Transform is assumed...periodic. That is, you may be passing a voice recording of "One Two Three" to the FT to recover its harmonic component but the FT assumes that the signal, in reality sounds like "...One Two Three One Two Three One Two Three One Two Three..." extending to infinity both forwards and backwards. For more information about this please see: http://en.wikipedia.org/wiki/Discrete_Fourier_transform (Periodicity property).

The signal that you call "Time Signal", contains a brief period of silence in the beginning and ends abruptly, that is, it would be discontinuous if it were to be a periodic signal. By the way, you may also be interested in the technique of windowing, for more information please see: http://en.wikipedia.org/wiki/Window_function#Symmetry_and_asymmetry

I almost can't extract any information fro the phase plot. What am I missing?

The phase spectrum is simply the angle between the imaginary and real components of the signal spectrum. However, having seen how the imaginary and real components of the FT look like, their phase content is not expected to be a smooth and well behaved function, except perhaps for very simplistic cases.

By the way, you may want to look up the concept of phase unwrapping, to improve your phase plots. For more information, please see https://ccrma.stanford.edu/~jos/fp/Phase_Unwrapping.html

Hope this helps.


Signal-to-Noise Ratio (SNR) needed to discern two superimposed signals Say I am trying to measure a exponentially decaying signal ($S_y$), that constitutes a portion $\gamma$ (about 5% so $\gamma=\frac{1}{20}$) of total signal measured($S=S_y+S_b$ and $S_b=S_0 \beta e^{\frac{-t}{C_b}}$ and $S_y=S_0 \gamma e^{\frac{-t}{C_y}}$). Where $\gamma+\beta=1$ and $C_Y$ is about $\frac{1}{4}$ of $C_b$. I want to determine the least Signal-to-Noise Ratio (SNR) needed to be able to discern that signal. The noise in our signal is additive white noise. Should the SNR be $>\gamma=\frac{1}{20}$ or is there a better more formal way to decide the SNR limit? $$S(t)=S_0(\gamma e^{\frac{-t}{C_y}}+\beta e^{\frac{-t}{C_b}}) + n(t)$$

So my question is how do I determine the minimal SNR when I measure signal $S$ over time $t$, but the thing I am interested in is fluctuation (max 10%) in $S_y$. All constants are known to me, however i only know the mean of $C_y$, as it fluctuates and so $S_y$ does too. Thanks for any help. If I need to clarify any details please tell me and I will do so.

This is a work-in-progress until more information is added to the question.

First, I'm assuming that you know the magnitudes of your decaying exponentials, but you don't know their time constants (either $C_x$ or $C_y$) or your signal level $S_0$.

Second, you say that there is Gaussian noise, but not whether it's additive. I'm going to assume it's additive. That means your signal is really:

$$
S(t)=S_0(\frac{1}{20}e^{\frac{-t}{C_y}}+\frac{19}{20}e^{\frac{-t}{C_x}}) + n(t)
$$

where $n(t)$ is additive, white, Gaussian noise with zero mean and variance $\sigma^2$.

Then, to set up the problem, we need to postulate the estimated values of $\hat{C}_x$, $\hat{C}_y$, $\hat{S}_0$, and $\sigma^2$ in $\hat{S}(t)$:

$$
\hat{S}(t; \hat{C}_x, \hat{C}_y, \hat{S}_0)=\hat{S}_0(\frac{1}{20}e^{\frac{-t}{\hat{C}_y}}+\frac{19}{20}e^{\frac{-t}{\hat{C}_x}})
$$

Then a least squares approach to solving it would define the error, $E$, as:

$$
E(\hat{C}_x, \hat{C}_y, \hat{S}_0) = \int_{I} \left| S(t) - \hat{S(t)}\right|^2 dt
$$

where $I$ is your time interval of interest, and minimize this with respect to $\hat{C}_x$, $\hat{C}_y$,  and $\hat{S}_0$.

I interpret your question as asking: what value $\frac{S_0^2}{\sigma^2}$ does this work over?

The answer is: it depends on what you mean by "work".  You can get an estimate for any noise level.  The question is, what error can you tolerate in that estimate.  Also, it will depend somewhat on your integration length ($I$).

EDIT

OK, I see you've changed notation so it's $C_b$ and $C_y$ now.

I've written a short scilab script to try to see how things change with noise level. I've assumed EVERYTHING is precisely known except the $C_y$ parameter.

The error in estimating the parameter versus the value of $\sigma$ is shown in the plot below.  Code below generates it.

 S0 = 1;
 Cb = 201;
 Cy = 100;
 sigma = 0.1;
 mb = 1/20;
 my = 19/20;

 T=1000;
 t=[0:T-1];
 Sb = mb*exp(-t/Cb);
 Sy = my*exp(-t/Cy);

 clf
 subplot(311)
 plot(Sb)
 plot(Sy,'g')
 plot(S,'r')

 ERRORS = [];
 sigmaRange = [0.0 0.1 0.2 0.5 1 2 5 10];
 for sigma=sigmaRange,
 CyhatRange = 90:.01:110;
 estimates=[]
 NRuns = 100;
 for n=1:NRuns
     S = S0*(Sb + Sy) + sigma*rand(1,T,'normal');

     ERR= [];
     for Cyhat = CyhatRange,
         S_hat = S0*(Sb + my*exp(-t/Cyhat))

         ERR = [ERR; sum((S-S_hat).^2)];
     end

     [mx,ix] = min(ERR);

     CyhatEst = CyhatRange(ix);
     //disp(CyhatEst);
     estimates = [estimates; CyhatEst];

     subplot(312);
     plot(CyhatRange,ERR)
 end

 disp(mean((estimates-Cy).^2))

 ERRORS = [ERRORS; mean((estimates-Cy).^2)];
 end

 subplot(313)
 plot(sigmaRange,ERRORS);


Given 1 recording of my voice, how can I process it to sound like 2 independent recordings of my voice played in unison? I have one recording of my voice $x[n]$ singing a song. One thing I could do is independently record me singing the song again to obtain $y[n]$. Then, $z[n] = \frac{1}{2}\left(x[n] + y[n]\right)$ would sound like a "chorus" of some sorts (assuming that both recordings were done in perfect sync).

My question is, without doing a bunch of independent recordings, is there some decent way to go straight from a single recording $x[n]$ to the 2-fold chorus $z[n]$? More generally, I'd like to be able to go to a $k$-fold chorus (equivalent to averaging $k$ independent recordings).

My gut tells me one could simulate independent recordings by randomizing the phase in $x[n]$'s STFT somehow, but I'm not sure of specifically what should be done. In this regard, the only relevant method I'm aware of is Paul stretching, which doesn't really fit the bill, I think.

Also, is there a standard name for this particular effect I'm trying to create?

To make it sound natural you typically introduce some small random pitch, speed & gain modulations. That's not a trivial amount of work, especially if you want something that sounds natural and good and maintains the original phrasing.

This is a pretty common plug-in in audio processing. It's typically called "vocal doubler" or "voice doubler" and there are many off-the-shelf options available. Here is an example: https://www.izotope.com/en/products/vocal-doubler.html. No endorsement or commercial interest, but you can't beat the price :-)


Poisson noise and curve fitting - denoise first? If I have an image that is severely corrupted by Poisson noise, and I want to fit a function to the image, is it "better" to attempt to denoise the signal first before fitting, or should I move straight to the fitting stage?

In the example below, a 2D Gaussian function has been corrupted by Poisson noise. Should I fit a 2D Gaussian to the noisy data, or to a denoised version?

Denoising images is often good for qualitative reasons, but I'm curious to know about the quantitative case, for example where the volume of the Gaussian is important?

By "denoising", I'm thinking along the lines of the techniques such as Non-local PCA rather than median filtering etc.

If the example images you've given are at all representative of your application, you may want to consider thinking about the problem a little differently. Instead of thinking of the image as "corrupted by Poisson noise", think of the observed data as a limited number of photons sampled at random from the latent image intensity map. The photon counts you get are providing incomplete information about the latent image, not corrupt or noisy information.

From this perspective, if you are truly fitting a curve or formula of few parameters to your data, this is a classical parametric statistical estimation problem (lots of samples, few parameters). Such problems have a good estimation theory which guarantees you will accurately learn the fit parameters without needing too many sample points. If you apply a denoising step before parametric fitting, the theoretical guarantees will no longer apply because the error in the denoised image is some complicated composite of residual fit variance and bias from the denoiser's implicit signal model.

So there's a theoretical justification for fitting the function straight to the data. Just make sure that you use the correct Poissonian likelihood function when you do so, or use the EM algorithm if you don't want to worry about such things.

EDIT: actually you don't even need Poissonian likelihood here (unless your data is quantized/binned). Just take the approach described in http://en.wikipedia.org/wiki/Maximum_likelihood. Here the $x_1, x_2, ... x_n$ are a collection of $n$ vectors giving the 2-D coordinates of your photon observations, and $\theta$ is the parameter vector for the function you want to fit. In the case of a Gaussian you have $\theta = (\mu,\Sigma)$ where $\mu$ is the mean (centroid in 2D) and $\Sigma$ is the covariance matrix. Then in the notation of the wiki article, your likelihood is

$$ \prod_{i=1}^n f(x_i | \theta) = C \prod_{i=1}^n \exp( -\tfrac{1}{2} (x_i - \mu)^T \Sigma^{-1} (x_i - \mu)) $$

where $C$ is some constant with square roots of pi in it and such :) If your data is quantized/binned then you consider the count value in each bin as a Poisson random variable and go from there. Weirdly, the Poisson and the $n$ observation approaches are nearly equivalent from an optimization standpoint.

I don't have a great reference for you offhand because I learned this stuff in lab classes or wherever years ago. But the wiki article and refs therein seem good. This paper seems ok too: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.74.671&rep=rep1&type=pdf


How to display the output as the input with a fractional delay? I am trying to introduce a fractional delay in my input vector x[n] so the output is something like y[n]=x[n-M/N]. I am not supposed to use the MATLAB fracdelay function.So, I used sampling and delay. The three steps I followed are:

1)Upsample the input by M.

2)Delay the upsampled version.

3)Downsample the output obtained in step 2 to get y[n].

However I can't get the correct fractionally delayed output. Can someone please help? Here is my MATLAB code:

x=[1 2 3 4 1 3 4 5 6 3];
up=upsample(x,2);

z=zeros(1,10);

for n=1:10;
   k=n+1;
   z(1,k)=up(1,n);
end

y=downsample(z,2,1);
subplot(411);
stem(x)
subplot(412)
stem(up)
subplot(413);
stem(z)
subplot(414);
stem(y)

Here is how my output looks like:

The problem is that you are using "upsample", which merely inserts zeros in between the original samples, instead of "interp", which upsamples and then low-pass filters.


Convert measured spectrum to impulse response I am still learning, so please bear with me:
My audio measurement tool exports the measured spectrum of the DUT in the following format:

Frequency [Hz]
Magnitude [dB]
Phase [degree]

It is stored in a simple text file.
I managed already to read the file in Matlab, split the lines into vectors for each column. So I have now a vector for the frequency points, a vector for the magnitude values and a vector for the phases.
Now I want to generate the impulse response from it. I managed already to interpolate the values on a FFT frequency grid e.g. for a 4096-points FFT and that the second half (the negative frequencies) can be obtained with conj(X) and mirrored. That's all clear. So I can run ifft on it.
What is not clear to me is, how to handle the phases.
Can I just write:
X = 10.^(mag/20) .* exp(j*phi/360*2*pi)

with

mag the vector with imported magnitude values in dB and
phi the vector with imported phase values in degree?

Or do I need to apply another conversion to phi? Something that takes the sampling frequency  into account?

Can I just write: X = 10.^(mag/20) .* exp(jphi/3602*pi) ?

Yes. If your software returns the phase in degrees, that's exactly right.
The tricky part is, however, the interpolation. The phase wraps round at 360 degrees. For example, if you have 1 degree and 355 degree, the correct interpolation would 0 degree (since 355 is the same as -1) but regular interpolation will give you 180 degrees which is about as wrong as you can get.
In this case unwrapping the phase often helps. You would start somewhere in the middle of the frequency spectrum where you have good signal to noise ratio and than unwrap up and down from there. Unwrapping means to look for phase jumps of more than 180 degrees and add/subtract multiple of 360 degrees to the phase until the jump is less than 180 degrees.
Another potential problem is bulk delay in the measurement which is fairly common for acoustic measurements. Bulk delay leads to very steep phase gradients which makes both unwrapping and interpolation difficult. This can be addressed by "eye balling" the bulk delay, take it out of the measurement. Then do the interpolation and inverse FFT and put the bulk delay back into the impulse response.


Finding the friction coefficient by using experimental data I have the differential system equation with some unknown parameters (friction coefficients)
Also I have experimental data. What I want is to determine the friction coefficients which best fit the experimental data. How to do it using the Mathematica, I know there NDSolve and FindFit functions, but they can't be used here? And what are the general approaches of this problem?

Edit 1:
Here is the model:

$$k_1 (x_3[t]-x_1[t])+a_1 \left(x_3'[t]-x_1'[t] \right)+M_1 x_1''[t]=0$$
$$k_2 (x_3[t]-x_2[t])+a_2 \left(x_3'[t]-x_2'[t]\right)+M_2 x_2''[t]=0$$

$$\begin{eqnarray}
& k_3 x_3[t]-k_1 (x_3[t]-x_1[t])-k_2 (x_3[t]-x_2[t])
\\& + a_3 x_3'[t]-a_1 \left(x_3'[t]-x_1'[t]\right)-a_2 \left(x_3'[t]-x_2'[t]\right)
\\& +C_1 x_3''[t]=0
\end{eqnarray}$$

I can only propose a crude starting point based on a more modest work I did. My goal was to determine the bleaching coefficient of a fluorophore when submitted to varying illumination intensity over time. I had the theoretical law for the bleaching over time as a function of intensity (depending on unknown parameters), the amount of fluorophore over time, and the starting intensity and intensity at the end of the experiment. The goal was to retrieve the parameters for the bleaching law and the whole intensity over time.

Anyway. What I did was to construct (using MATLAB) a function to minimize. In your case, given that you know $x(t)$ experimentally, I suggest to build a function $O(x, y)$ such that:

$$ O(x, y) = \sum_t |x(t) - y(t)|^2 $$

with $x(t)$ being your experimental data and $y(t)$ being the solution of your differential system for some parameters $k, M, a$. Put this function as the objective in a simplex optimizer that will find its minimum, allowing $k, M, a$ to vary.

This is done through numerical solution. So you have a function to minimize, say $F(k,a,M)$ that returns the value of $O(x,y)$ (which is scalar btw) for the current numerical resolution with the given $k,a,M$.

At each step, the whole system will be solved again, and the solution will be matched against the experimental data until they fit.

I have found this to be working quite well in my case, as long as you are sure you have the correct model to describe the experimental data. Otherwise, expect some seriously irrelevant results.


How to implement an improper discrete-time transfer function into Simulink I am trying to model a discrete-time control system in Simulink. The control system contains an inner loop that utilizes the inverse of the plant model. The plant in question takes the form of:

$P(s) = \frac{3.332\times s + 4.679}{s^2+30.97\times s+45.5}$

s = tf('s');
P = (3.332*s + 4.679)/(s^2+30.97*s+45.5);

I inverted the plant, $P(s)^{-1}$:

P_inv = inv(P);

Now I have an inverted model however you cannot explicitly implement an improper function into Simulink. I followed this discussion that has been referred to multiple times scattered throughout forums.

[r,p,k]=residue(P_inv.num{1},P_inv.den{1});
sys1 = tf(r(1),[1,-p(1)]);
sysk1 = k(1)*s;
sysk2 = k(2);

Once the improper function is broken down, I discretized each component separately.

Ts = 0.005;
d_sys1 = c2d(sys1,Ts)
d_sysk1 = c2d(sysk1,Ts,'matched') % <--- problem

A problem now exists in d_sysk1. It results in yet another improper function taking the form of:

d_sysk1 =
    59.99 z - 59.99
Sample time: 0.005 seconds
Discrete-time transfer function.

In my simulation I want to implement a discretized improper transfer function, P(s). Is there something I did incorrectly?

I think you missed the point in the link you posted. Implement the improper function as a differentiator. You cannot get rid of the proper part but only replace it with differentiators. Hope this helps!


Correlation between Sin and Square wave I would like to know how $\frac T2$ came for time period in the image below.Because I think the frequncy is 1 so time period should be 1.Is that correct? or can anyone help me..

The article from where I got this image is here

Note that the second half of the period of the sinusoid is simply equal to the first half times $-1$. Since the same is true for the square wave, the multiplication of the two is the same in both intervals $[0,T/2]$ and $[T/2,T]$. That's why the integral over the whole period equals twice the integral over the first half of the period.


Low analog lowpass + ADC vs High analog lowpass + ADC + low digital lowpass I'm doing EEG data collection and my amp forces use of an analog lowpass filter before the ADC. I'm fairly sure that I won't be interested in phenomena associated with frequencies of 30Hz or higher, so it strikes me that the most disk-space efficient thing to do is to analog lowpass at 30Hz (it uses a Bessel filter) and tell the ADC to sample at 60Hz. However, disk space considerations aside, I'm curious whether there is any difference between this strategy and, say, setting the lowpass to something really high like 250Hz and sampling at 500Hz, then once digitized and doing analysis, applying a digital lowpass at 30Hz. Would the double filtering induce artifacts? (i.e. Should I opt for the low analog lowpass for reasons other than disk space?)

I'm curious whether there is any difference between this strategy and, say, setting the lowpass to something really high like 250Hz and sampling at 500Hz, then once digitized and doing analysis, applying a digital lowpass at 30Hz

Good idea.  This is how oversampling ADCs work, sampling at a high rate using a poor quality converter (even a 1-bit converter), and then using digital filtering to downsample to the desired value with better accuracy.  Digital filters can have a much sharper cutoff than analog filters for essentially no cost, reducing aliasing that would otherwise get into your recording.  It can also improve your signal-to-noise if the analog signal isn't already restricting this.

  Would the double filtering induce artifacts?

Not unless you design the filters poorly.  :)


Difference equation when transfer function expressed as poles and zeros The transfer function $H(z)$:

$$
H(z) = \frac{Y(z)}{X(z)} = \frac {b_0 + b_1 z^{-1} + b_2 z^{-2}} {1 + a_1 z^{-1} + a_2 z^{-2}}
\tag{1}
$$

Has difference equation:

$$
y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] - a_1 y[n-1] - a_2 y[n-2]
\tag{2}
$$

How would you compute the difference equation when $H(z)$ is expressed as poles and zeros without expanding brackets:

$$
H(z) = g \frac {(z - d_0)(z - d_1)} {(z - c_0)(z - c_1)}
\tag{3}
$$

I'm happy to assume real poles and zeros for the sake of simplicity. I know that when performing a cascade on a computer that the output of the first section is fed to the second and second order sections are used. I'm interested in the maths though.

What I'm ultimately trying to do is get the difference equation for this expression direct without expansion:

$$
H(z) = g \prod_{k=1}^{K} \frac {b_{k0} + b_{k1} z^{-1} + b_{k2} z^{-2}} {1 + a_{k1} z^{-1} + a_{k2} z^{-2}}
\tag{4}
$$

If I understood your question right, you would like to obtain a Linear Constant Coefficient Difference Equation (LCCDE) representation of a given system from its given System Transfer Function, $H(z)$, when it's expressed in a pole-zero product form such as:
$$
H(z) = g \frac {(z - d_0)(z - d_1)} {(z - c_0)(z - c_1)}
\tag{1}
$$

instead of a, more directly apparent, power series form such as:
$$
H(z) = \frac {b_0 + b_1 z^{-1} + b_2 z^{-2}} {1 + a_1 z^{-1} + a_2 z^{-2}}
\tag{2}
$$

As known, the latter form can be directly converted to its corresponding LCCDE representation from the coefficients ${b_k}$ and ${a_k}$ as follows:
$$y[n] + a_1 y[n-1] + a_2 y[n-2] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2]$$
Where the pattern of conversion speaks for itself, without a need of further description.

Now, we'd like to know whether one can obtain a similar conversion from the pole-zero form into the LCCDE without expanding the product?

Well yes and no! Depending on how you consider the process of computing the coefficents. If your intention is on the fact that you do not want to algebricaly perform the paranthesis expansion operation, but instead use a machanism that computes each coefficient separately by a set of arithmetic operations; then yes you can devise such a method to, for example, compute $a_3$ or $b_5$ without computing any other terms. And No, eventually you must always compute the coefficient $a_k$ that multiplies $z^{-k}$ so as to find the coefficient that multiplies $y[n-k]$, hence you need to always expand the brackets and compute multipliers of $z^{-k}$ in one or the other way.

Then I want to present here a simple method to find those coefficients :
given a polynomial in the product of zeros form: $$p(x) = (x-a)(x-b)(x-c)$$ we can perform the following operation to get its coefficients in the expanded form $$p(x) = d_0 x^3 + d_1 x^2 + d_2 x + d_3$$

$$d[k] = (1 -a)\star(1 -b)\star(1-c)$$ where $k =0,1,2,3$ in this particular case.

As clearly seen this is a discrete convolution operation, where each multiplier bracket is turned into a respective convolution operand.  For which a simple MATLAB line to compute $d[k]$ would be:

 d = conv([1 -a],conv([1 -b],[1 -c]));

where a, b, and c would be replaced by their numeric values or if you're competent in the symbolic math capability, you could instead obtain a pure symbolic result for the coeffients $d_k$, which would therefore give you a way to compute LCCDE coefficients from pole-zero products of $H(z)$ in a rather roundabout way.

Finally, the computations for a power of $z^{-1}$ products is basicly the same if you describe the product as $$P(z) = (1 - az^{-1})(1-bz^{-1})(1-cz^{-1})$$ and,
$$P(z) = d_0 + d_1 z^{-1} + d_2 z^{-2} + d_3 z^{-3} $$

where again $$d[k] = (1 -a)\star(1 -b)\star(1-c)$$


Physical meaning of average values of random signals This question might be a bit stupid, anyway, i'll risk it, since i want to get better understanding of this subject.

Let's consider random signal x(t), and let's say that we know that it is ergodic by mean, and by autocorrelation (this implies that it is wide sense stationary).

This means that time averaged mean of this signal is equal to it's statistic averaged mean, same goes for autocorrelation.

So we have  the following:

$\mu_X=E(x)=\frac{1}{T} \int_T x(t)dt=\langle x(t)\rangle  \\ R_{xx}(\tau)=E[x(t)x(t+\tau)]= \lim_{T->\infty} \frac{1}{T} \int_T x(t)x(t+\tau)dt=\langle x(t)x(t+\tau)\rangle$

Then, as i found in one of the books i have, it says following:

$\langle x(t) \rangle$- corresponds to DC level of given signal

$\langle x(t) \rangle ^2$ - corresponds to normalized power of DC component

$\langle x^2(t) \rangle$ - corresponds to total average normalized power

$\overline{\sigma_x^2}= \langle x^2(t) \rangle - \langle x(t) \rangle^2$-corresponds to average normalized power in AC component of the signal

$\overline{\sigma_x}$-corresponds to rms of AC component

My question is, is there any way to prove this, or to gain intuitive understanding why this is true. Any help appreciated!

Dilip: I'm confident that what you said is correct but the weirdo sample path that you refer to ( the one with probability zero ) can be used to estimate the autocorrelation at different lags, namely $\hat{R}_{xx}(\tau)$. This is true even though it's only one realization and this is the case only because of the ergodicity assumption. So, to a person coming in new, my opinion is that your answer doesn't give enough significance to the one realization. It's a weirdo but a useful weirdo for sure  !!!!!


how to process audio files in Matlab I have project in which I'll work on audio files using Matlab, I'll be listing to different  machine and should get thei rotational speeds, which is kind of periodic, my question how should I analyse this should I  work with the spectrum and its magnitude using the fft or should I use the PSD , which is the result of pwelch.
I'm not asking for code! thanks for your help

Spectrum, FFT magnitude, and pwelch are pretty much the same thing, i.e. a metric of energy density as a function of frequency. They all can be configured in different ways: window length, window shape, overlap, FFT size, lin/log representation, smoothing, different averaging methods, etc.
All these parameters need to be chosen carefully depending on the requirements of your application: how fast do you need to track, what frequency resolution do you need,  what's the  noise spectrum look like etc.
Choosing the parameters properly is actually more important than the core algorithm (which is more or less the same anyway).


Conversion from Floating point representation to Fixed point representation I have written a matlab function for a constellation mapper and I want to create a Fixed-point version of that code for implementation on FPGA. I know matlab has a specific toolbox for fixed-point conversion but i want to do this manually. My specific question is:

1) which parts of the code would require this conversion? i mean do all inputs, outputs, temporary  variables need to be converted to a fixed point format ?

2) What are the steps that I should follow to do this conversion? i read that converting from floating point to fixed is no more, no less than scaling numbers. but I'm not quite sure how I go about doing this. Also, I want this fixed point representation to be of a custom length, not necessarily equal to standard lengths like 8, 16, 32-bits .

PS: I don't want to compare the performance of the fixed point code and floating point code in the time being

I think you need to use the floor() or round() functions in MATLAB, to emulate fixed-point variables and operations.  So you have to know the range of the fixed-point value and the precision of it.  The ratio of the range to the precision is the dynamic range and you get 6.02 dB and one bit of word width every doubling of that ratio.

If you convert:

x_fixed = precision*round(x_float/precision);
if x_fixed > fixed_max
    x_fixed = fixed_max;
elseif x_fixed < fixed_min
    x_fixed = fixed_min;
end;

then fixed_max - fixed_min = (2^N - 1)*precision  for N being the number of bits in your fixed-point word.  If you do this after every mathematical operation, you will be emulating fixed-point arithmetic.

there might be better ways to round x_float than i shown.  sometimes the floor() is better than round().  and "fraction saving" is a better way to round an intermediate result with more bits to your N-bit fixed-point word, but requires an error state for each variable.

for every variable and do your arithmetic with their fixed-point values and round the result to a fixed-point value with code like above.

that's how to convert floating-point MATLAB to emulating fixed point math.