PSD of modulated signal I know this question has been previously asked just saw it and tried to do myself but reaching an expression which is different from expression of PSD which I remember for modulated signal

Actual answer $$S(f-f_c)+S(f+f_c)$$

My answer $$S(f-f_c)+S(f+f_c)+2X(f-f_c)X(f+f_c)$$

where $S(f)$ is PSD of $x(t)$ and $X(f)$ is its Fourier transform while it is modulated by $\cos$ signal of $f_c$ frequency.

What I have followed is to find fourier transform of modulated signal and then square it.Where i have done wrong.

If $x(t)$ is a finite-energy signal with Fourier transform $X(f)$, then $x(t)\cos(2\pi f_c t)$ is also a finite-energy signal with Fourier transform $\left.\left.\frac 12 \right[X(f-f_c) + X(f+f_c)\right]$. This is just the modulation theorem of Fourier transform theory.

The energy spectral density of $x(t)$ is $S_x(f) = |X(f)|^2$ while the energy spectral density of $y(t) = x(t)\cos(2\pi f_c t)$ is
\begin{align}S_y(f) &= \bigg|\left.\left.\frac 12 \right[X(f-f_c) + X(f+f_c)\right]\bigg|^2\\
&= \left.\left.\frac 14 \right[ |X(f-f_c)|^2 + |X(f+f_c)|^2 + 2\Re\big(X(f-f_c)X^*(f+f_c)\big)\right]\\
&= \left.\left.\frac 14 \right[S_x(f-f_c)+S_x(f+f_c)+2\Re\big(X(f-f_c)X^*(f+f_c)\big)\right]\tag{1}
\end{align}
which has a vague resemblance to what the OP wrote.  Now, if $x(t)$ is a low-pass signal with bandwidth $W$ Hz and $W < f_c$, then $X(f-f_c)$ occupies the frequency band from $f_c-W > 0$ to $f_c+W$ while $X(f+f_c)$ occupies the frequency band from $-f_c-W$ to $-f_c+W < 0$. We see that there is no value of $f$ for which $X(f-f_c)$ and $X(f+f_c)$ both are nonzero, and thus $X(f-f_c)X^*(f+f_c)$ equals $0$ for all $f$. We are left with
\begin{align}
S_y(f) &= \left.\left.\frac 14 \right[S_x(f-f_c)+S_x(f+f_c)\right]\\
&= \left.\left.\frac 14\right[|X(f-f_c)|^2+|X(f+f_c)|^2\right].\tag{2}\end{align}
But when $f_c$ is smaller than $W$, then, as the OP peevishly has informed us all that he has placed no restrictions on the value of $f_c$, he must learn to live with the more complicated $(1)$.


Does the use of a sparse basis in Compressed Sensing imply the need to have access to all the information beforehand? According to literature, the CS framework operates on the knowledge that most natural signals are sparse in some domain given by a sparsifying transform operation $\Phi$ (Fourier, Haar, WHT, etc.).

Following a single-pixel camera example, if we are trying to image a phenomenon described by a signal $x$, we use a sampling matrix $A$ to obtain a compressed measurement $y=A\Phi x$ of a representation of $x$ in the $\Phi$ domain, which is sparse.

My question is, assuming we have some kind of physical single-pixel device (using a DMD to project linear combinations of $x$ (given by the sampling matrix $A$) onto a detector to obtain a measurement vector $y$), where does the sparsifying transform operation $\Phi$ manifest in this device?
Do the DMD patterns correspond not to the sampling matrix $A$, but to the $A\Phi$ product?
Must we provide the DMD with the signal in the sparse domain, $\Phi x$? In that case, don't we require previous knowledge of $x$ in its entirety in order to obtain its transform in the $\Phi$ domain? It is my understanding this is not the case, but I fail to grasp the mapping between the mathematical framework and hardware implementation in this matter.

It seems to me you have a little misassumption here. During the sampling only $\Phi$ matrix is applied to the signal $x$, which is resulted in measurements vector $y$,
$y=Φx$ . Later, in reconstruction we use sparsity assumption to reconstruct the signal through solving problems of this form (there are other forms of solution):

$$\min|x|_1   s.t.   |\Phi x-Y|_2 < \mu$$

with $\mu$ being a small figure. This means find an $x$ with minimum norm 1 length, that satisfies the second term. Till now, we assumed $x$ to be sparse in sampling domain (in case of images it is space, in case of signals it is time and in case of video it time-space). But what if the signal is not sparse in the sampling domain and it is sparse in transform domain. Do we have to get the transform of the signal/image prior to measuring it through $Φ$ matrix? No! we don't. We can still sample the sparse signal in time or space domain but modify the reconstruction problem so as to find the sparse signal of transformed signal. So, if the input signal has an sparse representation in a transform domain with operator of matrix $\Psi$, we can write $x=\Psi^{-1}a$ or $a=\Psi x$ and substitute it in the reconstruction problem as follows:

$$\min |a|_1   s.t.   |\Phi\Psi^{-1}a-Y|_2 < \mu$$

So it is on reconstruction phase where we enter the sparsity and transform matrix to the problem and there is no need to have the transformed signal during the sampling.

Having said that, in response to your question:

  Does this mean the DMD pattern represents the AΦ product?

No, DMD pattern only represent $\Phi$. The $A\Phi$ product appears later on in reconstruction phase, and it is not $A\Phi$ if by $A$ you meant transform matrix . It rather would be $\Phi A$


Is there such a thing as band-limited non-linear distortion? So if you generate a square wave by just switching a signal between two values, at sample boundaries, it produces an infinite series of harmonics, which alias and produce tones below your fundamental, which is very audible.  The solution is Band-Limited Synthesis, either using additive synthesis or band-limited steps to produce waveforms that are the same as if you had band-limited the ideal mathematical square wave before sampling it:

http://flic.kr/p/83JMjT

But I just realized that if you apply large amplification to a digital sine wave and then clip it digitally, it will produce the same square wave shape, without the Gibbs phenomenon ripples.  So it's then also producing aliased distortion products, right?  So any non-linear distortion in the digital domain that produces harmonics outside of the Nyquist limits will produce aliased distortion products?  (Edit: I've done some tests and confirmed that this part is true.)

Is there such a thing as band-limited distortion, to simulate (in the digital domain) the effects of distorting (in the analog domain) before band-limiting and sampling?  If so, how do you do it?  If I search for "bandlimited distortion" I find some references to Chebyshev polynomials, but I don't know how to use them or if they only work for sine waves or what:

  This instrument does not attempt to generate band-limited distortion. Those interested in band-limited distortion should investigate the use of Chebyshev polynomials to generate the effect. Hyperbolic Tangent Distortion

  "Chebyshev polynomial" -- shaping functions with the important property that they are intrinsically band-limited i.e. they do not introduce spurious spectral harmonics due to overlapping etc. Wave Shaper

Applying a non-linear function will always introduce harmonics, and mixing non-linear functions with sampled versions of continuous signals does add the wrinkle you note above (where high-frequency harmonics are aliased to low frequencies.)

I can think of a few ways to proceed:

You can use an oversampling factor high enough to capture the extra harmonics (up to some arbitrary precision, e.g. your noise floor),
You can use a "softer" clipping function (see, for example, here) that has harmonics that die out sooner than the hard clipper. This is easier to model, but introduces its own distortion at low frequencies.
Building on the approach you suggested above, interpolate your sampled signal (e.g. using a Lagrange or Chebyshev interpolator) to construct a continuous-time model. Then, apply the hard clipper and low-pass in a simulated continuous-time domain. Sample the result.

You can combine (1) and (2). The third approach is complex, but gives you the best control over how much distortion to admit, and will probably scale better to very high fidelity requirements.

For non-linear functions that admit a series expansion (e.g. Taylor/Maclaurin), you can get a decent intuition for how fast the harmonics decay. The Maclaurin expansion of a function $f(x)$ is:

$$f(x)=\sum_{n=0}^\infty \left[ \frac{f^{(n)}(0)}{n!} x^n\right]$$

In your case, $f(x)$ is the clipping function. (You can't do this with a hard clipper, at least not naively!) If you consider the substitution $x=g(t)$, where $g(t)$ is your input signal, the $x^n$ becomes $g(t)^n$, which you can consider to be the convolution of your input signal with itself $n$ times. Thus, for low-pass signals, the $n$th term of the infinite summation has a bandwidth $n$ times that of your signal. To complete the picture, you need to figure out the amplitude associated with each term and decide how many terms in the summation are relevant.

(With a little thought, you might also be able to use this form directly to approximate the filtered non-linearity. That would require a good series representation for the clipper.)


Nyquist with only poles at the origin I need to solve this three exercises:
$$G_1(s) = \frac{1}{s}$$

$$G_2(s) = \frac{1}{s^2}$$

$$G_3(s) = \frac{1}{s^3}$$

This is what I did:

For the single pole ($1/s$), magnitude of A is $\infty$ and phase is 0

For the single pole ($1/s$), magnitude of B is $\infty$ and phase is -45

For the single pole ($1/s$), magnitude of C is $\infty$ and phase is -90

For the 2 poles one ($1/s^2$), magnitude of A is $\infty$ and phase is 0

For the 2 poles one ($1/s^2$), magnitude of B is $\infty$ and phase is -90 ($=2*(-45)$)

For the 2 poles one ($1/s^2$), magnitude of C is $\infty$ and phase is -180 ($=2*(-90)$)

For the 3 poles one ($1/s^3$), magnitude of A is $\infty$ and phase is 0

For the 3 poles one ($1/s^3$), magnitude of B is $\infty$ and phase is -135 ($=3*(-45)$)

For the 3 poles one ($1/s^3$), magnitude of C is $\infty$ and phase is -270 ($=3*(-90)$)

MATLAB can't do these Nyquist analysis in a correct way. That's why I'm asking here.

Thanks

Your answer is correct. Each pole gives 90 degrees phase lag. Furthermore, here can you download a file lnyquist1 which is able to give a more accurate result when you have a pole at zero.

num = 1;
den = [1 0 0 0];
lnyquist1(num,den);

In addition, WolframAlpha is also able to give an accurate result, see here.


Detect valleys of a signal I have raw values that are coming directly from a sensor (on the fly). These values are coming from a wearable sensor installed on a cows neck. When the cow ruminates (meaning that it chews food that brought back from the 2nd stomach of it), generates some kind of acceleration values that have some kind of specific amplitude and frequency (the chew duration is mostly between 800 msec to 1200 msec). A rumination period is between 30 to 120 seconds every time and the cow ruminates around 10 hours per day.

I do a quick calculation on these values and I end up with a graph like the bellow image.

I would like to be able to detect every time when the signal falls (like the highlighted peaks) as it is a representation of a rumination period. That fall in the graph represents the bolus event of the cow (meaning that the cow just shallow-ed the chewed food and brought back a new one [from the stomach] to chew again).

I tried to do it with a threshold approach (if the signal falls bellow a set threshold), but it fails to operate always correctly. This is because:

1) The signals intensities change slightly during the day and between different sensors eg. On the image, the negative peaks do not fall bellow the black horizontal line which is my threshold. As such I have false detection (not detecting them at all). In case I increase the threshold, I do wrong detections since a lot of other points in the signal might pass it.

2) When the signal is slightly noisy I do wrong detections or no detections at all. As an example, the 9th yellow point of the signal should be somehow detected as valley but it is well above the threshold. Of course it could be mistakenly detected as it has the same amplitude as some nearby points.

Does anybody know how I can detect such points (dips) in time?

Note:
These data are coming real-time and they are not stored somewhere. As such I must somehow be able to detect dynamically the fall of the signal without having any knowledge of the future values of it (well I can have some small history of course - lets say last 200 values or something).

Extra info:
These data values are actually coming from a 3-axis accelerator sensor. I calculate the magnitudes of them (sqrt(x^2 + y^2 + z^2)) which they deviate around 1g and then I get the Moving Standard Deviation of the last second (62.5 values).
That is what is displayed on the graph above. The yellow points are at around 0.01g deviation.

Further:
After being asked by some users, I also attach an image of the original signal (only a tiny part of it). With yellow are marked the areas of interest (on the magnitude). The rest 3 graphs are the X,Y,Z (legend exists on the top-right of the graph).

Also, I attach a hyperlink where you can download sample signals in CSV format. The format is as following:
1st, 2nd, 3rd columns are the X, Y and the Z axis accelerometer values.
The 4th column is the magnitude of them (sqrt(x^2 + y^2 + z^2)) which you can also manualy calculate if you want so.

https://1drv.ms/f/s!Au5-DmkSvYQzgZlLo2KJ_UBAVeYDWw

There are files with 4 quite good signals. The Good3s.csv contains some small junk at the start of it (some 10000 first values can be discarded). The Good4e.csv has the same at the end. There are also 2 more signal files that are quite noisy and you can parse them if you like hardcore situations.

First, you'll probably have better luck posting this on dsp.stackexchange.  That's a more specialized group that does stuff like this all the time.

In terms of your problem, here's a couple of options.

One is a machine learning approach.  e.g. create a training set by taking a bunch of data and hand marking the points that are good versus bad (like you've done above).  Then take an Artificial Neural Network, feed it the training data and teach it to recognize the pattern you want.

Another approach might be to take a Statistical process control approach.  Lots of tools there for detecting when something gets out of line.

Another way might a floating threshhold.  Compute the mean of the last X samples, and if the signal drops more than Y below that mean, call it bad.  could also compute the standard deviation of the last X samples, and if the signal drops more than Z std devs then call it back.

or, potentially just high-pass filter the signal to remove the slowly drifting mean and then you can use a fixed threshhold.


transmitting the real/imaginary parts of pass-band signal As I know, the transmitted pass-band signal whose carrier frequency is $f_c$ can be written as follow:
$s_1(t) = \text{Re}(a(t)e^{j2πf_ct})$ , where $a(t)$ is the complex modulated base-band signal. It's clear to notice that at the receiving end,   we can recover the complex modulated symbol $a(t)$ based on the real transmitted part $s_1(t)$.
But the question which I couldn't get a sure reply for it, what's about if we transmitted the imaginary part of our signal? .. it means that  $s_2(t) = \text{Imag}(a(t)e^{j2πf_ct})$. Will we be able to recover the complex modulated symbol based on $s_2(t)$ too?  How can we prove that mathematically?
thank you in advance

Will we be able to recover the complex modulated symbol based on $s_2(t)$ too?

Yes, no information is lost between the two.

How can we prove that mathematically?

One way to think about it is to write out each individually and see what you have. Let $a(t)=a_I(t)+ja_Q(t)$, expanding $s_1(t)$ and $s_2(t)$ we have:
\begin{align}
s_1(t) &= a_I(t)\text{cos}(2\pi f_ct)-a_Q(t)\text{sin}(2\pi f_ct) \\
s_2(t) &= a_I(t)\text{sin}(2\pi f_ct)+a_Q(t)\text{cos}(2\pi f_ct)
\end{align}
The in-phase and quadrature parts are present in each so we just need to bring them back to baseband to get $a(t)$ back. Further, whatever workings you have to do this for $s_1(t)$ you can reuse for the processing of $s_2(t)$ by first applying a simple $+\frac{\pi}{2}$ to it ($\text{cos}(\theta + \frac{\pi}{2})=-\text{sin}(\theta)$).


How to build a realtime simple audio spectogram? I'm trying to build a spectogram view of history time of 4 seconds and I'm stuck at the part where you load/draw FFT vectors into the 2D bitmap. I just don't know how to do that, but i got the FFT in realtime.

I'm doing it in swift(iOS)

Thank you!

The Discrete Fourier Transform of a signal is a series of $\mathbb{C}$omplex numbers. When performed in short time segments and depicted in a spectrogram, what is really depicted is a spectrum.

That is, the strength of a particular frequency component, which is a $\mathbb{R}$eal number.

There are various different expressions for the strength of a component of the spectrum such as Amplitude, Power, Power Spectral Density and so on.

The most straightforward one is the Amplitude spectrum. To calculate the Amplitude of a given spectrum component, you need to get the absolute value of the complex number that represents a frequency component in your spectrum.

Programming languages such as C, Python, Julia, MATLAB / Octave and others, have native complex data types. In that case, the output from an fft() kind of function is an array of type complex and each element of that array contains a function (e.g. abs()) to obtain the absolute value of that complex number. If your programming language does not have a complex data type, then the output of your fft() is probably two arrays of real numbers. One for the real part of the complex number and one for the imaginary part. In that case, to obtain the spectrum, all that you have to do is spc[m] = sqrt(real_part[m]^2 + imaginary_part[m]^2) where ^ denotes exponentiation and m is the $m^{th}$ element of your DFT spectrum.

Once you do that, you will have an array of real numbers. These are going to vary wildly, because of the way the DFT sums work and for this reason you usually apply a logarithmic transformation. This transformation keeps large numbers large and makes small numbers large. You would typically apply this with something like spc[m] = log(1 + spc[m]).

At this point, you can depict spc[m] in grayscale by normalising it to 255 (usually) tones of gray with something like rendered_spc[m] = int(round(spc[m]*255.0)). int denotes an integer data type, you are likely to use something like a BYTE or whatever is more convenient for your bitmaps. The other thing you can do is to map the value of spc[m] within a palette of N colours and "paint" the pixels with that colour.

Finally, it would be probably useful to emphasise here that the spectrogram does not depict one "real time FFT" but many short time fourier transforms. So, if you want to depict 4 seconds worth of a histogram, you would have to accumulate a number of FFTs, get their amplitude spectra, transform them logarithmically, map them to colours and then depict them. Think of each FFT amplitude spectrum as one column of your bitmap.

This is a brief outline on rendering spectrograms practically. Depending on the desired accuracy, you might have to add some more complex features, such as "overlapping" DFT segments, interpolation (for smoother bitmaps with less data) and others.

Hope this helps.


How to understand the sum of all "fourier frequencies"? I got the data of an acceleration sensor to analyze. It consists of special terms of 30 Hz, 60 Hz and 120 Hz. In the following you can see in the first plot the 60 Hz data and in the second one the spectrum:

The plotted values are the results of a FFT. My question now is: How to understand the spectrum?
The 60 Hz plot is clear. A FFT was operated and all the values were extracted which correspond to the 60 Hz signal. This is (or should be) the upper plot.
The spectrum is allegedly the value across the whole frequency range - this I do not quite understand. How can I make use of this or how can it physically be interpreted? Is it finally kind of a total amplitude? Is it the origin signal (I don't think so..)? Is it even possible without an inverse transformation?
Or is it just what it is: The sum of all (fourier) frequencies without any higher meaning? I'm not sure how to name the signals/peaks in the frequency domain.

Let me put a practical answer with the following Matlab / Octave Code :

L = 2*1000;     % signal sample count
n = 0:L-1;      % discrete-time index

Fs = 44100;     % sampling frequency
am = [1, 0.5, 2, 0.5, 0.3, 0.6, 0.1, 0.2];                % magnitudes of 8 components
fm = [882, 2646, 4410, 6615, 8820, 10000, 13230, 15876];  % frequencies

% Time domain signal :
% x = a1*sin(2*pi*f1*n) + a2*sin(2*pi*f2*n) + ... + am*sin(2*pi*fm*n)
x = sin((2*pi/Fs)* n'*fm)*(am');

% Frequency domain signal :
X = abs( fft( x, L ) );

% Display :
figure,subplot( 2,1,1 );
plot( n/Fs , x );
title(' signal as a sum of sinusoidal frequencies ' );
xlabel(' time [s]');
subplot( 2,1,2 );
plot(linspace( 0 , Fs/2 , L/2 ) , (2/L)*X(1:L/2) );
title(' Fourier (frequency) spectrum magnitude of the        signal' );
xlabel(' Frequencies from 0 to Fs/2 [Hz]' );
ylabel(' FFT Magnitude ');

The output will be :

The output speaks for itself. The upper plot is the time-domain signal as a superposition of (M=8) different sinusoidal signals each with a particular frequency and amplitude. The lower plot shows its frequency spectrum magnitude which displays a peak per frequency component in $x[n]$. The peak's amplitude is proportional to the time-domain component amplitude.


Why does $\mathrm{sinc}[k]= \delta[k]$ hold? On page 5 of this document I read that

$$\mathrm{sinc}[k] = \delta[k]$$

Could anyone please explain why this is true?

EDIT: After editing according to the comments, the real question is if and why the following equality holds:

$$\text{sinc}[k]\doteq\frac{\sin(k\pi)}{k\pi}=\delta[k],\quad k\in\mathbb{Z}\tag{a}$$

The answer is, yes, of course it holds. Since $\sin(k\pi)=0$ for any integer $k$, the value of the left-hand side of the above equation must be zero for any integer $k\neq 0$. The case $k=0$ is different because here also the denominator is zero. It can be shown by a limit argument that for $k=0$, the function $\sin(k\pi)/(k\pi)$ equals $1$. Another way to see it is to use the equation that you quoted in your comment:

$$\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{j\omega k}d\omega=\ldots=\text{sinc}[k]$$

If you evaluate this equation separately for $k=0$ you get

$$\frac{1}{2\pi}\int_{-\pi}^{\pi}d\omega=1$$

So we get a sequence which is $1$ for $k=0$ and zero for all other (integer) values of $k$. Consequently, this sequence can be rewritten using the discrete-time unit impulse $\delta[k]$, as shown in Equation (a).

(Below is the original answer explaining that $\delta[k]$ can't be the DTFT of any sequence.)

As mentioned in the comments, this can't be true. Be careful to distinguish between the Discrete-Time Fourier Transform (DTFT) and the Discrete Fourier Transform (DFT). The first is defined by

$$\begin{align}\text{DTFT:}\quad X(e^{j\omega})&=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega}\\x[n]&=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(e^{j\omega})e^{j\omega n}d\omega\end{align}\tag{1}$$

The DTFT $X(e^{j\omega})$ of a sequence $x[n]$ is a $2\pi$-periodic function that is continuous in $\omega$.

The (length $N$) DFT is defined by

$$\begin{align}\text{DFT:}\quad X[k]&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\\x[n]&=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\end{align}\tag{2}$$

Here, $X[k]$ is also a discrete sequence.

The discrete unit impulse $\delta[k]$ can be the DFT of some time domain sequence, but it cannot be the DTFT of any sequence (because it is itself a sequence, not a continuous function). To see which sequence $x[n]$ corresponds to the DFT $X[k]=\delta[k]$, simply use the formula given in $(2)$:

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}\delta[k]e^{j2\pi nk/N}=\frac{1}{N}\tag{3}$$

From $(3)$, it is the constant sequence $x[n]=1/N$ that corresponds to the unit impulse $X[k]=\delta[k]$.


system function $H(\omega)$ relationship to odd and even components of h[n] What qualities of $h[n]$ are necessary for:

$$
H(e^{j\omega}) = DTFT\{h_{even}[n]\} + j\ DTFT\{h_{odd}[n]\}
$$

Do all real / causal h[n] have the property that:

$$
H(e^{j\omega}) = DTFT\{h_{even}[n]\} + j\ DTFT\{h_{odd}[n]\}
$$

where:

$$
h_{even}[n] = \frac{1}{2}(h[n] + h[-n])
$$

$$
h_{odd}[n] = \frac{1}{2}(h[n] - h[-n])
$$

The DTFT relationships

$$x_{even}[n]=\frac12\left(x[n]+x^*[-n]\right)\Longleftrightarrow\textrm{Re}\left\{X(e^{j\omega})\right\}$$

and

$$x_{odd}[n]=\frac12\left(x[n]-x^*[-n]\right)\Longleftrightarrow j\,\textrm{Im}\left\{X(e^{j\omega})\right\}$$

hold for any sequence $x[n]$ for which the DTFT exists. There is no assumption about $x[n]$ being real-valued or causal (note the complex conjugation $^*$ in the definition of even and odd signals). If $x[n]$ is real-valued you can leave out the conjugation.

Note that the DTFT of the odd part $x_{odd}[n]$ equals $j$ times the imaginary part of the DTFT $X(e^{j\omega})$, so you have

$$X(e^{j\omega})=\textrm{DTFT}\{x_{even}[n]\}+\textrm{DTFT}\{x_{odd}[n]\}$$

(without a $j$ on the right-hand side).


Audio fingerprinting - 2D (time - frequency) filter / Energy differences I am trying to implement the algorithm described in this paper. However, I am having trouble understanding the last part where I need to find the energy differences of bands.

In another paper, this process is described as shown below. Also, I am not sure I understand that either. What exactly I should do in this part of the algorithm ? Any hints or tips would be really appreciated.

If you split the signal in 33 Frequency bands for every frame you will have a vector band containing all subbands, all right, and now to compute the Energy Difference from each Band you just need to do:

$$E(x)=\sum_{k=1}^{32}band(k) - \sum_{k=1}^{32}band(k+1)$$

For the first frame compute the equation above, and for all others frames block, You can do something like:

$$E2(x)=\sum_{k=1}^{32}band(k) - \sum_{k=1}^{32}band(k+1)$$

$$F(n,m)= 1 if(E2 -E)>0$$
$$F(n,m)= 0 if(E2 -E)<=0$$

$$E(x) = E2(x)$$

This will generate a 32 bits array (0's and 1's) in $F$, Your end result should be a matrix of dimension 32x256, and it is your fingerprint !

EDIT

Answering your other question, if you are following the article exactly as described, the frame rate used in the paper is 5000hz, then each sub-fingerprints 32x1 is equivalent at 0,0116ms of audio its give you an analyze of 58 Frames = 0,0116 * 5000, it does not contain enough information to make a fingerprint, you'll need something comes close to 3 seconds and 256 comes very close look, 0,0116 * 256 = 2,9696, Ok you need do one fingerprint for ever 256 sub-fingerprints and it will create a lot of 32x256 matrix for your audio :-)


Algorithm to compute the Variance of the signal? I have 6.4 sec recording of a discrete signal on 40Hz frequency (256 values).

How can I compute the variance of the signal?
Please, share an algorithm in pseudo-code or in an actual programming language.

You want the sample variance:

$$
s^2 = \frac{1}{N-1} \sum_{n=0}^{N-1} (x_n - \mu_x)^2
$$

where the mean, $\mu_x$, is given by

$$
\mu_x = \frac{1}{N} \sum_{n=0}^{N-1} x_n.
$$

The $\frac{1}{N-1}$ factor should be used to account for the bias that is otherwise there if $\frac{1}{N}$ is used. One way to think about this is the sample variance is computed from $N-1$ differences between samples, rather than $N$ samples.


Numerically generating noise time series from spectral density I'm trying to figure out a way to construct (numerically) time series of random variables given some spectral density $S_{xx}(f)$. I'm interested in various types of spectral densities (flat, lorentzian, 1/f) so if there is a general method that starts from an expression of $S_{xx}(f)$ and generates a suitable time series that would be ideal of course, but I am not sure if there is.

Honestly I am not so sure where to start. My guess is that you compose your random signal as a sum of different components with amplitudes related to your spectral density at a certain frequency. I understand that if you want your noise to have a complete frequency range you need infinitely many terms, so I suppose you cut off your noise at some frequency $f_{max}$, giving you some maximum amount of terms $K$. You'll also need some step size $\Delta f$ for the different frequency components, but after that I get kind of stuck. How do you incorporate your spectral density into this?

I understand that the question is perhaps a bit too general, so I'm also very much interested in advice on some sources on the subject. The ones I have been able to find are very technical, while for now I'm simply trying to create some time series in Mathematica given some spectral density as input, after which I want to see if I can reconstruct the spectral density from the output. That would get me started quite well, and I can think about what features are essential from there on.

As mentioned by @MBaz in comments, you could use the fact that filtering white noise with unitary power spectral density (for example from independent Gaussian samples zero mean and unitary standard deviation) with a filter having a frequency response $H(f)$ gives you colored noise with power spectral density $S_{xx}(f) = |H(f)|^2$.

Then comes the question of designing such a filter. A flexible approach would be to use a Frequency sampling-based FIR filter design. Note that this is readily available as FrequencySamplingFilterKernel in Mathematica (or similarly fir2 in matlab for future readers), using an evenly spaced sampling of $|H(f)| = \sqrt{S_{xx}(f)}$ up to the Nyquist frequency (i.e. $f_{max}$ which would correspond to half your time-domain sampling rate). Once you have the FIR filter coefficients, you can use those to filter your white noise with ListConvolve.

Note that the steeper the slopes in your power spectral density, the more frequency samples (resulting in a higher number of FIR filter coefficients) you'll need to get an accurate representation of your power spectral density.


Camera image pixel structure origin I am using machine vision for industry task with opencv. I have images from a camera monitoring the process. The detailed look at the small part of the image follows:

Enlarged:

In the detail, it is possible to see pixel chessboard-like pixel structure. I can get rid of it in many ways (any kind of blur or averaging will do the trick). However, I would like to know where this structure origins. My ideas:

camera use colors by default - can the grayscale setting use the RGB sensors in the way that cause this?
could this be produced by compression or storing of the images? I have no information about how it is done - but the images are in PNG format when I get them.

As I said, I can homogenize the image, but I wonder, maybe there is something fundamentally wrong what we do and that should be fixed.

It looks like that is a color camera with the Bayer pattern, and the color camera data is directly used as grayscale data, as if it was coming from a grayscale camera. If it is a color camera, the color camera data should be de-Bayerized first (and smoothed etc processing) to get color RGB data pixels as usual, and then the RGB data can be converted to grayscale if grayscale is needed.


What would be the effect of repeating samples in time domain on frequency domain? Consider the operation of up-sampling a signal n times by repeating each sample n times in the time domain. What would the frequency domain response of the resultant signal be like given the frequency response of the original signal?

It's quite benefical to view the mathematical expression of the process as well.

In discrete-time interpolation, a mere $N$ times repetition of a single sample (pixel) into the new enlarged block is an interpolation with zero-order hold kernel. This is mostly applied for simple fast but crude image zooming.

In 1D, the kernel is $$h_0[n] = \begin{cases}{ 1 ~~~,~~~ 0 \leq n \leq N-1 \\ 0 ~~~,~~~ \text{ otherwise } } \end{cases} \tag{1}$$

The interpolated signal can be expressed as :

$$y[n] = \sum_{k} x[k] h_0[n-kN] \tag{2}$$ this can be expressed as a convolution like

$$y[n] = \left(\sum_{k} x[k] \delta[n- Nk] \right) \star h_0[n] \tag{3}$$

where the signal in parenthesis is the expanded version of $x[n]$:

$$x[n] \rightarrow \boxed{ \uparrow N} \rightarrow x_e[n] = \sum_{k} x[k] \delta[n- Nk] \tag{4}$$

Applying DT Fourier transform on Eq(3) yields:

$$ Y(e^{j\omega}) = X_e(e^{j\omega}) \cdot H_0(e^{j\omega})  \tag{5} $$

where $X_e(e^{j\omega}) = X(e^{j\omega N}) $ due to the expansion operation in Eq(4) and $H_0(e^{j\omega})$ is the DTFT of the rectangular pulse $h_0[n]$ which is:

$$H_0(w) = e^{-j\omega (N-1)/2} \frac{ \sin( \omega N /2 ) } { \sin( \omega / 2 ) } $$

then the result becomes:

$$ Y(e^{j\omega}) = X(e^{j\omega N}) \cdot e^{-j\omega (N-1)/2} \frac{ \sin( \omega N /2 ) } {\sin( \omega / 2 ) }  \tag{6} $$

Whose magnitude is :
$$ \boxed{ |Y(e^{j\omega})| = |X(e^{j\omega N})| \cdot \left| \frac{ \sin( \omega N /2 ) } {\sin( \omega / 2 ) } \right| }  \tag{7} $$

Note that the $X(e^{j \omega N})$ is composed of repeated images of $X(e^{j \omega})$. And the effect of the interpolation kernel $h_0[n]$ is to remove the images as much as possible with that crude frequency response $H_0(e^{j \omega}) $.

The following Octave/Matlab code demonstrates the operation:

clear all; close all; clc;

N = 8;    % interpolation
L = 32;   % original signal length
x = 5*( sin(0.3*pi*([0:L-1]-L/2))./(pi*([0:L-1]-L/2))).^2;
x(isnan(x))=5*0.3*0.3;

h0 = ones(1,N);

xe = zeros(N*L,1);  % expanded signal
xe(1:N:N*L) =  x;

yi = conv(xe, h0);  % the interpolated signal: y[n] = xe[n] * h0[n]

figure,subplot(2,1,1);
stem(x);
title('original signal');
subplot(2,1,2);
stem(yi);
title('zero-order hold interpolated output');

figure,subplot(2,1,1);
plot(linspace(0,2,1024),abs(fft(x,1024)),'b');
title('Spectrum magnitude of x[n]')

subplot(2,1,2);
plot(linspace(0,2,1024),abs(fft(xe,1024)),'b');
hold on
plot(linspace(0,2,1024),abs(fft(h0,1024)),'r');
plot(linspace(0,2,1024),abs(fft(yi,1024)),'g','linewidth',2);
legend('|X_e(w)|','|H_0(w)|','|Y_i(w)|');

with the following outputs:

and the spectrum:

as you can see, the nulls of the periodic sinc (the Dirichlet kernel) are at the centers of the repeated spectra (the images) of $X(w)$. However this is not enough to suppress them all unlike an ideal lowpass filter (or one of its adequate approximations) which would completely eliminate the image spectra and yield only the one inside the main lobe; hence yielding a perfect interpolation compared to this crude zero order hold kernel giving.


Does there exist a simple regression method to fit a single sinusoid period to a data set? I'm currently working on a project where I want to fit a single sinusoid period to a data set. Essentially I have very good control/knowledge of the signal's dominant frequency, so I'm only sampling for a single period. What this leaves me with is the problem of determining the amplitude and phase offset in the presence of the sensor's noise.

I'm curious if anyone is aware of a easy to implement method for estimating these two parameters that is optimal in some sense. Most of my research results have been complex algorithms implemented in MATLAB which is not suitable for my application (I'm trying to measure this in real-time and don't have a license).

A simple solution is to calculate the first coefficient of a DFT of appropriate length, using the summation formula instead of FFT. To get amplitude and phase, transform the result to polar coordinates.


How to fix audio clipping in Adobe Audition? I am editing audio from video that was captured from the sound board. The audio was recorded too loud since someone turned up the board and and the bass and vocals sound scratchy. I have brought the levels down a bit to normalize them but still have the scratchy sound with the vocals. Is there an effect in Adobe Audition to help this?

In Adobe Audition 3 there is Effects... Restoration... Clip Restoration (process).


Video stitching for static scene (infinite loop video) I have a video that is a panoramic video of a static scene. Just to be clear, I am not talking about a spherical 360° degree video but a video that starts in some direction and then does a full 360 horizontal tour either clock or counterclock wise (direction is known). The video then cuts approximately off when a complete tour is completed (~10° error). The video is most likely between 10 to 20 seconds long. I would like to modify the beginning and end of the video of in such a way that if I play the video in a loop the user has the impression he sees a never ending video (end and beginning should have a seamless transition). The main problem is that beginning and end might be slightly vertically shifted and tilted with respect to one another. Otherwise I would only need to find a frame that matches.

My question is therefore I do best solve this problem. Should I look into image stitching or are there more specific algorithms for this type of video stitching. When I search for video stitching I mostly come across 360° video stitching. Pointers to the appropriate literature are also very welcome.

EDIT: You might be wondering why I don't just create a panorama picture as the scene is static. I want to take advantage of the camera that adapts to the lightning changes.

Lets assume your camera axis has polar and azimuthal angles and you want to align your frames only regarding the azimuthal angle and set the polar angle free.

First you have to find the geometric transform between sequential frames of your video, then by multiplying these sequential geometric transforms between your reference frame and the current frame you could obtain the geometric transform between your current frame and the reference frame. To obtain geometric transform between sequential frames you could use automatic point matching algorithms like SIFT or SURF (considering all your frames don't have any intersections but sequential frames have intersection I suggested to use it sequentially).

After finding the geometric transform you have to modify the geometric transform to set the polar direction of camera free. To do this You have to decompose the geometric transform to find external camera parameters which contain the position and the pose of camera. Pose of camera is a Rotation matrix which rotate the axis of camera between your reference and current frame. After obtaining rotation matrix you have to decompose the rotation matrix to obtain the rotation angles which are yaw, pitch and roll and set Yaw angle to zero (I'm not sure which one is the polar angle of your camera but any of these angles which is polar angle of camera must be zero to set the polar direction free), then you have to obtain the modified geometric transform from the modified rotation matrix.

Now you could align all your frames using the inverse geometric transform.

All parts of suggested method have intrinsic MATLAB functions, also if your camera has nonlinear distortion like fish-eye you have to take care of that before point matching part.


Finding the Correlation Matrix of Wiener Filter using Matlab I have the problem below and need to solve it using Matlab. I not a big expert in Matlab and just don't know how to implement it.
Basically, it says that
\begin{equation}
R\:=\:\left[X\left(n\right)X^T\left(n\right)\right]
\end{equation}
but just don't know how to do that in Matlab. I searched a lot on the internet and didn't find anything relevant to what am looking for.
Can somebody please help me implement this. I would really appreciate it.

Here is a simple example of Wiener filter in MATLAB. First, the input data, denoted as $\mathbf{x}(n)$. Note that we are considering only at time instant $n$ and an extension for $N$ time samples can be made to incorporate more time domain samples by concatenating these vectors together to get a matrix $\mathbf{X}=[ \mathbf{x}(n) \ \ \mathbf{x}(n+1) \ \ \cdots \ \mathbf{x}(n+N-1)]$ (that is a vector for each time instant). But for the single time instant case, lets say that $\mathbf{x}(n)$ is IID Gaussian with some input mean and varaince, then this input vector can be generated in MATLAB as:
x = sqrt(inputVariance)*randn(10, 1) + inputMean;

The 10 comes from your problem description. Now that you have some example input data, you can go ahead and compute the correlation matrix as you defined $\mathbf{R} = \mathbb{E}\big[ \mathbf{x}(n) \mathbf{x}(n)^H \big]$:
R = x*x'/size(x, 2); % size(x, 2) = # of time domain samples

Since our example input data is easy to understand (Gaussian), we can know what to expect to find in the correlation matrix. We generated the input data to be uncorrelated, so the off diagonal elements should be small and the diagonal elements should be close to the variance. One useful way to look at these matrices in MATLAB is as an image using:
figure
imagesc(R)

After doing this, you might be confused though. Because you expect to see something like this:

But instead you might see something closer to the next figure. This is because you are only taking a single time domain sample! Not enough samples for the statistical properties to emerge from the correlation matrix. In fact, it is an estimate of the true correlation matrix and you can't expect an estimator to be that good with only one sample (the first image was generated by following the steps but using x = randn(10, 100) to generate 100 time domain samples).


Signal vs. noise at a particular frequency I collect samples (time series) at a frequency.  The sample rate is arbitrary (I can change it to whatever I want.)  I want to determine if what's at this frequency is a signal, or (normally distributed) noise.

I've tried using a normality test (from scipy); that didn't seem to work - there were no big differences in the results when what was at the frequency was a man-made signal (center frequency of an FM broadcast station, for example) vs. noise.  I tried using auto-correlation using Pandas - Pandas always thought it was noise.

Are my techniques off?  Perhaps the techniques are good, but I need to sample different (I've been sampling about every two seconds, for, say, 100 samples)?  Something else?

The usual scenario is where you want to find if the received signal is noise: $$r(t)=n(t),$$ or if it is a signal plus noise: $$r(t)=s(t)+n(t).$$ It is generally assumed that you know $s(t)$. For example, you could have a radar application where you transmit a pulse and need to find if it was reflected back; or a digital communications receiver, where you want to see if the signal received was $p(t)$ (say, indicating a bit 0) or $q(t)$ (a bit 1).

The way to do this is by calculating the correlation between $r(t)$ and $s(t)$. The idea is that if the signal is present, you have
\begin{align}
\text{corr}(r(t),s(t)) &= \text{corr}(s(t)+n(t),s(t))\\
&=\text{corr}(s(t),s(t))+\text{corr}(s(t),n(t))
\end{align}
which is large. If the signal is not present, then the correlation is small.

This scheme also works in the digital domain, assuming that you sample fast enough. In general you want to calculate the correlation of a large number of samples, so that (by the law of large numbers) the probability of being unlucky and getting "bad" noise samples is small.

In your question, you don't say if you know what the signal is. If you don't, but you know some of its properties, you may be able to look for it. For instance, if the noise is zero mean and the signal is not; or if you know that the signal's energy is concentrated in a certain band (the noise is flat). However, a specific answer is impossible without more information.

If you don't know anything about the signal, detecting it may well be impossible.


DFT and perfect reconstruction of a square wave on a digital computer I know that in theory, when reconstructing a square wave from its Fourier coefficients, unless we have an infinite amount of them, the resulting reconstruction will have Gibbs ringing artifacts due to lack of enough harmonics.
On a computer, we can take the Fourier transform X = fft(x) of a square wave x, and reconstruct it without artifact with x_rec = ifft(X), maybe with some rounding error of the order of 1e-17 or something but no visible ringing.
I don't have a satisfying answer for that? I guess there has to be something to do with the fact the "the square wave" x is a digitized version of a continuous wave, and my Fourier basis vectors (complex exponentials are also discretized of course since we are in a computer...) but still... how would you justify the absence of Gibbs ringing artefacts from Fourier reconstruction of the Fourier transform of a digital square wave ?
%%%%%%%%%%%%%%%%%%%% Tought experiment proposed by  Dan Szabo
fs=10;%sampling frequency
t=0:(1/fs):1-(1/fs);

s = [1 1 1 1 1 0 0 0 0 0];
sTr = imtranslate(s,[0.5 0])

sTr =

    0.5000    1.0000    1.0000    1.0000    1.0000    0.5000         0         0         0         0

First of all, it's a misunderstanding that the Gibbs phenomenon disappears if you use infinitely many Fourier series coefficients to reconstruct a discontinuous periodic function, such as a square wave. It doesn't. The reason is that generally the Fourier series doesn't converge point-wise, but it converges in the mean, i.e.,
$$\lim_{N\to\infty}\int_{0}^{T}\left|x(t)-\sum_{k=-N}^{N}c_ke^{j2\pi kt/T}\right|^2= 0\tag{1}$$
if $x(t)$ is a $T$-periodic function, and $c_k$ are its Fourier coefficients.
Taking the discrete Fourier transform (DFT) of a finite sequence of numbers just corresponds to a matrix multiplication:
$$\mathbf{y}=\mathbf{Ax}\tag{2}$$
and as long as the matrix $\mathbf{A}$ is invertible you can compute $\mathbf{x}$ from $\mathbf{y}$:
$$\mathbf{x}=\mathbf{A}^{-1}\mathbf{y}\tag{3}$$
This has nothing to do with the Fourier series of a continuous periodic function, and it has nothing to do with the Gibbs phenomenon.


Prove that, in DFT, Upsampling in time domain is equal to replication in frequency domain Given:
$x[n]$ is an $N$-point sequence whose DFT is $X[k]$
$$x[n]\xrightarrow{\mathcal{DFT}} X[k]$$
then,
Prove that: DFT of the same sequence after insertion of $(M-1)$ zeroes between successive samples is $M$ times repetition of DFT of original Sequence .
$$\begin{aligned}x\left[\frac{n}{M}\right]\xrightarrow{\mathcal{DFT}}\{&X[0],0,0,\dots,0,\\&X[1],0,0,\dots,0,\\&X[2],0,0,\dots,0,\\&\vdots\\&X[N-2],0,0,\dots,0,\\&X[N-1],0,0,\dots,0 \}\end{aligned}$$
My attempt:
$$ \mathcal{DFT} \left\{x\left[\frac{n}{M}\right] \right\}=\displaystyle\sum_{n=0}^{n=MN-1}x\left[\dfrac{n}{M}\right]W_{N}^{nk}$$
replacing $n\to{ nM}$
$$\text{LHS}=\displaystyle\sum_{n=0}^{n=N-\frac{1}{M}}x[n]W_{N}^{Mnk}\neq \big(X[k]\big)^M$$
I'm doing some mistake that is why i'm not getting $\text{RHS}=\big(X[k]\big)^M$.
Can anyone help me in proving this property?
Note:
$x[n]$  have period $ N$ so, $x[\frac{n}{M}]$ will be having period $MN$.
$W_{N}$ is twiddle factor $=N$ th root of unity$=\exp\left(-j\tfrac{2\pi}{N}\right) $

We start with
$$X[k] = \sum_{n=0}^{N-1} x[n] \cdot W_{N}^{nk}$$
Let's define a new sequence of $y[r]$ of length $R=M \cdot N$, where $R, M, N \in \mathbb{Z} > 0$

$$y[r] = \begin{cases}
y[n \cdot M] = x[n] \qquad & r = n \cdot M \\
\\ 0 \qquad & \text{otherwise} \\
\end{cases}$$

Transform it
$$Y[k] = \sum_{r=0}^{R-1} y[r] \cdot W_{R}^{rk}$$

Now most of the samples of $y[r]$ are zero, so we can eliminate them from the sum

$$Y[k] = \sum_{r=0}^{R-1} y[r] \cdot W_{R}^{rk} = \sum_{n=0}^{N-1} y[n \cdot M] \cdot W_{R}^{nMk} = \sum_{n=0}^{N-1} x[n] \cdot W_{R}^{nMk} $$

Next we take a look at the twiddle factor
$$W_{R}^{nMk} = e^{-j\frac{2 \pi nMk}{R}} = e^{-j\frac{2 \pi nk}{N}} = W_{N}^{nk} $$
 and put this back into  the Fourier Transform

$$Y[k] =  \sum_{n=0}^{N-1} x[n] \cdot W_{N}^{nk} $$

That's pretty much it. Comparing to the first equation we see directly that
$$Y[k]=X[k]$$

We also see that $Y[k]$ is periodic in $N$ and not just in $R$. Since both $x[n]$ and $W_{N}^{nk}$ are periodic in $N$, so must be $Y[k]$.


A basic question in FIR filter design using MATLAB I was reading from Internet about FIR filter design using MATLAB, and I saw parameters like passband frequency, stopband frequency, passband attuenation, stopband attuenation, etc. What made me confused was the parameter sampling frequency.

I don't understand the significance of sampling frequency with respect to a filter. Is it the sampling frequency of input signal that this filter can work on or is it something else?

All digital filter frequency parameters (passband begin frequency, passband end frequency, and stopband begin frequency) are stated in terms of an input signal sequence's Fs sampling frequency. For a lowpass filter example, if I said a lowpass filter's passband width (it's "cutoff" frequency) is 0.2, I'm saying that the cutoff frequency is 0.2 times Fs. So if an input signal's Fs was 1.0 kHz that lowpass filter's cutoff frequency is 200 Hz. However, if an input signal's Fs was 20 kHz then in this situation that same lowpass filter's cutoff frequency is 4 kHz.

For a different scenario, if you said, “Rick, please use Matlab and design a digital lowpass filter whose cutoff frequency is 300 Hz.” It would be impossible for me to do that unless you told me what was the filter’s input sequence’s Fs sample rate. But if you then said “Fs = 1000 Hz”, then I’d know the filter’s ‘normalized’ cutoff freq was 300/1000 = 0.3 enabling me to design the lowpass digital filter.


Are oversampling and downsampling the same thing? According the books I have on DSP,

is oversampling, sampling a signal at a much higher samplingrate then normally required, so to reduce the sample-rate (and increase the resolution of the signal) afterwards
is downsampling (decimation), lowering the samplerate of a signal, e.g. to reduce the CPU-requirements of signal-processing blocks later down the chain

So can I conclude that these two techniques are basically the same?

And does this mean that -as a concequence- downsampling will also increase the signal-resolution, just as with oversampling?

Or am I missing something?

Kristoff

OK. it looks like I got things mixed up.

The book on DSP I have mentions oversampling and that it is / can be followed by downsampling/decimation, ... and therefor I (wrongly) assumed that "oversampling" included both steps.

Mystery solved.
Thanks for the correction.

Kristoff


Determine filter type (HPF or LPF) by its coefficients I want to determine optimal initial states for a FIR/IIR filter to get rid of oscillations and shifts at the initial interval.

I've found that these initial states depend on filter type in the following way (up to a constant):

for k=1:order % (Direct Form II Transposed)
   z_hp(k) = sum(1*b(1+k:end)- 0*a(1+k:end)); % for HPF:  x=const, y=0
   z_lp(k) = sum(1*b(1+k:end)- 1*a(1+k:end)); % for LPF:  x=const, y=const
end

Then I should choose z_hp or z_lp - how can I do this given only filter coefficients?

By now the solution is empirical - passing 10 points of constant value to the filter and compare the shift in output:

x = ones(1,10);
y_hp = filter(b,a,x,z_hp);
y_lp = filter(b,a,x,z_lp);

if (1-y_lp) > y_hp
    z = z_hp;
else
    z = z_lp;
end

I wonder if I there is a simple formula for this.

I am not sure if my answer will correctly suits your question as it seems a little confusing. Anyway, I understand that you want to eliminate the initial transition of your filter. What I propose is to initialize your filter state so that the output equals the input at $t=0$

I had to solve that problem in the past. What I did was to assume that all delays were initialized to the same value. Doing that, you simplify your problem a lot. I don't have a lot of theory to back this technique, but that seems to work quite good on my side.

Let's start with Direct Form II

If we consider that all delays are set to the same value $K$, we find this output equation.

$$ y = (-a_1k-a_2k+x)b_0+b_1k+b_2k $$
A little algebra :
$$ y = b_0k(-a_1-a_2-...)+b_0x+k(b_1+b_2+...)$$
$$ y = -kb_0\sum_{n=1}^{N_a}{a_n} + b_0x+k\sum_{n=1}^{N_b}{b_n}$$

Now, since we know we want the output to be equal to the input, we can say $y=x$

$$ x = -kb_0\sum_{n=1}^{N_a}{a_n} + b_0x+k\sum_{n=1}^{N_b}{b_n}$$

Then

$$ x-b_0x = -kb_0\sum_{n=1}^{N_a}{a_n}+k\sum_{n=1}^{N_b}{b_n}$$
$$ x(1-b_0) = -kb_0\sum_{n=1}^{N_a}{a_n}+k\sum_{n=1}^{N_b}{b_n}$$
$$ \frac{x(1-b_0)}{-b_0\sum_{n=1}^{N_a}{a_n}+\sum_{n=1}^{N_b}{b_n}} = k $$

So there you go. Initialize all your delays with $K$ and your filter will start with an output equal to the input.

You can do this for other forms as well.

Hope that helps


Why does over-modelling an adaptive AR NLMS filter fix sharp spikes? I just simulated an auto-regressive second-order model fueled by white noise and estimated the parameters with normalized least-mean-square filters of orders 1-4.

As the first-order filter under-models the system, of course the estimations are weird. The second-order filter finds good estimates, although it has a couple of sharp jumps. This is to be expected from the nature of NLMS filters.

What confuses me is the third- and fourth-order filters. They seem to eliminate the sharp jumps, as seen in the figure below. I can't see what they would add, as the second-order filter is enough to model the system. The redundant parameters hover around $0$ anyway.

Could someone explain this phenomenon for me, qualitatively? What causes it, and is it desirable?

I used step size $\mu=0.01$, $10^4$ samples, and the AR model $x(t)=e(t)-0.9x(t-1)-0.2x(t-2)$ where $e(t)$ is white noise with variance 1.

The MATLAB code, for reference:

% ar_nlms.m
function th=ar_nlms(y,order,mu)
N=length(y);
th=zeros(order,N); % estimated parameters
for t=na+1:N
    phi = -y( t-1:-1:t-na, : );
    residue = phi*( y(t)-phi'*th(:,t-1) );
    th(:,t) = th(:,t-1) + (mu/(phi'*phi+eps)) * residue;
end

% main.m
y = filter( [1], [1 0.9 0.2], randn(1,10000) )';
plot( ar_nlms( y, 2, 0.01 )' );

What appears to be happening is, as you start over-modeling, the error signal becomes less and less white.

I modified your code to return the error signal (part of the residue term).

This plot shows the off-zero-lag coefficients of the xcorr of the error for order = 2 (blue), 3 (red), and 4 (green).  As you can see, the close-to-but-not-zero lag terms are getting larger in magnitude.

If we look at the FFT (spectrum) of the xcorr of the error, then we see that the lower frequency terms (that cause the large drifts) are getting smaller (the error is containing more high frequencies).

So it seems the effect of over-modeling in this case is to high-pass filter the error, which (for this example) is beneficial.

function [th,err]=ar_nlms(y,order,mu)
eps = 0.000000001;
N=length(y);
th=zeros(order,N); // estimated parameters
err = zeros(1,N);
for t=order+1:N
    phi = -y( t-1:-1:t-order, : );
    err(t) = y(t)-phi'*th(:,t-1);
    residue = phi*( err(t) );
    th(:,t) = th(:,t-1) + (mu/(phi'*phi+eps)) * residue;
    size(residue)
end


IIR filter implementation in Direct Form 2 I have implemented Direct Form 2 IIR filter. The input is the Kronecker delta function. I have written code for the response of the filter for Kronecker delta input. The code is:
clc;
clear all;
close all;
x=[1,zeros(1,150)];%input sequence

y=[];%output sequence

b=[1,-1.5511722144889886,1];%numerator coefficients

a=[1,-1.7473520798913149,0.79762351115754226];%denominator coefficients

bs=0.099825504845559299*b;
%...floating point.....%
w=[];

w(1)=x(1);

w(2)=-a(2)*w(1)+x(2);

y(1)=bs(1)*w(1);

y(2)=bs(1)*w(2)+bs(2)*w(1);
 for n=3:1:150

w(n)=x(n)-a(2)*w(n-1)-a(3)*w(n-2);

y(n)=bs(1)*w(n)+bs(2)*w(n-1)+bs(3)*w(n-2);

end

%freqz(y);

%

%

%.....fixed point conversion....%

%.....considering precision of 16 i.e q16 format....%

for n=1:3

    a_f(n)=int64(a(n)*pow2(16));%.....................converted into fixed......%

end

 for n=1:3

     b_f(n)=int64(bs(n)*pow2(16));

 end

for n=1:3

    a_fix2float(n)=double(a_f(n))*pow2(-16);

end

for n=1:3

    b_fix2float(n)=double(b_f(n))*pow2(-16);

end

w_f(1)=(x(1)*pow2(16));

w_f(2)=((-a_f(2)*w_f(1))*pow2(-16));

for n=3:150

    w_f(n)=((-a_f(2)*w_f(n-1))-(a_f(3)*w_f(n-2)))*pow2(-16)+x(n);

end

y_f(1)=(b_f(1)*w_f(1))*pow2(-16);

y_f(2)=((b_f(1)*w_f(2))+(b_f(2)*w_f(1)))*pow2(-16);

for n=3:150

 y_f(n)=((b_f(1)*w_f(n))+(b_f(2)*w_f(n-1))+(b_f(3)*w_f(n-2)))*pow2(-16);

 end

...converting fixed to float...%

% format long

y_fix2float(1)=double(y_f(1))*pow2(-16);

y_fix2float(2)=double(y_f(2))*pow2(-16);

for n=3:150

    y_fix2float(n)=double(y_f(n))*pow2(-16);%...convert to float..%

end

 for n=1:150

     w_fix2float(n)=double(w_f(n))*pow2(-16);

 end

   %to manually plot frequency response..%

  HH = abs(fft(y));

 yy=abs(fft(y_fix2float));

 plot(HH(1:50),'.r');

 hold on

 plot(yy(1:50),'g');

 xlabel('time');

ylabel('response');

legend('floating','fixed');

title('directform2 q16');

I have written MATLAB code in both floating and fixed point. In fixed point I considered
precision of q16.16. My response is:

My question is:  Is the response I got correct?

It seems to be wrong. The first sample of the impulse response is always equal the b0 coefficient (b(1) in Matlab, so in this case  that should be 0.099825504845559299 which doesn't seem to be the case.

The correct answer can be checked with yref = filter(bs,a,x); and it's substantially different from your answer. You have shown that the fixed point version matches the floating point version, however the floating point version is wrong to start with.

Finally, Direct Form II is a horrible choice for fixed point filters. The transfer function between input and state variables is given by the pole-only transfer function. Even in this harmless example this is about 27 dB of gain. It's almost impossible to prevent clipping in Direct Form II. Use Direct Form I or Transposed Form II instead


Why is $y(t)=x(t/2)$ a non-causal system? I was going through my signal and system notes.they say $y(t)=x(t/2)$ is a non causal system? As non causal system depend on future inputs. how $t=t/2$ is future value of time? i could not understand how? could somebody please explain it to me?

Because for negative values of $t$ you have, for example, $y(-2) = x(-1)$ which depends on a future value of $x(t)$ at $t=-1$ for the current value of $y(t)$ at $t=-2$.

Note that $t=-1$ represents a future time for $t=-2$.


Detect circles in image I have an image that has pixelated circles in it (a cross section of wires). I want to identify the wires, their location and count them.

The circles never overlap and the wires are all within a 95% radius. However the spacing might not be exact as there can be gaps between the wires.

What is a good way to find the geometric mean of the group of pixels that represent a wire?

The above picture is what I would like to do, to fit a circle to a group of pixels (which I did using paint). I would think that the first step is to find the geometric mean of the group.

I did try using a tophat kernel and threshold to emphasize the centers, but I end up with a few pixels and no center pixel. I need to identify the center and end up with only one pixel which is difficult to do with thresholding since the pixel values differ in brightness.

Raw picture:

1) Normalize your image to range $[0,255]$.

2) Select a threshold and threshold the image. For your image, what worked is: $\tau=[140-150]$.

3) Compute a Euclidean distance transform.

4) Apply watersheds segmentation.

If I apply this procedure, here is what I get:

Not perfect, but maybe a good start. The result looks similar to performing a Voronoi diagram on a smoothed/noise filtered version of the image.

Now to each region, one could fit a circle, using any algorithm (edge extraction and fitting or Hough transform etc.) For now, I used a simple circle fitting to the edges:

Again, not the best, but could be a sufficient starting point.

I also think that you would find the following useful:

  Segmentation Methods for Light Microscope Image

  Detecting Spheric Particles


How to find the Freqs bins Ok, so I'm kind of moving forward with implementing the STFT. I realise that I was taking the entire FFT and not the first half.

The problem that I'm having now is, understanding how to get the Freq bins, what calculation should I use? For example, I have been using this formula:

(10 * log10(sqrt(v[i][j].re * v[i][j].re + v[i][j].im * v[i][j].im)));

Which, I realise that this is for decibels. Currently, my graph looks like the following:

I know this makes sense to me now, this is finally starting to come together!

Would this algorithm work?:

I take the magnitude of each of the elements inside the bin `mags = (re * re + im + im) I then find the maximum index of the bin, multiply it by the Fs and divide by the bin size?

If I did this, however, would this only give me 1 result for each block?

Your Y axis likely should be scaled to frequency (Hz), not FFT bin index [0..127].  The frequency scale would be something like the bin_index * sample_rate / FFT_length .


Discrete Cosine Transform (DCT) Coefficient Distribution I have two images :

Original Image

Binarize Image

I have applied Discrete Cosine Transform to the two images by dividing the 256x256 image into 8x8 blocks. After, I want to compare their DCT Coefficient Distributions.

import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
import numpy as np
import os.path
import scipy
import statistics

from numpy import pi
from numpy import sin
from numpy import zeros
from numpy import r_
from PIL import Image
from scipy.fftpack import fft, dct
from scipy import signal
from scipy import misc

if __name__ == '__main__':
    image_counter = 1

    # Opens the noisy image.
    noise_image_path = 'noise_images/' + str(image_counter) + '.png'
    noise_image = Image.open(noise_image_path)

    # Opens the binarize image
    ground_truth_image_path = 'ground_truth_noise_patches/' + str(image_counter) + '.png'
    ground_truth_image = Image.open( ground_truth_image_path)

    #Converts the images into Ndarray
    noise_image = np.array(noise_image)
    ground_truth_image = np.array(ground_truth_image)

    #Create variables `noise_dct_data` and `ground_truth_dct_data` where the DCT coefficients of the two images will be stored.
    noise_image_size = noise_image.shape
    noise_dct_data = np.zeros(noise_image_size)
    ground_truth_image_size = ground_truth_image.shape
    ground_truth_dct_data = np.zeros(ground_truth_image_size)

    for i in r_[:noise_image_size[0]:8]:
         for j in r_[:noise_image_size[1]:8]:
              # Apply DCT to the two images every 8x8 block of it.
              noise_dct_data[i:(i+8),j:(j+8)] = dct(noise_image[i:(i+8),j:(j+8)])
              # Apply DCT to the binarize image every 8x8 block of it.
              ground_truth_dct_data[i:(i+8),j:(j+8)] = dct(ground_truth_image[i:(i+8),j:(j+8)])

The above code gets the DCT of the two images. I want to create their DCT Coefficient Distribution just like the image below:

The thing is I dont know how to plot it. Below is what I did:

#Convert 2D array to 1D array
noise_dct_data = noise_dct_data.ravel()
ground_truth_dct_data = ground_truth_dct_data.ravel()

#I just used a Histogram!
n, bins, patches = plt.hist(ground_truth_dct_data, 2000, facecolor='blue', alpha=0.5)
plt.show()

n, bins, patches = plt.hist(noise_dct_data, 2000, facecolor='blue', alpha=0.5)
plt.show()

image_counter = image_counter + 1

My questions are :

What does the X and Y-axis in the figure represents?
Are the value stored in noise_dct_data and ground_truth_dct_data, the DCT coefficients?
Does the Y-axis represents the frequency of its corresponding DCT coefficients?
Is the histogram appropriate to represent the DCT coefficient distribution?
The DCT coefficients are normally classified into three sub-bands based on their frequencies, namely low, middle and high frequency-bands. What is the threshold value we can use to classify a DCT Coefficient in low, middle or high frequency band? In other words, how can we classify the DCT coefficient frequency bands radially? Below is an example of the radial classification of the DCT coefficient frequency bands.

The idea is based from the paper : Noise Characterization in Ancient Document Images Based on DCT Coefficient Distribution

Probably the value of the DCT coefficient and its frequency (i.e., a histogram plot). However, it's impossible to say that with certainty given the information you have provided. If you found that plot in the linked paper, then the definitions of the x- and y-axes are in all likelihood given therein. The paper is behind a paywall, so it's not helpful to anyone without an IEEE account.
Yes. Check the scipy documentation for the dct() function.
See 1.
This is impossible to answer. What does appropriate mean in this context?
To classify a coefficient as either low, middle, or high frequency, compute its distance from the DC bin. For bin (m,n), its distance will be r = sqrt( (m-1)^2 + (n-1)^2 ). Judging by the picture you included, it looks like low frequencies are those with r < M/2, middle frequencies are those with M/2 <= r < M, and high frequencies are those with r >= M.

EDIT:

From the linked paper:

  Furthermore, it has been explained that noise addition to an image
  changes its DCT coefficient distribution [9], as we confirm in Fig. 2
  for a document image. Part (a) of the figure is an example of
  noiseless document image, part (b) is the noisy version of image (a)
  due to Gaussian blurring and Gaussian noise addition, and part (c)
  shows the higher frequency DCT coefficient distribution of of part(a)
  and (b) respectively. A lower-peaked histogram belongs to the
  ground-truth (GT) image, while a sharp-peaked histogram belongs to the
  noisy document image.

So part (c), which is the image you included, shows the histograms of the "higher frequency" DCT coefficients. Thus, the x-axis is the value of the DCT coefficient for the high frequency band, and the y-axis is its normalized histogram value. (The normalized histogram value is sometimes called "frequency", but that can get confusing in this case since DCT coefficients also represent frequencies. The normalized histogram value is just the number of occurrences of a certain value divided by the total number input samples to the histogram.)


How do I efficiently evaluate a convolution integral between different-sized domains? This is a question that I've previously asked over on math.stackexchange, and I have yet to receive a useful answer. It was suggested that I post this here.

The problem itself originally comes from simulation of wave propagation, and involves large two-dimensional arrays of complex data. I'll explain it here as a simplified one-dimensional problem, however.

Suppose I have the following integral relationship involving a convolution:
$$
f(x) = \int_{-\infty}^{+\infty}dx'\, g(x')\, K(x - x')
$$
Assume that the functional form of the kernel $K$ is known, and suppose that the function $g$ has compact support on $[-a/2, +a/2]$ for some width $a>0$. Assume further that I have data about the function $g$ in the form of $N$ equally-spaced samples $g_n$:
$$
g_n = g(x_n)\, ,
$$
where
$$
x_n = -\frac{a}{2} + n\frac{a}{N}\, , \qquad n = 0,..., N-1
$$
My goal is to generate some $N$ evenly-spaced samples of the function $f$.

If we suppose that we want to evaluate $f$ on the same domain on which $g$ lives, and for which we have samples of $g$, then it's easy. I discretize everything in the usual way, and the whole thing becomes a discrete convolution:
$$
f_n = f(x_n) = \sum_{m = 0}^{N-1} g_m K_{n - m}\, ,
$$
where
$$
K_{n - m} = K(x_n - x_m) = K\left(\frac{a}{N}(n-m)\right)\, .
$$
Everything is simple since I can use the discrete convolution theorem: I FFT the $\{g_n\}$ and $\{K_n\}$ sequences, multiply them, do an inverse FFT, and I'm done. I have my $\{f_n\}$ data. The whole thing will be an $O(N\log N)$ process.

Suppose instead that I have to evaluate $f$ on a different-sized domain. For example, perhaps the $g$ data is narrow, but the kernel $K$ is very wide. Let's say I have to generate $N$ samples of $f$ on the domain $[-L/2, +L/2]$, for some domain width $L > a$. $L$ could even be much larger than a.

If I try to discretize everything in a similar way, I'd have to first define a sequence of points $\{X_n\}$ in the wider domain:
$$
X_n = -\frac{L}{2} + n\frac{L}{N}\, , \qquad n = 0,..., N-1
$$
Then the data I seek is $\{f_n\}$ where $f_n = f(X_n)$. The original integral is still over the smaller domain, so we have:
\begin{align}
f_n
&= \int_{-a/2}^{+a/2}dx'\, g(x')\, K(X_n - x') \\
&\approx \sum_{m = 0}^{N-1} g_m\, K(X_n - x_m) \qquad \text{(discretize)}\\
&= \sum_{m = 0}^{N-1} g_m\, K\left(\frac{L-a}{2} + \frac{a}{N}(\alpha\, n - m)\right)\, ,
\end{align}
where $\alpha = L/a > 1$ is the ratio of the two domain sizes.

This is no longer a discrete convolution, since the second term is no longer a function of $n-m$, but instead of $\alpha\, n - m$. So I can no longer use the convolution theorem. I can imagine evaluating the sum above by brute force, but that would be an $O(N^2)$ process, which is prohibitive for my problem.

Is there any way of efficiently calculating the discrete version of this convolution integral between two different domain sizes? Thanks.

Essentially, to my understanding this problems boils down to a fast-convolution problem:

You have a function $g(x)$, which has a small domain and you have a function $K(\tau)$, which has a big domain and you want to calculate

$$f(x) = \int_{-\infty}^{\infty}g(x')K(x-x')dx'$$

in a discretized fashion, where you want to calculate as efficient as possible. Especially, you are interested in $f(-L/2+nL/N)$ for $n=0,\dots,N-1$.

When considering the discretization you need to take care of the bandwidth of $g$ and $K$ and you need to sample with the according Nyquist frequency at least. Let's assume that $g$ can be sufficiently sampled with $g(-\frac{a}{2}+n\frac{a}{N})$ (i.e. the sampling is fast enough regarding Nyquist frequency). Let's also assume, that $K(\tau)$ has the same or smaller bandwith, i.e. the sampling rate for $g$ is also sufficient for $K$.

Then, let's further assume $L=qa, q\in\mathbb{N}$, i.e. L is a multiple of $a$. Define the following functions:

$$g[n]=g(-\frac{L}{2}+n\frac{L}{qN}), \quad n=0,\dots,qN-1$$
and
$$K[n]=K(-\frac{L}{2}+n\frac{L}{qN}), \quad n=0,\dots,qN-1$$

Note, that $g[n]$ is zero for most $n$, but that is fine here.

Now, you have

$$f[n]=\sum_n'g[n']*K[n-n'], \quad n=0,\dots,qN-1$$

Which is a higher-resolution version of your original $f_n$ (It has qN instead of only N samples within the interval $\pm\frac{L}{2}$). So, what you actually want to have is a downsampled version of this $f[n]$ given by

$f_n=f[qn]$ (where $f_n$ is the $f_n$ from your post, i.e. with $n=0,\dots,N-1$.

So, this is also a discrete convolution problem, and it can be solved via FFT (note that the FFT actually considers circular convolution, if you don't zero-pad). However, since your $g[n]$ has a much smaller domain that $K[n]$, you can resort to fast-convolution algorithms (Overlap-save or Overlap-add).

To summarize: The trick is actually just to sample both signals with the same sampling frequency (accepting that $g[n]$ contains a lot of zeros) and then understanding this problem as a simple discrete convolution.


Testing whether a process is a Wiener process Ideally I would like links to code implementations (eg. Matlab ) or book references, but I would appreciate suggestions on various methods.

We start with sampled process $X_{t}$.

A straightforward way is through repeated normality tests:

a. Test for normality of $X_{t}$ using KS test

b. Independence of increments: we test whether the product $X_{T/2}(X_{T}-X_{T/2})/\sigma^{2}$ is close to a standard normal, where $\sigma^{2}=(T/2)^{2}$ .

c. Joint normal for increments $Y_{n}=(X_{t_2}-X_{t_1},...,X_{t_n}-X_{t_n-1})$: we divide the possible values of $Y_{n}$ into $m$ cells and denote $O_{j}:=$# samples that fall into cell $j$. Then the statistic $\sum \frac{(O_{j}-E[O_{j}])^{2}}{E[O_{j}]}$ should be approximately a $\chi^{2}-$distribution.
Maybe testing for quadratic variation=$t$ and martingale property (Levy characterization)?
Some spectral characterization for WP?
The other link is testing for fractional WP, which is more involved.
Turning it into a test for White noise since it is the "derivative" of WP. So maybe taking finite difference for WP and doing white noises tests for it:
$$(X_{t+\Delta t}-X_{t})/\Delta t.$$

That's not a good question to ask (you are ok and I'll vote up your question in a sec).
BM induces a measure on functions, called the wiener measure. Since this is a measure on functions, it is infinite dimensional. The thing with measures on infinite dimensional spaces is that they tend to be mutually singular: events that happen w.p. 1 wrt one measure may have zero probability wrt the other. For example, take BM $B_t$ and BM with some tiny constant added $C_t= B_t + \epsilon$.
Then $P(B_0 = 0)=1$ but $P(C_0=0)=0$. This situation is fairly generic.

What this means to you is that unless your process really is exactly BM (highly unlikely), as you take more and more data points your measure will converge to something that is mutually singular to BM and your hypothesis test will give a "no". This may happen even if your process is pretty "similar" to BM.

That's why this is not a good question to ask. Measures over function spaces are tough to work with.


The error probability of 16QAM The signal constellation for a communication system with 16 equiprobable symbols is shown as below. The channel is AWGN with noise power spectral density of $N_ 0/2.$Using the union bound, ﬁnd a bound in terms of $A$ and $N_0$ on the error probability for this channel.

Does the union bound means error probability?the solution said
$$P_e \le 15Q\left(\sqrt{\frac{d^2_{min}}{2N_0}}\right)=15Q\left(\sqrt{\frac{2A^2}{N_0}}\right)
$$

but why is

\begin{align}
&2\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)+3\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)+3\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)+2\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)\\
&+3\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)+4\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)+4\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)+3\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)\\
&+3\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)+4\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)+4\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)+3\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)\\
&+2\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)+3\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)+3\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)+2\cdot\frac{1}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)\\
&=\frac{48}{16}Q\left(\frac{A}{\sqrt{\frac{N_0}{2}}}\right)\quad{?}
\end{align}

In general, the symbol error probability is given by

$$P_E=\sum_{i=1}^MP[s_i]P[E|s_i]\tag{1}$$

where $P[s_i]$ is the probability that symbol $s_i$ was sent, $P[E|s_i]$ is the error probability given that $s_i$ was sent, and $M$ is the number of symbols. If all symbols are equally likely, $(1)$ simplifies to

$$P_E=\frac{1}{M}\sum_{i=1}^MP[E|s_i]\tag{2}$$

For $16$-QAM it is possible to derive an exact expression for the error probabilities $P[E|s_i]$ (and hence for the symbol error probability $P_E$). You can find the derivation in most digital communication text books. A good approximation (better than the union bound) is given in this answer.

The union bound approximation gives an upper bound on $P[E|s_i]$:

$$P[E|s_i]\le\sum_{j\neq i}P[s_j\text{ chosen}|s_i\text{ sent}]\tag{3}$$

Under the given assumptions we have

$$P[s_j\text{ chosen}|s_i\text{ sent}]=Q\left(\sqrt{\frac{d^2_{ij}}{2N_0}}\right)\tag{4}$$

where $d_{ij}$ is the distance between symbols $s_i$ and $s_j$. With $(3)$ and $(4)$ the union bound approximation of $(2)$ becomes

$$P_E\le\frac{1}{M}\sum_{i=1}^M\sum_{j\neq i}Q\left(\sqrt{\frac{d^2_{ij}}{2N_0}}\right)\tag{5}$$

Note that generally you will get a sum of $Q$ functions with different arguments because the distances $d_{ij}$ are usually not all the same. This is also the case for $16$-QAM.

The solution you provided is still more conservative (hence, worse) than the union bound. That solution uses the fact that the $Q$ function is monotonically decreasing, which implies

$$Q\left(\sqrt{\frac{d^2_{ij}}{2N_0}}\right)\le Q\left(\sqrt{\frac{d^2_{min}}{2N_0}}\right)\tag{6}$$

with $d_{min}=\min_{i\neq j}\{d_{ij}\}$ being the smallest occurring distance between symbols. Plugging $(6)$ into $(5)$ we finally get the bound

$$P_E\le (M-1)Q\left(\sqrt{\frac{d^2_{min}}{2N_0}}\right)\tag{7}$$

which, as mentioned above, is even less tight than the union bound.


explanation of correlation of stationary stochastic processes I have some doubts about correlation in stationary stochastic processes.

I know that the expectation of a random variable is
$$E(x)=\int_{-\infty}^{+\infty} a f_x(a)da$$
and the mutual correlation of two variables is $E(x_1 x_2)$.

In two stationary stochastic processes the mutual correlation should be $E(x_1(t+\tau)x_2(t))$ now when and why I can translate this formula in:
$$
\int_{-\infty}^{+\infty} x_1(t+\tau)x_2(\tau)d\tau
$$
I have not understand where this last formula comes from.

You are mixing up two different notions.

Your random process is a collection of random variables $\{X(t)\colon -\infty < t < \infty\}$, one random variable for each time instant. The autocorrelation function of the process is
$$R_X(t, \hat{t}) = E[X(t)X(\hat{t})], -\infty < t, \hat{t} < \infty$$
and is a two-dimensional function (meaning it has two arguments,
and if you try to visualize it, you should imagine a surface above
(and possibly
below) the plane with coordinate axes $t$ and $\hat{t}$. Definitely above
because $R_X(t,t) = E[(X(t))^2]\geq 0$ for all $t$ (and processes for
which $R_X(t,t) = 0$ for all $t$ are kinda boring).

For a stationary process, $R(t,\hat{t})$ depends only on the
difference $\hat{t}-t$ of the two arguments, and it is common
practice to write it as a function of one variable $\tau = \hat{t}-t$.
Thus,
$$\mathcal R(\tau) = R(t,\hat{t}) = E[X(t)X(t+\tau)].$$

The other notion is that if the experiment is performed, then
each random variable takes on a numerical value. Thus, $X(t)$ takes
on a specific numerical value that we name as $x(t)$ while
$X(\hat{t})$ takes on some other numerical value that we name
as $x(\hat{t})$ and so on. Thus, we have a real-valued function
(or signal) $x(t)$ (nothing random about that!) for which the
autocorrelation function is defined as
$$R_x(\tau) = \int_{-\infty}^\infty x(t)x(t+\tau)\,\mathrm dt.$$
Note that $x(t)$ is called a sample path or realization
of the random process: it is what you observe on your oscilloscope
if you lift up your head from the math a little bit
and look at the real world.

It is every DSPer's fondest hope that $R_x(\tau)$ is the same function
as $\mathcal R(\tau)$. Unfortunately, nothing in what has been
said above justifies this hope. But, since people like to assume
that it is indeed true that $R_x(\tau) = \mathcal R_X(\tau)$, mathematicians have come
up with a notion called ergodicity which (for our purposes) says that
for an ergodic process, the ensemble means (expectations)
equal the corresponding time-averages.  As a practical matter,
all model-builders assume that the processes that they want
to model are indeed both stationary and ergodic,
and decide on what the mean and
the autocorrelation function etc should be by making observations of
actual sample paths etc.

For more than what you probably want to know about random processes,
see this answer of mine, and
for an example of a (boring) stationary random process that is
not ergodic, see this other answer of mine.


How to do sliding window multiplication in frequency domain I want to perform the following computations with the impulse response of length L of a LTI system:

Multiply the impulse response with a window w(n) of length M << L.
Take the FFT of the result in 1 and save as a vector Y2.
"Move" the window to the "right" by M/2 samples, i.e. take w(n-M/2), and multiply the impulse response with this new window.
Take the FFT of the result in 3 and save as a vector Y2.
Move the window to the right again by M/2 samples, i.e. take w(n-2*M/2), and multiply the impulse response with this new window.
Take the FFT of the result in 5 and save as a vector Y3.
Etc... (keep doing this process until the end of the impulse response, I hope you get the idea).

I think this is easy to do with a for loop in Matlab. However, I already have immediate access to the frequency response (the FFT of the impulse response). I was therefore wondering if there is a way to get my vectors Y1, Y2, Y3, etc... by computations in the frequency domain instead? I know that multiplication in the time domain corresponds to convolution in the frequency domain. So to get for instance the vector Y5 with the above notation, would it be possible to do as follows?
Take the FFT of the window w(n-5*M/2), and simply convolve it with the frequency response? And how do I know if it's faster/slower than just taking the inverse FFT of the frequency response and doing the steps described above?

Is there a way to get my vectors Y1, Y2, Y3, etc... by computations in the frequency domain instead?

Yes. Rough outline: you take a full length FFT if the zero padded window. For each frame you can implement the time shift a multiplication with $e^jn\omega$ and then do a circular convolution and downsample. If you choose the window size well, you can do a down-sampled convolution instead.

And how do I know if it's faster/slower than just taking the inverse FFT of the frequency response and doing the steps described above?

It will almost always be slower. Convolution has computational complexity of $N^2$ whereas multiplication has only $N$. Saving a few short FFTs will not compensate of this.


Processing an ECG signal with a median filter I have read in a couple of papers that the noise from an ECG signal can be removed via median filter. One such example I found was on stackoverflow, where multiple methods were suggested and one of them being the median filter. The following image is taken from the post on stackoverflow.

The main idea of the median filter is to run through the signal entry by entry, replacing
each entry with the median of neighboring entries. The pattern of neighbors is called the
"window", which slides, entry by entry, over the entire signal.

What I do not understand is, why doesn't the QRS complex get removed as well? When I tried it using matlab, the bigger I set the window, the better the signal got filtered, should't it be the opposite?

You are certainly doing something very wrong. You should upload your data in order to get better responses. You can upload your data into any upload sites you wish and provide its link here, practically.

Median filter is a highly highly nonlinear filter (it re-orders the sample positions!). The output of the median filter at a position $n$ is the median of the values that reside in the window scope; i.e., it's the value that resides in the middle when the samples are sorted in order. Hence median filtering requires sorting for each computation. This makes it quite slow as well (a deailed answer actually depends on the architecture...)

Median filter is mainly used for speckle or salt and pepper noise removal, in essence these are local noise samples whose frequency domain filtering is not possible without degrading the whole signal. Such local (in time) peaks will have wide band frequency spectrums which inhibit frequency domain attacks to remove them, therefore, leaving only the time domain (or time-frequency domain) approaches possible.

Median filter has a tendency to preserve edges, therefore quite preferred in certain image enhancement operations. However it also has the side effect of washed out results (texture details are lost and only strong edges remain) which indicates that their use should be performed with care.

In principle, the longer the window size, the stronger will be the washed out effect. So it's customary to use as short as possible window sizes (unless otherwise dictated by the particular application)

Coming to your example plot. Using a window size larger than 50 samples will wash out local details that might be important to you, so you should use a window size less than 50, I guess (based on the plot you provided).


What is the relationship between angular frequency and normalized angular frequency This is a slide from my lecture notes:

My professor used the following words " We denote digital frequencies with capital letters and analogue frequencies with lower case letters"

The problem I have is $$\Omega=\omega T_s$$

Wherever I am looking it is actually the other way round i.e.

$$\omega = \Omega T_S$$

For example  here

Is my proffesor wrong?

Again there is no wrong or right here. In the Alan Oppenheim's Discrete-Time Signal Processing book, the notation is as follows:

when there are only continuous-time signals we use $\omega$ for radians per second frequency.
when there are only discrete-time signals we use $\omega$ for radians per sample frequency
when both types of signals are present, (as in sampling), we use $\Omega$ for the continuous-time radians per second frequency, and $\omega$ for the discrete-time radians per sample frequency.

Furthermore the relationship between the two frequencies because of sampling normalization is:

$$ \Omega = \frac{\omega}{T_s} $$

or equivalently

$$ \omega = \Omega ~ T_s $$

However, your instructor seems to prefer the opposite notation for frequencies. That's why most literature seems to be reciprocal of your instructor's.


Reconstruction of a windowed signal (specifically Hamming windowed) Wondering about the resconstruction of the input signal when an overlapped windowing function has been applied.

I'm specifically applying a Hamming window with 50% overlap and I wish to reconstruct the input signal.

Here http://www.dsprelated.com/showthread/comp.dsp/105321-1.php it's given that summing the 50% overlapped windows gives a sum of 1.08.

In order to reconstruct the input signal (with the same exact amplitudes), then can I merely subtract the extra ~".08" from all the samples? Or, equivalently, divide by 1.08?

The overlapp-add does not result into the exact same signal though, but there's some miniscule variation for some reason. By testing with an array of 1.0s I get back varying values that are of the form 0.999 ...

I guess reconstructing doesn't always lead to the same exact input signal, depending on the windowing function and the amount of overlap?

did you ask this at comp.dsp?  appears so.  if you do what hotpaw says, your Hamming window will be

$$ w(t) =
\begin{cases}
0.426 \cdot \cos(\pi t) + 0.5 & \text{if }|t| < 1 \\
0.037 & \text{if }|t| = 1 \\
0 & \text{if }|t| > 1
\end{cases} $$

still a Hamming window, but it overlap-adds to 1 (with 50% overlap), rather than 1.08.


How do I quantify the amount of information being lost by filtering a digital signal? Let's say I'm doing feature extraction from some sensor reading using a DTW algorithm. I want to extract and catalogue a set of features.

Before I send my signal over to the DTW extractor, I have several tools to remove what I consider to be noise from the signal. I can use bandpass/notch filters or DFTs.

When doing this interactively, it's easy for me as a human to spot which part of the spectrum is interesting and which part of it is low entropy noise.

For instance, the following signal:

has three signals superimposed.

a low frequency square pulse (400 samples long - 2.5mHz)
a mid range sine wave (33 samples long - 30mHz)
and a high frequency noise.

To me, this is trivially visible. I say visible, but this is also audible  if I were to listen to the signal.

The ground noise happens to be the RF interference that a brushless fan makes when it sporadically comes on. So I could have the noise there or not there depending a completely random variable.

It's not clear whether the middle frequency sine wave is simply a superposition on the square wave, or if it is being generated by the same fundamental underlying phenomenon. If it is not the same phenomenon, there's a possibility that their respective phases will drift, and the rising edge of the square signal will only sometimes align with the rising edge of the sine wave. This means that for what ends up being a square wave signal, I could have dozens and dozens of different looking DTW features, and this is something I'd rather avoid.

Having described the above scenario, my questions are as follows:

what metric can I use (or even theory or branch of analysis if the answer isn't a simple one) to determine the amount of filtering necessary? For any given analysis window of say 500 samples, how do I determine programatically how to preprocess my signal before trying extract features?
given a sample window with a particular spectral distribution, does an optimal cutoff for a bandpass filter exist? (or a notch frequency)
if so, does it always exist?
if an optimal value exists, how do I obtain it?
what is the relationship between this and the amount of information in the signal?

Notes: the vagueness of my question is part of my problem. I don't know where to start looking for answers.

To be clear: I'm not attempting to get more information than there is. I'm wanting to degrade my signal as little as possible while removing obvious low entropy noise.

I've looked into things like adaptive filters, but these seem to often time have a separate measurement to aid in the feedback loop.

what metric can I use (or even theory or branch of analysis if the answer isn't a simple one) to determine the amount of filtering necessary?

You can use what you already have, that is, the "strengths" of the components that result from the analysis. You can sort the coefficients in decreasing order and then obtain their cumulative sum. Finally, you can set a threshold and stop adding components once the cumulative strength of the chosen components goes (for example) above 90%. In this way you would be selecting the $N_{sc}$ strongest components and rejecting components that have higher signal-to-noise ratios (or in other words are closer to some background noise level). If this has to be a low pass filter, then after you have isolated the first $N_{sc}$ components that contribute to 90% (for example) of the total signal strength then you can set the low pass filter cut off frequency at the component with the highest frequency (or bin).

  For any given analysis window of say 500 samples, how do I determine programatically how to preprocess my signal before trying extract features?

Please see previous answer, given the frequency (or scale) limits imposed by the length of the window (here 500 samples), you might have to adjust thresholds and cut offs. Please note that if you try to apply the above algorithm per frame (of 500 samples each), then, inevitably, you will get some sort of modulation between frames which, depending on your application you will have to manage somehow. This is because the signal strength will be changing between windows of 500 samples, which means that the potential cut off of the low pass filter will be changing as well, which means that each frame will bear different amounts of filtering. The easiest way to deal with this would be to increase the window overlap (at the expense of more frames and therefore more calculations).

  Given a sample window with a particular spectral distribution, does an optimal cutoff for a bandpass filter exist? (or a notch frequency)

Only under specific assumptions about the spectral distribution of the noise. In the absence of this, you can study your signal and try to create a Matched Filter that is tailored to the spectral characteristics of the signal and therefore anything outside that would be filtered out.

  if so, does it always exist?

Please see previous answer

  if an optimal value exists, how do I obtain it?

Please see previous answer

  What is the relationship between this and the amount of information in the signal?

This "Amount of information" is a bit unclear. Information has a specific meaning, so are you refering to the amount of bits that are expected to be able to encode completely your "clean" waveform? Or entropy derived by the spectrum maybe? Or spectral flatness? Or, are you talking about the actual information that this pulse sequence would carry if it was to be seen as a communications signal (?). As a general comment, the low pass filter would lower the entropy of the signal in its output...but this is too general a statement to be useful in your context perhaps (?).

Hope this helps. Happy to ammend the response if more data becomes available.


How Exactly Does MATLAB Zero Pad Signal? I am doing some massive number crunching in MATLAB which involves millions of PSD estimations. Each data segment has length 41. So I have been using the multi-taper method with nfft=41 (details below in the code). My last simulation took more than eight hours. The code is optimized and vectorized as I can make it so I thought maybe if I change nfft=64 that might speed up. The profiler confirms that pmtm is the slowest part by far.

So, looking at the MATLAB help pages, I see that in the documentation it says nothing more than if nfft is larger than the length of the data segment then it is zeropadded. My question is what exactly is MATLAB doing? How exactly is the signal padded?

I can't really find anything googling and I ran my own experiments and I am not getting what I'm supposed to. Here I present my code and the result.

close all
clear all

dseg = [ 395.9640
  392.3630
  388.8230
  385.2820
  381.6130
  378.2600
  374.7470
  371.1860
  367.9920
  364.7830
  361.6630
  358.5660
  355.5620
  352.5820
  349.5530
  346.8110
  343.8130
  340.9010
  338.1460
  335.5110
  332.6300
  330.0470
  327.4020
  324.0930
  321.4680
  318.8760
  316.4470
  314.0850
  311.5310
  309.0600
  306.7570
  304.6690
  302.5820
  300.3890
  298.2480
  296.3650
  294.1660
  291.7610
  290.3300
  287.9700
  285.9850];

% Original signal with length=41
[ffty1,fftf1] = pmtm(dseg,4,41,1/30);
semilogy(fftf1,ffty1,'-o')
hold on

% Matlab "zeropads" to length=64
[ffty2,fftf2] = pmtm(dseg,4,64,1/30);
semilogy(fftf2,ffty2,'-or')

% Original zeropaded on both sides
dsegpad = [zeros(11,1); dseg; zeros(11,1); 0];
[ffty3,fftf3] = pmtm(dsegpad,4,64,1/30);
semilogy(fftf3,ffty3,'-ok')

% Original zeropadded at the end
dsegpad = [dseg; zeros(23,1)];
[ffty4,fftf4] = pmtm(dsegpad,4,64,1/30);
semilogy(fftf4,ffty4,'-og')

% Original signal just repeated
dsegpad = [dseg; dseg(1:23)];
[ffty5,fftf5] = pmtm(dsegpad,4,64,1/30);
semilogy(fftf5,ffty5,'-oc')

I compare the original PSD with MATLAB's padded PSD and they are close. Then I zeropad myself instead of letting MATLAB do it. I zeropadded on both sides of the signal. I tried zeropadding all on one side. I even tried periodically repeating the signal to pad the length. But in all three of the cases my padded PSDs are way larger than the original and MATLABs PSDs. Here is what I see.

I tried looking at pmtm code but got lost pretty fast. Anybody knows what MATLAB is doing and why does it obfuscate it so?

Lets say you have a vector $ x = {\left[ 1, 2, 3, 4 \right]}^{T} $.
You want to have a look on its DFT transform then you apply DFT on it and have the 4 points DFT transform of the data.
In MATLAB it would be something like:
vXDft = fft(vX);

Yet, what happens if you want the 64 point DFT of the same data?
In MATLAB it would be something like:
vXDft = fft(vX, 64);

So what happens?
MATLAB just adds zero elements at the end of the vector until its size is 64:
vXDft = fft([vX; zeros([60, 1])]);

This will match what happens above (I assumed Column Vector).
In practice you get a denser grid in the Frequency Domain (64 Points).
Where the points are interpolated by Dirichlet Kernel (The Discrete equivalent of the Sinc function).


Can deep learning be applied to video compression? H.265 is based on incremental gains over time. The paradigms have largely remained the same. What are the odds of deep learning providing huge performance gains? Here is the main metric:

Perceptual quality remains the same at 2x improvement in bit-rate.
Assume unlimited compute power at encoder.
No need for speed.

Of course, I understand this is a really ill-posed problem because there are so many factors involved. But if it can be done, what would give the performance gain? Is it the following?

Larger CTUs instead of the 64x64 one in H.265. Perhaps even go towards object-based coding.
Better block matching algorithms that map longer dependencies.
Perceptual loss for block matching instead of MAD or MSE.
Etc...

[EDIT: addition of a March 2017 preprint]

Deep learning already has many applications in video, like enhancement (Deep Convolutional Neural Network for Decompressed Video Enhancement) or semantic analysis.

Recently, there have been some announcements related to video compression, for instance:

Generative Compression (Preprint, March 2017)

  Traditional image and video compression algorithms rely on
  hand-crafted encoder/decoder pairs (codecs) that lack adaptability and
  are agnostic to the data being compressed. Here we describe the
  concept of generative compression, the compression of data using
  generative models, and show its potential to produce more accurate and
  visually pleasing reconstructions at much deeper compression levels
  for both image and video data. We also demonstrate that generative
  compression is orders-of-magnitude more resilient to bit error rates
  (e.g. from noisy wireless channels) than traditional variable-length
  entropy coding schemes.

Video: Deep Learning reconstructing Blade Runner for new video compression techniques
Wyner-Ziv Video Coding using Hadamard Transform and Deep Learning

How important is deep learning here, and what performances are obtained is not clear to me yet.