What is the distinction between ergodic and stationary? I have trouble distinguishing between these two concepts. This is my understanding so far.

A stationary process is a stochastic process whose statistical properties do not change with time. For a strict-sense stationary process, this means that its joint probability distribution is constant; for a wide-sense stationary process, this means that its 1st and 2nd moments are constant.

An ergodic process is one where its statistical properties, like variance, can be deduced from a sufficiently long sample. E.g., the sample mean converges to the true mean of the signal, if you average long enough.

Now, it seems to me that a signal would have to be stationary, in order to be ergodic.

And what kinds of signals could be stationary, but not ergodic?
If a signal has the same variance for all time, for example, how could the time-averaged variance not converge to the true value?
So, what is the real distinction between these two concepts?
Can you give me an example of a process that is stationary without being ergodic, or ergodic without being stationary?

A random process is a collection of random variables, one for each time instant under consideration.  Typically this may be continuous time ($-\infty < t < \infty$) or discrete time (all integers $n$, or all time instants $nT$ where $T$ is the sample interval).

Stationarity refers to the distributions of the random variables. Specifically, in a stationary process, all the random variables have the same distribution function, and more generally, for every positive integer $n$ and $n$ time instants $t_1, t_2, \ldots, t_n$, the joint distribution of the $n$ random variables $X(t_1), X(t_2), \cdots, X(t_n)$ is the same as the joint distribution of $X(t_1+\tau), X(t_2+\tau), \cdots, X(t_n+\tau)$.  That is, if we shift all time instants by $\tau$, the statistical description of the process does not change at all: the process is stationary.
Ergodicity, on the other hand, doesn't look at statistical properties of the random variables but at the sample paths, i.e. what you observe physically. Referring back to the random variables, recall that random variables are mappings from a sample space to the real numbers; each outcome is mapped onto a real number, and different random variables will typically map any given outcome to different numbers.  So, imagine that some higher being as performed the experiment which has resulted in an outcome $\omega$ in the sample space, and this outcome has been mapped onto (typically different) real numbers by all the random variables in the process: specifically, the random variable $X(t)$ has mapped $\omega$ to a real number we shall denote as $x(t)$. The numbers $x(t)$, regarded as a waveform, are the sample path corresponding to $\omega$, and different outcomes will give us different sample paths. Ergodicity then deals with properties of the sample paths and how these properties relate to the properties of the random variables comprising the random process.

Now, for a sample path $x(t)$ from a stationary process, we can compute the time average
$$\bar{x} = \frac{1}{2T} \int_{-T}^T x(t) \,\mathrm dt$$  but, what does $\bar{x}$ have to do with $\mu = E[X(t)]$, the mean of the random process? (Note that it doesn't matter which value of $t$ we use; all the random variables have the same distribution and so have the same mean (if the mean exists)).  As the OP says, the average value  or DC component of a sample path converges to the mean value of the process if the sample path is observed long enough, provided the process is ergodic and stationary, etc.  That is, ergodicity is what enables us to connect the results of the two calculations and to assert that
$$\lim_{T\to \infty}\bar{x} = \lim_{T\to \infty}\frac{1}{2T} \int_{-T}^T x(t) \,\mathrm dt$$
equals
$$\mu = E[X(t)] = \int_{-\infty}^\infty uf_X(u) \,\mathrm du.$$  A process for which such equality holds is said to be mean-ergodic, and a process is mean-ergodic if its autocovariance function $C_X(\tau)$ has the property:
$$\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^T C_X(\tau) \mathrm d\tau = 0.$$

Thus, not all stationary processes need be mean-ergodic. But there are other forms of ergodicity too. For example, for an autocovariance-ergodic process, the autocovariance function of a finite segment (say for $t\in (-T, T)$ of the sample path $x(t)$ converges to the autocovariance function $C_X(\tau)$ of the process as $T\to \infty$. A blanket statement that a process is ergodic might mean any of the various forms or it might mean a specific form; one just can't tell,

As an example of the difference between the two concepts, suppose that $X(t) = Y$ for all $t$ under consideration.  Here $Y$ is a random variable. This is a stationary process: each $X(t)$ has the same distribution (namely, the distribution of $Y$), same mean
$E[X(t)] = E[Y]$, same variance etc.;
each $X(t_1)$ and $X(t_2)$ have the same
joint distribution (though it is degenerate) and so on.  But the process is not
ergodic because each sample path is a constant.  Specifically, if a trial
of the experiment
(as performed by you, or by a superior being) results in $Y$ having
value $\alpha$, then the sample path of the random process that corresponds
to this experimental outcome has value $\alpha$ for all $t$, and the DC
value of the sample path is $\alpha$, not $E[X(t)] = E[Y]$, no matter how
long you observe the (rather boring) sample path.  In a parallel universe,
the trial would result in $Y = \beta$ and the sample path in that
universe would have value $\beta$ for all $t$.
It is not easy to write mathematical specifications to
exclude such trivialities from the class of stationary processes, and so this is a very minimal example of a stationary random process that is not ergodic.

Can there be a random process that is not stationary but is ergodic? Well, NO, not if by ergodic we mean ergodic in every possible way one can think of: for example, if we measure the fraction of time during which a long segment of the sample path $x(t)$ has value at most $\alpha$, this is a good estimate of $P(X(t) \leq \alpha) = F_X(\alpha)$, the value of the (common) CDF $F_X$ of the $X(t)$'s at $\alpha$ if the process is assumed to be ergodic with respect to the distribution  functions.  But, we can have random processes that are not stationary but are nonetheless mean-ergodic and autocovariance-ergodic.  For example, consider the process
$\{X(t)\colon X(t)= \cos (t + \Theta), -\infty < t < \infty\}$
where $\Theta$ takes on four equally likely values $0, \pi/2, \pi$ and $3\pi/2$.
Note that each $X(t)$ is a discrete random variable that, in general, takes on four equally likely values $\cos(t), \cos(t+\pi/2)=-\sin(t), \cos(t+\pi) = -\cos(t)$ and $\cos(t+3\pi/2)=\sin(t)$,  It is easy to see that in general $X(t)$ and $X(s)$ have different distributions, and so the process is not even first-order stationary. On the other hand,
$$E[X(t)] = \frac 14\cos(t)+
\frac 14(-\sin(t)) + \frac 14(-\cos(t))+\frac 14 \sin(t)  = 0$$ for every $t$ while
\begin{align}
E[X(t)X(s)]&= \left.\left.\frac 14\right[\cos(t)\cos(s) +  (-\cos(t))(-\cos(s)) + \sin(t)\sin(s) + (-\sin(t))(-\sin(s))\right]\\
&= \left.\left.\frac 12\right[\cos(t)\cos(s) + \sin(t)\sin(s)\right]\\
&= \frac 12 \cos(t-s).
\end{align}
In short, the process has zero mean and its autocorrelation (and autocovariance) function depends only on the time difference $t-s$, and so the process is wide sense stationary. But it is not first-order stationary and so cannot be
stationary to higher orders either. Now, when the experiment is performed and the value of $\Theta$ is known, we get the sample function which clearly must be one of $\pm \cos(t)$ and $\pm \sin(t)$ which have DC value $0$ which equals $0$, and whose autocorrelation function is $\frac 12 \cos(\tau)$, same as $R_X(\tau)$, and so this process is mean-ergodic and autocorrelation-ergodic even though it is not stationary at all. In closing, I remark that the process is  not ergodic with respect to the distribution function, that is, it cannot be said to be ergodic in all respects.


Under what conditions do the phase margin and Nyquist criteria give the same results? When designing feedback systems, I often evaluate stability by thinking about phase margin: the closed loop system $$T(s) = \frac{L(s)}{1+L(s)}$$ is stable if $L(s)$ has positive phase margin, i.e.,
$$\angle L(j\omega) > -\pi/2$$ at all frequencies for which $$|L(j\omega)|>1$$. This is easy to visualize graphically using a Bode plot.

But for some systems -- typically involving right-half-plane poles and/or zeros in $L(s)$ -- this criterion doesn't work and we need to use the Nyquist criterion to check stability. That is, the number of unstable closed-loop poles is equal to the number of encirclements of $-1+j0$ by the polar plot of $L(s)$ when $s$ is evaluated along an appropriate contour. My understanding is that this Nyquist criterion is more general, and the "positive phase margin" rule is a simplification that only applies in certain cases.

I'm looking for a criterion that will let me decide when these two criteria (positive phase margin and the Nyquist test) will give me the same results. I suspect that it involves right-half-plane poles and/or zeros in $L(s)$, but I'm not sure.

As you have pointed out, the nyquist stability criterion is more general, moreover, is the only stability criterion. Nevertheless, the Bode plot give us the exact same information that the nyquist plot does, so the stability of a system can be determined with either the bode or the nyquist plot, some times its easier with the former and other times with the latter. for minimum phase system this is usually easier with the bode plot, since you only have to see that (as you said) $$\angle L(j\omega) > -\pi/2 $$
at all frequencies for which
$$|L(j\omega)| > 1 $$

Answering your question explicitly, this is only valid for minimum phase systems. For systems that are not minimun phase the gain and phase margins are not as straightfoward to find.


Speech dereverbaration via maximum kurtosis adaptive filtering I'm trying to code the algorithm described in Speech dereverbaration via maximum-kurtosis subband adaptive filtering by Gillespie, Malvar and Florencio, and the signal looks cleaner in when I plot it.
However there are 2 aspects of the results that are worrying me:

The sound level seems lower than the original when played.
And the sound comes out muffled.

I'm quite new to speech signal processing so I was wondering whether the latter issues occur commonly or is it just my bad coding?

Here is the part of the code that corresponds to the adaptive filter. I believe it should be the problematic part.

for n=1:500
 %Sig is the sum of the product of the FFT of the kusrtosis gradient and
 %the complex conjugate of the FFT of the LP residual of the reverberant
 %signal
    sig=0;
    ii=1;

    while ii<length(F)-L
        sig=sig+sum(F(ii:ii+L).*Yconj(ii:ii+L));
        ii=ii+L;
    end
 %H is the frequency domain representation of the filter and the following are the update equations

%Hpr is G' in the paper
    Hpr(n+1)=H(n)+(mu/M)*sig;

    if Hpr(n)==0 || isnan(Hpr(n))==1
       H(n+1)=0;
    else

        H(n+1)=Hpr(n)/abs(Hpr(n));
    end
 %getting the optimized signal
    Zt=Yleftres.*H;

    zt=ifft(Zt);
%updating the value of the kurosis gradient
    q=1:length(zt);
    while q<length(zt);

        secmoment=beta*secmoment+(1-beta)*zt(q:q+881).^2;
        fourthmoment=beta*fourthmoment+(1-beta)*zt(q:q+881).^4;

        f(q:q+881)=4*(secmoment*(zt(q:q+881).^3)-fourthmoment*zt(q:q+881))/(secmoment.^3);
        q=q+881;

        cleaner=isNaN(f(q:q+881));
        cleaner=cleaner-1;
        cleaner=abs(cleaner);
        f(q:q+881)=cleaner.*f(q:q+881);

    end
    F=fft(f);

   %   Hpr(n+1)=H(n)+(mu/M)*sig;
end

As was advised to me  by A_A in the comment section I will close this thread with a "self-answer".
The code present in the question is a code for a speech signal de-reverberation adaptive filter based on the Kurtosis signal of the signal's Lp residual. The original idea isn't mine (the reference is in the question). But the problem I faced was mainly due to the fact that my sensor did not convert to the desired result. I drastically increased the number of iterations the adaptive filter does and the result converged.(from n=500 to n=2000) There are others parameters that can be tweaked too (like mu) but I haven't had the time recently to thoroughly evaluate their effects.

Here's a link to download the code hope it helps. (expires by April 18th)

Code download link

Here's the link for the RIR code I used, I did not write the code though. I downloaded the file from matlab file exchange. The original author is Stephen G. McGovern, he also wrote a paper with the theoretical explanations for the RIR filter.

RIR code


Illustration of time domain aliasing cancellation (TDAC, MDCT, lapped transform) Can someone find or make a visual illustration of how time domain aliasing cancellation works in lapped transforms such as the MDCT?  I think a graphic would make it a lot easier to understand, and I can't find any.  For instance, this image clearly shows how the overlap-add process combines windowed chunks of a signal into the original signal:

But in MDCT, there is time-aliasing within each chunk?  And somehow overlapping and adding the chunks causes the aliasing to cancel out?

MDCT uses the DCT-IV, which inverts polarity outside of one boundary, so I assume it has something to do with that:

Edit: I have recently created two Jupyter Notebooks that illustrate this behaviour and let you play around with some actual matrices and actual signals.

I find understanding MDCT easiest if we define our transforms as matrix operations.

DFT

In case of a DFT such a matrix would look like

$$
F_M = \left( \sqrt{\frac{1}{M}} e^{-2 \pi i k n / M} \right)_{k,n=0 \dots M-1}
\in \mathbb{C}^{M \times M}
$$

$F_M \approx$

and our transform would become a matrix multiplication

$$
X = F_M x
\qquad{} x,X \in \mathbb{C}^M
$$

where $X$ and $x$ are column vectors of size $M$. The inverse of the matrix $F$ exists, so synthesis is possible

$$
x = F^{-1}_M X
\qquad{} x,X \in \mathbb{C}^M
$$

which means (note we are using the unitary matrix definition here)

$$
F_M F^{-1}_M = I
$$

Wide DFT

If we now just defined a "wide DFT" matrix to be larger in one dimension, lets say $M \times 2M$.

$$
G_{2M} = \left( \sqrt{\frac{1}{M}} e^{-2 \pi i k n / M} \right)_{k=0 \dots M-1 \quad{} n=0 \dots 2M-1}
\in \mathbb{C}^{M \times 2M}
$$

this is equivalent to having longer basis functions.

$G_{2M} \approx$

This means for $2M$ samples $x$ we get $M$ coefficients $\hat{X}$.

$$
\hat{X} = G_{2M} x
\qquad{} x \in \mathbb{C}^{2M}
\quad{} X \in \mathbb{C}^M
$$

Doing that, you might also notice is that there are now repititions in our matrix as our basis function $e^{-2 \pi i k n / M}$ starts to repeat after $M$ samples. It looks like two DFT matrices stacked next to each other, which we can represent as a folding matrix followed by DFT

$$
G_{2M} = F_M
\begin{bmatrix}
I_M & I_M
\end{bmatrix}
$$

$G_{2M} \approx$  $=$  $\times$

Where $I$ is the identity matrix. Obviously an inverse of $\begin{bmatrix}I_M & I_M\end{bmatrix}$ doesn't exist, so an inverse of $G$ doesn't exist either.

If we now look at the complete transform

$$
\hat{X} = F_M
\begin{bmatrix}
I_M & I_M
\end{bmatrix} x
\qquad{} x \in \mathbb{C}^{2M}
\quad{} X \in \mathbb{C}^M
$$

we see that always the sum of exactly two samples of $x$ are being transformed. This is called time domain aliasing.

We can do this for any size and the folding matrix will just repeat, e.g.

$$
G_{4M} = F_M
\begin{bmatrix}
I_M & I_M & I_M & I_M
\end{bmatrix}
$$

$G_{4M} \approx$  $=$  $\times$

and we will end up with the sum of more and more samples being transformed together, producing more and more aliasing.

DCT-IV

We can define DCT-IV as a matrix, too.

$$
C_M = \left( {\sqrt{\frac{2}{M}} \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{1}{2}\right)\right]} \right)_{k,n = 0 \dots M-1}
\in \mathbb{R}^{M \times M}
$$

$C_M \approx$

Interestingly, the matrix is orthogonal and its own inverse

$$
C_M = C^{-1}_M = C^T_M
$$

Wide DCT-IV

Again, we can define a wide DCT-IV of size $M \times 2M$ by just extending the basis functions.

$D_{2M} \approx$

However, this time the repititions in the basisfunction are different, and thus the folding matrix is:

$$
D_{2M} = C_M
\begin{bmatrix}
I_M & -J_M
\end{bmatrix}
\in \mathbb{R}^{M \times 2M}
$$

$D_{2M} \approx$  =  $\times$

where $J_M$ is the $M \times M$ antidiagonal matrix. If we repeat this further, we get

$$
D_{4M} = C_M
\begin{bmatrix}
I_M & -J_M & -I_M & J_M
\end{bmatrix}
\in \mathbb{R}^{M \times 4M}
$$

$D_{4M} \approx$  =  $\times$

after which our folding matrix will start to repeat. You will discover this is exactly the repetition pattern of the image on Wikipedia

We will need this pattern, so lets define it

$$
Q =
\begin{bmatrix}
I_M & -J_M & -I_M & J_M
\end{bmatrix}
\in \mathbb{R}^{M \times 4M}
$$

MDCT

Analysis

$$
X_k = {\sqrt{\frac{2}{M}} \sum^{2M-1}_{n=0} h(n) x_n \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{M}{2}+\frac{1}{2}\right)\right]} \qquad{} k = 0 \dots M-1
$$

If we look at the definition of the MDCT we note that we are transforming $n = 0 \dots 2M$ samples to $k = 0 \dots M$ coefficients, so we are using a "wide matrix" of size $M \times 2M$.

$$
E_{2M} = \left( {\sqrt{\frac{2}{M}} \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{M}{2}+\frac{1}{2}\right)\right]} \right)_{k = 0 \dots M-1 \quad{} n=0 \dots 2M-1}
\in \mathbb{R}^{M \times 2M}
$$

$$
X_k = E_{2M} x
\qquad{} x \in \mathbb{R}^{2M}
\quad{} X \in \mathbb{R}^M
$$

But also, if we compare the basis function to DCT-IV, we notice they are almost the same, but with an offset of $M/2$ for $n$. This means we are using a "shifted view" of DCT-IV matrix, shifted to the right:

$$
E_{2M} = \left( {\sqrt{\frac{2}{M}} \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{1}{2}\right)\right]} \right)_{k = 0 \dots M-1 \quad{} n=\frac{M}{2} \dots \frac{M}{2}+2M-1}
\in \mathbb{R}^{M \times 2M}
$$

We can express such a shifted view as a cyclic permutation matrix $P\in \mathbb{R}^{4M \times 4M}$, but only if our basis function was cyclic, too. We know $D_{4M}$ is expanded to the point where it becomes cyclic, so we can express $E_{2M}$ using $D_{4M}$ and using an extraction matrix and a permutation matrix:

$$
E_{2M} =
D_{4M}
P
\begin{bmatrix}
I_{2M} \\
0
\end{bmatrix}
$$

$E_{2M} \approx$  $=$  $\times$  $\times$

but since we know how to represent wide DCT-IV by a regular DCT-IV like

$$
D_{4M} = C_M Q
$$

we can write

$$
E_{2M} = C_M Q P
\begin{bmatrix}
I_{2M} \\
0
\end{bmatrix}
$$

$E_{2M} \approx$  $=$  $\times$  $\times$  $\times$

We don't need all of those shift matrices individually, so let's multiply them together into one matrix $R$:

$$
R = Q P
\begin{bmatrix}
I_{2M} \\
0
\end{bmatrix}
$$

$E_{2M} \approx$  $=$  $\times$

And end up with

$$
X = C_M R x
$$

This means MDCT is a permutation and superposition and then a DCT-IV of our input samples.

Synthesis

For reconstruction we need to form a linear combination of basis vectors, which can be written as

$$
\hat{x} = (C_M R)^T X = R^T C_M^T X
$$

And substituing our forward transform for $X$ we get

$$
\hat{x} = R^T C_M^T C_M R x
$$

we know

$$
C_M^TC_M = I_M
$$

but

$$
R^T R \neq I_M
$$

instead

$R^T R \approx$

But if we now overlap+add two frames with just the right hopsize

we get

perfect reconstruction in the overlapping part.

This means, if we first transformed our signal in windows with half a window overlap, all the aliasing produced in this overlap will exactly cancel out, making the entire transform perfectly reconstructing.

Caveats/Further reading

Windowing

All of this is obviously disregarding windowing. But if you look at the permutation matrix and look at what sample is superpositioned with what other sample it should be at least somewhat intuitively visible where the Princen Bradley condition comes from.

Nicer Permutation

Some authors prefer the MDCT to be defined as

$$
X_k = {\sqrt{\frac{2}{M}} \sum^{2M-1}_{n=0} h(n) x_n \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{M}{2}-\frac{1}{2}\right)\right]} \qquad{} k = 0 \dots M-1
$$

(note the minus). If you then solve all of the remaining math, you will find a permutation matrix that looks like

$\hat{R} \approx$

as opposed to

$R \approx$ .

This doesn't change any of the remaining explanations (TDAC still works etc.), but you can see that the samples in the center of your frame are now 'nicer': they're positive valued and transformed forwards, as opposed to negative valued and transformed backwards.

Math for overlap add

If you feel that last bit became a bit too "hand wavey", that's because it is. The traditional math of transforms or matrices starts to break down if you start talking about overlapping transforms. A very neat way of dealing with them are Polyphase Matrices in the z-Domain. They're hard to wrap your head around at first though.


16 bit PCM Audio Sounds Like Bit Crushing I'm developing an Android app which uses USB soundcard for input and output. The audio processes are in the native (C++) side.

The device's bit resolution is 16. I get the audio stream from the device as unsigned char* and have to convert them to float between -1.0f and +1.0f to process.

I do the conversion but it sounds like bit crushing effect if I increase the volume of input.

Here is what I do when converting:

For Input Conversion:

/*** buffer is an array of float ***/
/*** temp_buffer_in is the buffer coming from the device ***/

float divider = (float) std::pow(2, (BIT_RESOLUTION_IN - 1));
int BYTE_RESOLUTION_IN = BIT_RESOLUTION_IN / (sizeof(char) * 8);
int bytes = buffer_size * BYTE_RESOLUTION_IN; //Buffer size in bytes

for (i = 0, j = 0; i < bytes; i += BYTE_RESOLUTION_IN, j++) {

    UnionInteger16 union_int; // Union structure below
    union_int.i = 0;

    for (int k = 0; k < BYTE_RESOLUTION_IN; k++) {
            union_int.c[k] = temp_buffer_in[i + k];
    }

    buffer[j] = (float) union_int.i;

    buffer[j] = buffer[j] / divider;
}

For Output Conversion:

/*** temp_buffer_out is the buffer sent to the device ***/

float multiplier = (float) std::pow(2, (BIT_RESOLUTION_OUT - 1));
int BYTE_RESOLUTION_OUT = BIT_RESOLUTION_OUT / (sizeof(char) * 8);
int bytes = buffer_size * BYTE_RESOLUTION_OUT;

for (i = 0, j = 0; i < bytes; i += BYTE_RESOLUTION_OUT, j++) {

    if (buffer[j] >= 0) {
        buffer[j] = buffer[j] * (multiplier - 1.0f);
    } else {
        buffer[j] = buffer[j] * multiplier;
    }

    UnionInteger16 union_int; // Union structure below
    union_int.i = (int16_t) buffer[j];

    for (int k = 0; k < BYTE_RESOLUTION_OUT; k++) {
        temp_buffer_out[i + k] = union_int.c[k];
    }
}

I use Union for converting unsigned chars to int16_t and then to float. Here is the union structure:

union UnionInteger16 {
    int16_t i;
    uint8_t c[2];
};

Everything looks fine in the code but as I mentioned, the audio sounds like bit crushing when the input volume is increased and when the volume is decreased, the audio starts to disappear. I think I lost the dynamic range of the audio.

Sorry for the long question. I would be grateful if you help me. Thanks in advance.

In case anyone has the same problem, I want to tell that the above code works really well. My problem about the distortion of audio caused by wrong buffer size.

But while searching for a solution, I figured out that everyone says that using union may cause undefined behaviour in different CPU structures. It is said to use bitshift operations instead. So I decided to change the code for conversion as below:

Input Conversion

int16_t sample = 0;
    for (int k = 0; k < BYTE_RESOLUTION_IN; k++) {
    sample |= input_buffer[i + k] << (k * 8);
}

Output Conversion

int16_t sample_int16 = (int16_t) temp;
    for (int k = 0; k < BYTE_RESOLUTION_OUT; k++) {
    output_buffer[i + k] = (sample_int16 >> (k * 8)) & 0xFF;
}

And this conversion is made for Little-Endian system. If you have Big-Endian system, you have to change the k value properly.

Have a nice day


Plotting Power Spectrum (Matlab) I want to make a plot of the power spectrum for a .wav sound file, over the frequency range from -2000 to 2000 Hz.

Attempt:

I used the following code, but the resulting plot is one sided. It doesn't plot the negative frequencies. But the graph I need to produce has to contain negative frequencies. What do I need to do?

Any help would be appreciated.

Here is my code:

[s, Fs] = wavread('chord.wav');
nfft = 2^nextpow2(length(s));
Pxx = abs(fft(s,nfft)).^2/length(s)/Fs;
Hpsd = dspdata.psd(Pxx(1:length(Pxx)/2),'Fs',Fs);
plot(Hpsd)
axis([-2 2 -150 -35])  %change axis range

Produces:

First of all, wavread has been removed in the 2015 & 2016 versions of MATLAB. So, if you plan on using your code in a newer version of MATLAB, I recommend that you use the audioread function instead.

Now, to answer your question. MATLAB already have a built-in function for plotting the power spectral density (Pxx) of a given input signal - periodogram. I've taken the liberty to copy a piece of sample code from MATLAB's documentation that should help you plot the Pxx you want (it will contain negative frequencies as well).

periodogram(x,[],NFFT,Fs,'centered')

All you need to keep from your code snippet are the input values seen here. The rest, like the axis function, aren't necessary to get the plot you want if you use this periodogram function.
The above function will plot the power spectral density of a given input signal x with no specified window ([] which just means the function defaults to a simple rectangular window). The final argument is the one that I presume you are most interested in: freqrange. It can have several different inputs, including twosided and centered. You should try both options and see which output you find more useful.

I hope this helps you!


I Q sampling and baseband version of analytic signal Is it correct to say that if we have a radio signal $s(t)$ centered around the angular frequency of $\omega$ as $\omega\pm\omega_B/2$ (where $\omega_B$ is the bandwidth of the signal) and the corresponding analytic signal is $s_a(t)=s(t)+j \hat{s}(t)$ (where $\hat{s}(t)$ is the Hilbert transform of $s(t)$) then $I$ and $Q$ as received from an SDR receiver (RTL-SDR) can be represented as

$I=Re(s_a(t) e^{-j \omega t})$
$Q=Im(s_a(t) e^{-j \omega t})$

In other words is it correct that $I$ and $Q$ are the real and imaginary parts of the analytic signal shifted to baseband?

I'm trying to derive the above from the assumption that the $I$ component is a low-pass filtered version of $s(t) \cos\omega t$ and the $Q$ component is a low-pass filtered version of $s(t) \sin\omega t$.

Considering that

$s(t) \cos\omega t = (s_a(t)-j\hat{s}(t))\times {1\over{2}}(e^{-j\omega t}+e^{j\omega t})$
$s(t) \sin\omega t = (s_a(t)-j\hat{s}(t))\times {1\over{2}}j(e^{-j\omega t}-e^{j\omega t})$

and then expanding these and discarding high frequency components (both negative and positive) I arrived at the following

$I={1\over{2}}(s_a(t) e^{-j\omega t} - j\hat{s}(t)e^{j\omega t})$
$Q={1\over{2}}j(s_a(t)e^{-j\omega t} + j\hat{s}(t)e^{j\omega t})$

And I expected that $I+jQ$ should give me the shifted analytic signal $s_a(t)e^{-j\omega t}$ but when I substituted $I$ and $Q$ from above into $I+jQ$ I ended up with $-j\hat{s}(t)e^{j\omega t}$ which is the negative frequency replica of the signal shifted to baseband, which would be a mirrored version of $s_a(t) e^{-j \omega t}$ (i.e. of what I expected).

Is my initial statement about the relation between $I$ and $Q$ and the analytic signal correct and there's just a mistake in my math or am I making a wrong assumption somewhere along the way?

The problem in your derivation is the way you discard the high frequency components. Note that $\hat{s}(t)$, the Hilbert transform of the bandpass signal $s(t)$, has frequency components around $\omega_c$ as well as around $-\omega_c$. So you can't just discard one of the two terms $\hat{s}(t)e^{j\omega_ct}$ and $\hat{s}(t)e^{-j\omega_ct}$. The low pass filter filters out parts of both terms, but leaves other parts unchanged.

You could better represent $s(t)$ as

$$s(t)=\frac12\big[s_a(t)+s_a^*(t)\big]\tag{1}$$

where $s_a^*(t)$ is the complex conjugate of the analytic signal $s_a(t)$. This has the advantage that both components on the right-hand side of $(1)$ have only frequency components either at $\omega_c$ ($s_a(t)$) or at $-\omega_c$ ($s_a^*(t)$), but not at both frequencies.

So, for example, the in-phase component is easily derived as a low pass filtered version of $2s(t)\cos(\omega_ct)$:

$$\begin{align}I(t)&=\textrm{LP}\big\{2s(t)\cos(\omega_ct)\big\}\\&=\frac12\textrm{LP}\left\{\big[s_a(t)+s_a^*(t)\big]\left[e^{j\omega_ct}+e^{-j\omega_ct}\right]\right\}\\&=\frac12\left[s_a(t)e^{-j\omega_ct}+s_a^*(t)e^{j\omega_ct}\right]\\&=\textrm{Re}\left\{s_a(t)e^{-j\omega_ct}\right\}\tag{1}\end{align}$$

where $\textrm{LP}\{\cdot\}$ is the low pass filter operator.

In a similar way, the quadrature component $Q(t)$ can be shown to be a low pass filtered version of $-2s(t)\sin(\omega_ct)$.


Output power of a signal passing through an IIR filter suppose I have the following discrete time signal
$x\left(n\right)=0.7\sin\left(2\pi \times 623\:nT\right)$
where $N=4000$ and $T=\frac{1}{f_s}$, where $f_s=40$ kHz
I pass this signal through a simple IIR filter $ H\left(z\right)\:=\:\frac{0.2}{1-0.3z^{-1}}$
I calculate the power of the ouput signal $P_y\:=\:P_x\:\left|H\left(z\right)\right|^2_2$
basically am using the filter 2-norm in this case.
for that my calculations would be as follows:
$P_y\:=\left(0.245\right)\left(\frac{0.2^2}{1-0.3^2}\right)\:=\:0.0108$
However when I develop this in Matlab, I get the output power as 0.0199
Here's my code in Matlab
fs = 40000;
f1 = 623;
N = 4000;
n = 1:N;
b = 0.2;
a = [1 -0.3];
x = 0.7*sin(2*pi*f1*n/fs);
P_x = var(x);
y = filter(b,a,x);
P_y = var(y);

I don't understand what I was doing wrong, would someone please clarify for me? I would appreciate it. Thank you in advance.

The analog frequency $\Omega=2\pi\times 623$ and the digital frequency is $\omega = \Omega T = \Omega/f_s = 0.03115\pi$
The frequency response at this frequency is
$$H(z)|_{z=e^{j\omega}} = H(e^{j0.03115\pi}) = \frac{0.2}{1-0.3 e^{-j0.03115\pi}}$$
and its square of magnitude is
$$
|H(e^{j0.03115\pi})|^2 = \Big|\frac{0.2}{1-0.3 e^{-j0.03115\pi}}\Big|^2 = 0.0812
$$
Now you'll see $P_y = P_x |H(e^{j0.03115\pi})|^2 = 0.0199$


comparing 2 audio files in c# I have 2 audio .wav files recorded by the the same person with the same word. Each audio file has a format of 16bit, 44.1kHz and is a mono channel. I would like to know how to compare these 2 audio files if it is by the same person. What I know is to process the files to produce a floating point result of values between -1 and 1 and apply fft on it and compare the frequency. What I would like to know is how to compare the audio files, prefably with the code which shows the comparison and if possible the steps before that(processing to become values btw 1 and -1).Thanks. Im relatively new in this area so I would appreciate If u guys would help me out. I'm currently using a metro application in vs. I would prefer not to use any library like naudio or bass.net but I would still accept any answer.
A step by step solution starting from utilizing the .wav file is much appreciated. I have been trying to search topics related to what I'm doing but what I have found are small parts of it.

I'm thinking the signals could be compared by

$D = \mathcal{F}^{-1}(\mathcal{F}(A) - \mathcal{F}(B))$

Where A and B are the two audio files to be compared.

I guess after that you could come up with some criteria for deciding that the voices are different. Perhaps when $|D|$ goes over some threshold.

You could also use Cross correlation

Using some other library to help with this would save time. But I haven't done this with audio files.

AForge Looks like it has FFT although I haven't used it.


WiFi Channelization for 40, 80 and 160 MHz I have some questions on WiFi channelization in 802.11n and 802.11ac.  Here is the snip of List of WiFi channels (source: WiKi):

I am working on the L1 layer for 802.11n and 802.11ac model where I am transmitting and receiving some WiFi packets
Questions

Transmitting and Receiving WiFi Packets on 20MHz channels is pretty straight forward. If I want to use 40 MHz channel (Channel 38) , I need to do channel bonding, where channel 36 and Channel 40 are used.
Is the centre frequency going to be 5190 (channel 38) ? If that's the case, why we even bother to use the term primary and secondary channels?
What about the reception side? Do I need to tune my Rx to Channel 38 with samplerate > 2·40 MHz, or do I need to use any of the primary and secondary channels?

Your receiver does not know what BW you are using. Hence it will always be tuned to receive 20MHz BW (ch 36 or 40). So for the preamble part you need to transmit in channel bonding (2 x 20MHz) mode. Once the receiver has detected preamble it decodes the header to find out that data part is 40MHz. You need to turn on a bit in header for this (forgot the exact bit number) you can find it in standard. Same case for 80MHz BW. Preamble part is 4 x 20MHz.
Yes you need to continuously tune to 20MHz channels in the receiver. Otherwise no way of knowing what BW transmitter is sending. For fast synchronisation this is the feasible method. Once you detect preamble you need check for that header bit in both primary and secondary channels to check for 40MHz BW.


Understanding equal power crossfades I'm trying to implement an equal power crossfade between two uncorrelated signals ($A$ and $B$). From what I understand I need to setup the crossfade in the following way

Come up with a function $f(x)$ that determines the gain applied to a signal where $x$ is from $0$ to $1$. And $0$ represents being all the way faded to signal $A$,and $1$ represents all the way faded to signal $B$.

$f(x)$ should satisfy:
\begin{cases}\begin{align}
f(x)^2 + f(1-x)^2 &= 1\\
f(1) &= 1\\
f(0) &= 0
\end{align}\end{cases}

So I could use the square root cross-fade $f(x) = \sqrt{x}$:

Or I could use the $\cos$ cross-fade $f(x) = \cos\left(\frac{π}{2}x\right)$:

My question is where does this $f(x)^2 + f(1-x)^2 = 1$ restriction come from? I know that we interpret sound pressure on a logarithmic scale but I don't understand why a quadratic relationship is required between two signals to maintain equal perceived volume.

Power of a (real-valued) digital signal $s[n]$ is simply the $s[n]^2$ (or proportional to that, depending if you want to normalize or not).

That's why the constant sum power constraint is $f^2(x) + f^2(1-x) = 1$. That's simply the definition of constant sum of  powers, nothing more.

so, what you've not specified is "logarithmic scale" to what, or what the perceived quantity is.

And the answer is perceived volume is logarithmic to power. The fact that this relation is logarithmic doesn't really matter here – it's only important that to achieve the same perceived volume, you need to have the same power.

Your functions $f$ work on amplitude, and thus, the square comes from the amplitude-to-power relation, not from the power-to-perceived-volume relation.

Ok, this question has become kind of a moving target. I'll try to explain what's up here, in the order my process of thought allows me to explain it:

Perceived volume is a function of signal energy. That is, to get the same perceived volume, you need to produce the same energy.
Power of a signal is proportional to the square of its amplitude, so with $s(t)$ being our signal, $s^2(t) =  \text{const.}$ restraint is very intuitive; energy of that  signal within $\tau$ time units is $\int_{T_0}^{T_0+\tau}s^2(t)\,dt$.
Your signal processing happens on amplitudes; you add amplitudes, not powers. Thus, let $s(t) = a(t) + b(t)$. It follows directly that $s^2(t) = a^2(t) + 2a(t)b(t) + b^2(t)$. Notice the middle $ab$ term!
But the formula we have dictates that $a^2(t) + b^2(t) = 1$. This only makes sense if $\int_{T_0}^{T_0+\tau} 2a(t)b(t)\,dt\equiv 0\quad\forall T_0, \tau$. Which is exactly the case if $a$ and $b$ are uncorrelated. If $a$ and $b$ were correlated, the $ab$ crosscorrelation term could not be equal to 0.


How to use Expectation maximization to estimate poisson noise matlab I delved into the state of the art of algorithms for Poisson noise estimation in order to estimate the variance, I found that the Expectation Maximization algorithm is very used and it is very effective and easy to use,

But I have not found the matlab code, I found only one that is was used for classification,

can you give me the code or the algorithm for estimating and not for classification please.

I would be grateful if you answer me

Explanations of the Expectation Maximisation algorithm:

Nature.com ( Nature Biotechnology 26, 897 - 899 (2008) )

Wikipedia - Expectation–maximization algorithm

Stack Overflow - What is an intuitive explanation of Expectation Maximization technique?


How to implement a filter associated to a specific wavelet I am working on ECG signals, to eventually extract features in order to detect an arrhythmia and classify it.
I am using Discrete Wavelet Transform with biorthogonal wavelet bior6.8
During my research, I came to know that wavelet transform is the convolution of the input processed signal with the daughter wavelets to get approximation and detail

$$X(a,b)=\frac{1}{\sqrt{a}}\int_{-\infty}^{\infty} \Psi\left( \frac{t-b}{a}\right)x(t)dt $$
where $a$ is scaling and $b$ is time.

I can't find the expression of the mother wavelet anywhere also when it came to practice DWT is usually presented as a filter bank of high pass and low pass filter
The question is what is the difference between the different wavelets if it is always presented as a bank of filters

In my work, I used these two Butterworth high pass and low pass filters, but I still can't explain my choice, I read that Butterworth is the most used in signal processing and that it optimizes the frequency response in the passband, getting as much as you can from the wanted frequency
Still, I have no arguments why I shouldn't use any others and not sure whether it is the correct way to implement bior 6.8 wavelet as I do not know any other way to implement wavelets and I would implement Daubechies or any other the same which do not make sense

from scipy.signal import filtfilt , butter

def butter_highpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='high', analog=False)
    return b, a

def butter_highpass_filter(data, cutoff, fs, order=5):
    b, a = butter_highpass(cutoff, fs, order=order)
    y = filtfilt(b, a, data,padlen=0)
    return y

from scipy.signal import butter, lfilter

def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y

So basically and to resume, the question is what is the difference between wavelet transforms and how can we value that during implementation and then, the use of BF in this type of wavelet is it correct?

NB: the image inserted is from Wikipedia


Image processing: how do I erode on boundary and dilate inside boundary? I have an image that I am trying to process:

I want two seamingly contradictory effects:

I want to remove "holes" inside "the blob". That is I want all pixels inside the blob to be solid red.
I want to move the boundary of the blob inwards.

From what I gather the morphological operation dilatation would do 1. and erosion would do 2?

But dilatation and erosion counteract each other?

How would you suggest I go about to achieve the desired result?

You will want to use the morphological closing operation, which is the erosion of the dilation of an image. This means that you first dilate the image, removing the small 'holes' in the interior of the image, and then you erode the image, shrinking the boundary. If you want the final boundary to be 'smaller' than the original boundary, simply erode more than you dilate, ensuring that you dilate enough to completely fill the interior holes first.

The idea behind this is that by removing the holes via dilation, the information about those holes is lost from the image. When you then proceed to erode the image, there are no longer interior gaps to erode away from, so it only erodes away from the outer image. Therefore, despite being opposing operations, the erosion does not simply reverse the dilation, provided the feature sizes to be removed are 'smaller' than the size of the dilation.


How to Use the DFT (FFT) to Solve a Least Squares Regularization Problem (Inverse Problem)? Let $X$ and $K$ be an image and a Point Spread Function (PSF), respectively.

The blur image $B$ is obtained as follows

$$B = X * K$$

I want to solve the following general regularization problem

$$\min_X \left\|X * K - B\right\|_2^2 + \lambda \| f(X) \|_2^2$$

where $f$ is a regularization function. In some literature I have seen, the authors use the FFT to solve such a problem. However, I cannot find any resources that show the procedure. My questions are:

How one can use FFT to solve the above problem?
Is there any condition that must be satisfied to use FFT?

The question really depends on $ f \left( \cdot \right) $.
Yet in order to show how to use FFT we can even use 1D signals.

Let's rewrite the problem:

$$ \hat{x} = \arg \min_{x} \frac{1}{2} \left\| K x - b \right\|_{2}^{2} + \frac{\lambda}{2} \left\| f \left( x \right) \right\|_{2}^{2} $$

The derivative is given by:

$$ g = {K}^{T} \left( K x - b \right) + \lambda f' \left( x \right)^{T} f \left( x \right) $$

Now, $ K $ as a matrix is given by a Circulant Convolution Matrix.
Hence the operation $ K v $ or $ {K}^{T} v $ can be done in the Fourier Domain.
Under some circumstances it might accelerate the operation significantly.


SDR spectrum anomalies Is it normal to see periodic peaks going from 48.75 - 90 MHz? Also i was checking on my SDR dongle once and the whole band was raised above background offset for like a whole 50 MHz bandwidth and then it dropped to background again? Does anybody know what that means?

Periodic noise in an SDR spectrum is a common problem.  One likely source is EMI or RFI from nearby electronics.

The sources of EMI and RFI can be from computers, USB cables, networks cables and routers, power cords, wall plug or internal power supplies for various electronic devices, wireless chargers, LED lighting controls and voltage converters, and solar panels, etc., even radiated from neighboring buildings.  This noise might rise and fall as people or controllers in other locations turn on/off lights and other devices. The power and high edge rates generated by some of these electronic components can radiate a vast range of periodic harmonics into the RF spectrum.

Kill the circuit breaker to your house or building, or wait for a neighborhood wide power outage, and a vast amount of this periodic noise in the spectrum might disappear.  Or, if you are lucky, you can turn off or unplug things nearby, one at a time, and discover the dominant culprit.  If you are running the SDR from a laptop and have a portable antenna, try running the SDR well away from any buildings and utility lines (for instance: a nearby park or open field, etc.)


Possible to shift image pixel-wise with FFT? I know we can shift an image entirely with FFT as

https://stackoverflow.com/questions/25827916/matlab-shifting-an-image-using-fft

However, is it possible to shift the image pixel-wise? I just found a similar description in a recently published paper:
X. Wang, L. Li, and G. Hou, "High-resolution light field reconstruction using a hybrid imaging system," Appl. Opt.  55, 2580-2593 (2016).
They use a subband approach that really confused me.

The point that with "pixel-wise" you mean "differently for each pixel" indicates that a frequency-domain transform is a really bad approach for this. The DFT of a single element is that element (scaled, if at all modified), so you can't win with this.

You can use the FFT and similar algorithms to do elegant shifts, but the fact that it's used in light-field reconstruction stems from the fact that there's a certain, mathematically elegant model of what should happen to many of those pixels that relates nicely to frequency domain.

If you just want to shift/interpolate/distort your image pixel-wise, the normal distortion matrices are the way to go – ask google about that.

For example, OpenCV has a very short tutorial on camera calibration; maybe the theory contained in there is something for you?


What does clipping gaussian noise mean? I found the subject at here like this "Does clipping noise colour the spectrum"
When I researching for clipping noise.
But I can't understand because typically something has clipped means something has removed.
But he has some different point of view about this.
so I want to know what is the clipping noise and how to undertand the posing about "Does clipping noise colour the spectrum"?
update
Especially, I want to know the meanning  about clipping gaussian noise.
Especially this part
IV. CLIPPING (CENSORING)
A. Clipped observations model

If your noise has independent and identically distributed samples from a zero-mean distribution (for example Gaussian), it is white. Other definitions of white noise also require the distribution to be symmetrical, but that is not required for the spectrum to be flat. Clipping the samples of such white noise will only change the common probability distribution of the samples:

This will affect the mean of the distribution, but if both sides are clipped in a way that preserves the zero mean of the distribution, then there is no "coloring" of the spectrum by a 0 Hz peak and the noise remains white.

Example in Octave (MATLAB clone):

Create random variables from a Gaussian distribution:

N = 65536;
x = normrnd(zeros(N, 1), ones(N, 1));

Clip from above:

y = min(x, 1);

Clip symmetrically from below:

z = max(y, -1);

Plot values of Gaussian random variables (blue), clipped above (green), clipped symmetrically (red). Horizontal axis = index, vertical axis = value:

plot(1:N, x, 1:N, y, 1:N, z)

Plot spectrum of Gaussian random variables (blue), clipped above (green), clipped symmetrically (red). Horizontal axis = frequency, vertical axis = magnitude:

loglog(1:N/2, abs(fft(x)(1:N/2)), 1:N/2, abs(fft(y)(1:N/2)), 1:N/2, abs(fft(z)(1:N/2)));

The sequence that was clipped from above has a spectral peak at 0 Hz. Otherwise, the spectral envelopes are flat. This is not very easy to see because of the increasing density of plot points in the high frequencies in the logarithmic scale plot.


How to reconstruct an image from s-parameters in Matlab? I have a set of data which is the s-parameter (s11), consisting of 2 columns ($1001 \times 2$). The first column represents the frequency and the second column is the magnitude in dB. I need to use these data to reconstruct an image using Matlab but I have no idea how to do it.

I know that I have to use ifft function first to transform from frequency domain to time domain but the result I had consists of real and imaginary part and I do not know how to deal with this.

Each complex value can be described as $$A \cdot e^{i\cdot \phi}  $$.
Where A is the magnitude and $\phi$ the angle.
With this knowlege you can use the ifft as you allready mentioned.
If you consider that the image is a real signal you can take the real part by using real().

So overall it the part after reading the magnitude and angle would look like

IMG = uint8(real(ifft2(mag .* exp(i * angle))));

Where IMG is the resulting image, mag the array of magnitudes (in linear scaling), angle the array of angles.

The only problem with that is, that you have to know the dimensions of the image. The spectrum of an image has (after FFT) the same dimensions as the picture (113 x 200 pixel picture -> 113 x 200 matrix after FFT).
So your mag and angle must have the dimension of the resulting picture (e.g. 113 x 200).


Fourier convolution of a histogram Notation:

$\mathcal F\left\{a\right\}$ denotes applying the discrete forward Fourier transform to the histogram $a$. Similarly for $\mathcal F^{-1}\left\{a\right\}$ and the discrete inverse Fourier transform.
Capitals represent the discrete fast Fourier transform of the lower-case counterpart, e.g. $A = \mathcal F\left\{a\right\}$.
The convolution of two histograms $f$ and $g$ is denoted $f\star g$.

Say I have histograms $a$ and $b$ with $a=[1,2,3]$ and $b=[4,5,6]$ (the heights) and the edges given by $e=[0,1,2,3]$. I.e. For the histogram $a$, there is one entry in the interval $[0,1)$, two entries in the interval $[1,2)$, and three entries in the interval $[2,3]$.

In order to perform a linear convolution as opposed to circular, $a$ must be padded such that len($a$) = len($a$) $+$ len($b$) $-1$. Similarly for $b$.

The convolution theorem states:
$$\mathcal F\left\{a\star b\right\} = \mathcal F\left\{a\right\} \cdot \mathcal F\left\{b\right\}$$

So if $c = a\star b$ and $d = \mathcal F^{-1}\left\{A\cdot B\right\}$, then $c = d$.

On performing the Fourier convolution of the $a$ and $b$ mentioned above,  I obtain a result if $d = [4,13,28,27,18]$.

My questions:

Where does $c$ sit in terms of bins?
Do I need to expand the bin edges, $e$, in some way?
In which direction should I expand them?
Should the extra bins have the same width as the originals (in this case, $1$)?

Intuition tells me that there should be an extra bin on either side of $e$, of the same width.

Furthermore, if I was to crop $c$ such that it has the same length as the unpadded kernel or the signal, where would I make the cut? Is this even possible?

Actually, it does not really matter, if a or b are histograms or time domain sequences. However, I suppose you want to convolve to histograms to get the probability distribution of a random variable that is the sum of two variables. Then, its PDF is the convolution of the PDF of the original variables. But, considering the convolution of two time-domain sequences might make it easier to understand.
The convolution theorem states fft(ab)=fft(a)*fft(b), where * denotes the N-point circular convolution and a and b are of length N points. Your statement of the convolution theorem fft(ab) = fft(a) . fft(b) is not correct.
As I said, the convolution theorem in discrete domain corresponds to circular convolution. Most probably, you want to have a linear convolution. Then, you need to zero-pad a and b before applying the convolution theorem.
If the original bins range from 0-3 (left-side of the interval), then the bin edges after convolution would range from 0-6, because the (linear) convolution has a length of (N+N-1) samples, where N is the length of the sequence (which is 4 in your case). You can also understand this from the fact that the convolution of the histogram is the PDF of the sum of both random variables: The minimum sum would be 0, the maximum sum would be 3+3=6.

For more information on circular convolution, you might have a look at an article I wrote or the Wikipedia entry.


Discrete wavelet transform disadvantages I read in a paper that the  discrete wavelet transform (DWT) has two disadvantages
The first one is the shift variance property due to the downsampling process. Could you please help me understanding why downsampling leads to shift variance?
The second disadvantage is the lack of directional selectivity.  Why does DWT lacks directional selectivity?

A $2$-channel stage of a wavelet transform, combines two filters in parallel, followed by a down-sampling by two. The later is the cause for shift-invariance, as the filters are time-invariant. Signals $$x_0[n] = \{\ldots,0,1,0,1,0,1,\ldots\}$$ and  $$x_1[n] = \{\ldots,1,0,1,0,1,0,\ldots\}$$ which are shifted by only one sample, yield respectively  $$y_0[n] = \{\ldots,0,0,0,0,0,0,\ldots\}$$ and $$y_1[n] = \{\ldots,1,1,1,1,1,1,\ldots\}$$

Answers to down-sampling a thus not shift-invariant. This happens even with low-pass or high-pass filters. In an $L$-level wavelet decomposition, the overall wavelet filter bank is invariant to multiples of $2^L$ shift, not the intermediate one.

When applied in 2D, classical wavelet schemes, by simplicity, apply 1 DWT on rows and columns of the image, in a separate way. The resulting 2D wavelet filter is of rank one, very poor at separating directions other than vertical or horizontal. However, genuine 2D wavelets with better  directional selectivity exist. The paper A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity, Signal Processing, 2011, is devoted to that topic.


Why is a first/second derivative useful in spectroscopy? I'm currently working on Raman spectroscopy and while reading some literature I came across the first and second derivative of a spectrum. It's not clear to me why they are useful to look at since they look very noisy to me.

I tried graphing the first and second derivative of a spectrum I attained from liquid whey (see the following pictures, 1. accquired raman spectrum (smoothed), 2. first derivative, 3. second derivative):

What information can I read from the first/second derivative in this specific example and in general?

Derivative spectroscopy was at its peak in the 70-80s. It is still very useful. The only problem with derivatives is that they enhance noise like crazy i.e., they suppress low frequency signals and enhance high frequency signals.

Before you do anything with the derivatives, you have to apply a smoothing process on them, specially on the second derivative. This "problem" of noise with derivatives was very well known then and smoothing was necessary before anyone used them.

For smoothing, you have apply a filter in such a way that it does not shift the peak positions, and the simplest one is the centered moving average. Choose an odd number of data points, say 5 (or higher), so the data at point x_3 will have two data points above and two data points below. Take the average of all 5, at position number x_3. Continue doing this for the rest of the points. You will start to see features in the second derivative.

Basically, derivatives are used to reveal the fine structure of a band. Shoulders become clear. The second derivative can "act" like a high resolution signature of a given molecule.

A very nice paper on higher order derivatives is
High‐Resolution, Higher‐Order UV/VIS Derivative Spectrophotometry

And there are some other tricks with second derivatives to resolve overlapping signals by my co-author.

https://terpconnect.umd.edu/~toh/spectrum/Differentiation.html


Comparing a clean, noisy and enhanced signal I was doing some experiments on analyzing few audio samples and I'm stuck with this. Suppose we have a clean audio signal ($y_{clean}$), a noisy version of it ($y_{noisy}$) and an enhanced version ($y_{enhanced}$), by applying some speech enhancement algorithm to $y_{noisy}$.
What would be the right way to quantify the relative increase/decrease in the quality of $y_{enhanced}$ (compared to $y_{clean}$) in say, particular sections of the audio ? I'm not interested in SNR (or is it wrong to ignore it?)
I was thinking something like computing the spectrogram for both and compute the distance between them, say using $L_{\inf}$ or $L_2$ norm. I don't know if it's even correct way to do this. Any help would be much appreciated !
(PS: I'm from CS background, so I don't have much knowledge about signal processing. Apologies if this question seems too basic)
Thanks.

The answer to you question depends on how quality is defined in your specific application. Since you're working with audio signals, it is often beneficial to use the perceived quality to analyze the efficiency of an audio enhancement algorithm.
How to Determine the Perceived Quality
Methods to determine the perceived quality of an audio signal can be categorized into

subjective and
objective

methods. Subjective methods are based on human listeners that rate the audio signal according to some defined criteria. Objective methods use algorithms that use some measure, modeling parts of the human auditory system, to estimate the perceived quality.
Listening Tests
Usually, the most accurate way to determine the perceived quality is to perform a listening test and ask human listeners to rate the quality of different signals. Listening test methods comprise, e.g., the MUltiple Stimuli with Hidden Reference and Anchor and pairwise comparison approaches.
Objective Measures
A simple objective measure surely is the signal-to-noise ratio (SNR). In many cases, though, the SNR is only related to the actual perceived quality to a certain extent. Therefore, more involved objective measures have been created. Well-known ones are

the Perceptual Evaluation of Speech Quality (for speech signals) and
Perceptual Evaluation of Audio Quality (for audio/music signals)

algorithms.


Image Processing in MATLAB We try to multiplication operation of a black & white image with RGB image in same image but their class is not same: BW -> logical, RGB -> UINT8.

How can I convert from Logical to UNIT*? Or, can we do images which is in different class multiplication?

I guess you have a Binary (Logic) mask you want to multiply the image by.

It is easy:

mMask = double(mMaskBinary);

for ii = 1:numChannels

    mOutImage(:, :, ii) = mMask  .* mInputImage(:, :, ii);

end

In the latest MATLAB releases (From R2016b) expansion is done implicitly.
Hence if you want to do a multiplication between $ M \times N $ image to $ M \times N \times 3 $ image you can just multiply them element by element and the expansion to the 3rd dimension will happen automatically like in the loop above.

Enjoy.


Ideal Reconstruction of Upsampled Signal Problem: The signal $cos(2\pi14100t)$ is sampled at $F_s = 400 Hz$. It is then upsampled with a factor 3 and then reconstructed ideally with a new frequency $F = 500 Hz$.

I now want to find the new signal that is created.

Incomplete solution: With the frequency before sampling $\pm 141000 Hz$. After sampling the normalized frequencies:

$f = \pm {14100 \above 1 pt 400} \pm k = \pm (35 + {1 \above 1 pt 4})\pm k= \pm {1 \above 1 pt 4} \pm k $ per sample

Interpolation with factor 3 yields the normalized frequencies:

$f = {\pm {1 \above 1 pt 4} \pm k \above 1 pt 3} = $ 3 frequencies $\pm k?$

I know how to solve it past that to get the final $F= \pm f * F_s$. I.e. $y(t) = cos(2\pi f_1t) + cos(2\pi f_2t) + cos(2\pi f_3t)$.

So my question is how do I calculate the 3 frequencies? Would greatly appreciate a simple solution and not just an answer so that I actually understand this.

There are some mistakes in the question and your description.

When we sample a signal of frequency $f_m$ with a sampling rate $f_s$, the sampled signal contains the frequencies $f=f_m \pm nf_s$, where $n \in \mathbb{Z} $.
Here the frequencies available in the sampled signal are calculated as,

$f_m=14100 Hz$
 and $f_s=400 Hz$.
we know that,

$ f=f_m \pm nf_s $, where $n = 0, \pm 1, \pm 2, \pm 3, ... $

ie., $ f= ..., -1500, -1100, -700, -300, 100, 500, 900, 1300, 1700, ... $

When we up sample the signal by 3, the frequencies will become $f/3$.

So here the new frequencies are, $ f_u = ..., -633.33, -500, -366.67, -233.33, -100, 33.33, 166.67, 300, 433.33, 566.67, ... $

After upsampling the signal is reconstructed by using an LPF of cutoff frequency 500.
When we apply LPF to the above signal only the frequencies whose absolute value is below 500 remains and other frequencies will loss.

there for after applying LPF the frequency components available in the signal are,  $ -500, -366.67, -233.33, -100, 33.33, 166.67, 300, 433.33 $.

I hope this solves your problem.


How to generate AWGN with given power spectral density I am a little unsure about the generation of a additive white Gaussian noise, so I would appreciate to hear from you if I`m doing it correctly.

My goal is to generate the time-domain noise correspondent to a noise floor of -135 dbm/Hz

The procedure I`m adopting is the following:

1 - Convert the noise floor in dbm/Hz to Watts/Hz
2 - Get the average total power from the PSD.
3 - Generate the Gaussian noise with this average total power
4 - Take the FFT of the Gaussian noise and estimate the periodogram to verify whether it is indeed generating the noise with the correct PSD

Here is how I do the 4 steps above:

1 -
\begin{align} S_\text{watts} = \left(10^{\frac{S_\text{dBm}}{10}}\right)10^{-3} \end{align}
where $S_\text{watts}$ is the PSD in Watts/Hz and $S_\text{dBm}$ is the PSD in dBm/Hz.

2 - Average total power is given by $$ P_n = S_\text{watts} B ,$$ where $B$ is the Bandwidth

3 - Generate the randn MATLAB function multiplied by $\sqrt{P_n}$.

4 -
4.1 Take the FFT of the signal, let`s say, $X(k)$.
4.2 Get the Energy Spectral Density (ESD) as $\frac{\left| X(k) \right|^2}{N}$
4.3 Get the PSD by simply dividing the ESD by $N$ (the FFT length)
4.4 I do not simply plot the $10 \log_{10}\left(\text{PSD} 10^3\right)$, since I believe this would plot the PSD in dBm/subcarrier. Instead of this, I plot $10 \log_{10}\left(\frac{\text{PSD}}{\text{tone spacing}} 10^3\right)$, in attempt to get the plot as dBm/Hz.

The code below is what I am using. Please, indicate if I`m doing any mistake.

clear all
clc

bandwidth = 106e6;                          % in Hz
sampleFrequency = 2*bandwidth;              % in Hz
noiseFloor = -135;                          % in dBm/Hz
awgnPsd = (10.^(noiseFloor/10))*1e-3;       % in Watts/Hz
noisePower = awgnPsd*2*bandwidth;           % in Watts

% Generate noise signal:
N = 4096;           % signal length
nSymbols = 1000;     % number of symbols with length N to be generated
noise = randn(N,nSymbols)*sqrt(noisePower);

% Verify PSD
toneSpacing = sampleFrequency/N;        % DFT tone spacing
ESD = mean(abs(fft(noise,N)).^2)/N;     % Energy spectral density
PSD = ESD/N;

% plot PSD
figure
plot(10*log10((PSD/toneSpacing)*1e3))

Here is a Matlab snippet that might help. You should be able to apply the result to your favorite PSD routine to see the dbm/Hz in the frequency domain.

% create complex noise with a desired dBm/Hz
fs=1e6;             % whatever you want
N=1e5;              % # samples of noise
R=50;               % ohms
dbm_hz = -135;      % goal
bw = fs;            % complex noise with BW = sample rate
dbm = dbm_hz + 10*log10(bw);
dbw = dbm - 30;
watts = 10^(dbw/10);
vrms = sqrt(watts * R);
n = randn(N,1) + j*randn(N,1);
n = vrms * n / std(n);
meas_dbm_hz = 10*log10(sqrt(mean(n.*conj(n)))^2 / R) - 10*log10(bw) + 30


Equation that describes a gravity drained tank system I have a system like this:

The upper tank has an input ($q_i$) and an output, that is the input for the lower tank ($q_{12}$). The lower tank has an output ($q_o$).

Does anyone know how to get the equation that describes that system? I know that I will have two poles.

Gravity-fed outflow is proportional to the square root of the liquid level (see Bernoulli's equation ). Let's call $h_i$ the height of the upper (inflow) tank, and $h_o$ the height of the lower (outflow) tank. $q_{12} = k_{q_{12}} \sqrt(h_i) $ and $(dh_i/dt) = k_{h_i} (q_i - q_{12})$.

For the system outflow, $q_o = k_{q_o} \sqrt(h_o)$, and $(dh_o/dt) = k_{h_o} (q_{12} - q_o)$.

Valve assumptions: linear and fast. So the steady-state flowrate of $q_i$ is proportional to the valve opening via $k_{v_i}$, and the dynamics of the valve are much much much faster than the tank level dynamics, so no differential equation is needed - algebra only.

Under P-only control, $q_i = k_{v_i} k_c (h_{o_sp} - h_o)$

You might want to linearize those square root terms via Taylor expansion to make the Laplace transforms a bit less hairy, if you're working in the $s$ or $z^{-1}$ domain.

Of course, in industry (and assuming that the level of this tank is not critical - it rarely is) you'd just set it up to P-only control with a dimensionless gain $k_c = 1.67$. This way, that if the setpoint is 50%, the valve is wide open when the lower tank level is at 20%, and fully shut when the lower tank is at 80%. And then you'd move on to the next crisis.  :)


channel estimation for cdma I am looking for a method to make a multi-user channel estimation with CDMA.

Here is my question:

Is there a way to estimate multi-user channels with their spreading codes? If yes, what is the method?

Sure. The idea behind CDMA is that different users codes are orthogonal.

So, since your codes probably are pretty white, meaning that their autocorrelation function is essentially a dirac delta, you can simply encode an pulse for every user (or by every user), and use the standard correlation receiver to despread and thus recover each individual channel impulse response.


Strange noises when filtering audio signal I am using Naudio open source library and I am trying to do some simple filtering. The problem is that I hear some "clicks", not too loud. The library offers me the possibility to use at least two buffers, so the computing time doesn't introduce a delay between them.  Because in the most of the time I am dealing with a stereo signal, I have split it in two arrays and I compute each other independently. I would like to know if is there something special that I have to do with a filter when I used it on a buffer. I have first used a low pass biquad filter like the one below:

       //generate coeff
       //sincerely, I don't know what's up with q
       //I have taken into consideration some values
       //to see if the noise disappears
       double w0 = 2 * Math.PI * cutoffFrequency / _sampleRate;
       double cosw0 = Math.Cos(w0);
       double alpha = Math.Sin(w0) / (2 * q);
       _b0 = (1 - cosw0) / 2;
       _b1 = 1 - cosw0;
       _b2 = (1 - cosw0) / 2;
       _a0 = 1 + alpha;
       _a1 = -2 * cosw0;
       _a2 = 1 - alpha;
       for (int i = 2; i < length; i++)
       {
            output[i] = (float)((_b0 / _a0) * input[i] + (_b1 / _a0) * input[i - 1] + (_b2 / _a0) * input[i - 2]- (_a1 / _a0) * output[i - 1] - (_a2 / _a0) * output[i - 2]);
       }
       output[1] = (float)(
                (_b0 / _a0) * input[1] + (_b1 / _a0) * input[0] + (_b2 / _a0) * input[0]
                - (_a1 / _a0) * output[0] - (_a2 / _a0) * output[0]);
        output[0] = (float)(
                (_b0 / _a0) * input[0] + (_b1 / _a0) * 0 + (_b2 / _a0) * 0
                - (_a1 / _a0) * 0 - (_a2 / _a0) * 0);

I thought that  all my problems came from the first two samples (output 0:1), I've tried all combinations: output[-1]=0,output[-1]=output[0], but nothing worked. What values output[i-1], output[i-2] should have when "i" is 0 or 1?

I have encountered the same noise (clicks) when I used a LowPass Windowed-Sinc Filter, just like this:

//calculate coeff

        int i;

        int m = length;
        double PI = Math.PI;
        length=101;
        for (i = 0; i < length; i++)
        {
            if (i - m / 2 == 0)
            {
                _h[i] = 2 * PI * _cutOffFrecv;
            }
            else
            {
                //!=0
                _h[i] = Math.Sin(2 * PI * _cutOffFrecv * (i - m / 2)) / (i - m / 2);
            }
            _h[i] = _h[i] * (0.54 - 0.46 * Math.Cos(2 * PI * i / m));
        }
        //normalize the low-pass filter kernel for unity gain at DC
        double s = 0;
        for (i = 0; i < m; i++)
        {
            s = s + _h[i];
        }
        for (i = 0; i < m; i++)
        {
            _h[i] = _h[i] / s;
        }
        //convolve the input & kernel
        //_kernelSize=101
        //most often length is 6615 or 6614 for each channel
        //in these examples I compute only one channel
        for (j = 0; j < length; j++)
        {
            output[j]=0;
            for (i = 0; i < _kernelSize; i++)
            {
                if (j >= i)
                {
                    output[j] =(float)(output[j]+ _h[i] * input[j - i]);
                }
            }
        }

The problem surely is not from splitting the signal, or from combining channels, because I have tested this without any filter and everything is ok. I have also tried to simulate some delays created by a processing algorithm (but without changing the signal) and nothing went wrong. I am very sure that the problem comes from filtering.
Everything  I wrote is used on a buffer.

The context in which you use these functions is not clear, but it seems to me that your problem is "edge effects".

When you are evaluating the convolution or the biquad, you need to access samples which are outside the current buffer. Your two implementations evaluate these samples as zero. This is incorrect. For example, for the biquad, everytime you process a block of audio, you need to store the last 2 values of the input[] and output[]; and reuse them in place of input[-1], input[-2], output[-1], output[-2]. Even if the data you process comes in small chunks, you must process it as if it came in one single stream; so the state variables of your filters must not be reset to zero at the boundaries of each buffer.


Digital PLL loses lock every time mesage crosses zero I think my digital PLL is almost ready, but I've encountered this problem and I don't know what might be causing it, so I'd appreciate if you could help.

I'm using the DPLL in a Costas loop. It seems to work fine but whenever the message crosses zero, lock is lost. After that, it ends up locking again. You can see here:

Input to the Costas loop is:
$$
x(t) = m(t)cos(2\pi f_ct+\theta_{in})
$$
$$
m(t) = cos(2\pi f_mt)
$$
where $f_m << f_c$

Not sure if this is a known phenomenom (it probably is) but I haven't been able to find anything related. It'd be nice if you could point out what's going on and/or how to fix the temporary loss of lock?

One way I've managed to get it "working" consisted of squaring the message signal (although phase error rises a bit when message goes near zero), but I don't wish to square it so that wouldn't be of use to me.

This is the Costas loop I'm implementing for this test:

You are having trouble because that's not a Costas loop.  A Costas loop uses demodulated data in some form to change the phase that's expected from the signal.

You're just taking the I/Q demodulated signal and applying it to the atan2 function; that makes a sort of linearized extended phase detector, but without determining that the phase should have flipped by 180 degrees, your loop is naturally trying to follow the phase that's actually coming in.

You could make your loop into a simple Costas loop by changing the atan2 function into a multiply, and following it by a proper loop filter.


Compensating Loudspeaker frequency response in an audio signal I have been working on a project in which I was required to work on the audio signals recorded from the loudspeaker kept in front of a filter. So, to simply explain it:

$$\boxed{\rm LoudSpeaker} \longrightarrow \boxed{\rm Filter} \longrightarrow\boxed{\rm Microphone}$$

Now the project depends on finding how the filter reacts to the signals from the loudSpeaker. The loudSpeaker outputs a frequency sweep.

When I saw the datasheet of the loudspeaker I found that it had a certain frequency response which was necessary to compensate within the sweep. But now that I have already taken the readings into the microphone I have to subtract the loudSpeaker frequency response from the final spectrum.

A possible way to that is Deconvolution. But I can only explain that in theory.

Can someone help me how to implement these functions in MATLAB? Or a different way to solve this problem???

Yes, you can do this with an LMS equalizer which uses the Wiener-Hopf equation to determine the least squared solution to the filter that would compensate for your channel, using the known transmit and receive sequences. The channel is the unknown being solved, and the tx and rx sequences are known.
BOTTOM LINE:
Here is the Matlab function with error checking removed:
function coeff = equalize(tx,rx,depth,ntaps)
%Determines equalizer coefficients using the Wiener-Hopf equations
%TX = Transmitted (Desired) waveform, row vector, length must be > depth+2*ntaps
%RX = Received (Distorted) waveform, row vector, length must be >=depth
%DEPTH = Depth of solution matrix (recommend 10x ntaps but based on duration of stationarity)
%NTAPS = Number of taps for equalizer filter

%force row vectors
tx= tx(:)';
rx= rx(:)';

delay=floor(ntaps/2);
A=convmtx(rx(1:depth).',ntaps);
R=A'*A;
X=[zeros(1,delay) tx(1:depth) zeros(1,ceil(ntaps/2)-1)].';
ro=A'*X;
coeff=(inv(R)*ro);

USE:
Once the coeff for the FIR filter are determined using the function above, then the Matlab filter function can process the receive sequence:
tx_recovered = filter(coeff, 1, rx)
If you want to see the channel response of the filter use:
freqz(coeff)
If you want the solution to be the estimate of the channel instead of the compensation filter that undoes the channel response, simply swap tx and rx:
coeff = equalize(rx,tx,depth,ntaps)
DETAILS FOR THE VERY INTERESTED:
See my slides below giving a high level overview / derivation of the process, This in general form is the Normal Equation (http://mathworld.wolfram.com/NormalEquation.html) used for least squared curve fitting and other applications. I believe I was first introduced to this viewpoint in demonstrating how the Normal Equation is performing deconvolution from the book "Theory and Practice of Modem Design" by John A.C. Bingham.
In practice, I typically do a cross correlation first to determine the channel response time (delay spread) and initial time alignment, and then use an initial equalizer FIR length (# of taps) that exceeds the delay spread (not knowing if leading or trailing echos dominate I will typically start with 2x the delay spread for the FIR length). Once I see the result, the filter size can be reduced if desired based on insignificant magnitudes of the coefficients at the edges of the filter.  If the sequences are not exactly aligned, but still within the span of the filter, then the dominant tap will be offset accordingly- so not critical to align beforehand and this gives you insight into what happens if they are grossly misaligned.

Here is an interesting example of the equalizer function I used recently on a sound file from Dalen to equalize the waveforms received by the left and right channels as received by two microphones (treating left as transmit and right as receive and ignoring the actual third party transmitter for the two). The two channels are not recognizable prior to equalization, and completely aligned in amplitude, phase and characteristic after.
Here is a plot of the left and right channels prior to equalization:

Here is the same plot after equalization, right was filtered with the equalizer, and left was filtered with a simple filter just as long as the equalizer with a single unity gain tap in the center and zero elsewhere (to match the delay as the equalizer assumes nominal delay is in the center of the equalizer filter):

This is a zoom in plot of the waveforms after equalization showing how identical the two sequences have become:


Advantages of an orthogonal filterbank What are the advantages orthogonal filterbanks such as the DFT or DCT have over non-orthogonal filterbanks?

I have heard quite often that they have attributes that are desirable when modifying the intermediate signal, i.e. when performing quantization during audio coding.

Orthogonality provides an interesting backbone to the structure of the filter-banks (FB). First, from an analysis FB, the synthesis FB is very direct, so it can ease implementations. Second, the orthogonality often allows faster implementations, as there is "little redundancy" in computation. Third, orthogonality ensures that matrices are well-conditioned, this may reduce risks of error propagation (esp. with quantization). Fourth, energy is preserved, so operations performed in the dual domain can be measured (energetically) in the primal domain. Fifth: it's a great deal of simplification when you need to get tractable proofs for their design (energy compaction for instance) or performance (as for instance Gaussian noise remains Gaussian after an orthogonal transformation). Note that these five points are not completely orthogonal.

To provide a little more details:

Orthogonality, and more generally paraunitarity (cf. Paraunitary matrices by Barry and Ted Hurley), i.e. a square matrix $U(z)$ satisfying $U(z)U^∗(z^{−1}) = 1$, are practical, as the inverse is somehow a conjugate-transpose
For orthogonal bases in an $N$-dimensional space, the first orthogonal vector has about  $N-1$ degrees of freedom, the second one $N-2$, etc. So orthogonal matrices have a structure (see the Orthogonal group), which can be used to factorize them, hopefully with less operations (Fast multiplication of orthogonal matrices). For instance, $f 3 \times 3$ orthogonal matrices can be represented by  three Givens rotations, but other designs are possible
Determinants are $\pm 1$ (and of course the same of their inverses). When a matrix contains close vectors, projections split evenly of them, which has effects on 1) simple inversion 2) iterations (like in computerized tomography) where powers of matrices (like FFT)  over millions of cells/pixels/voxels may diverge.
When energy is preserved, a structured signal (in the appropriate basis or frame) may come out of unstructured noise. This is a rationale behind "sparsifying transformations", and how thresholding values can denoise. Somehow behind Stein estimates (Stein Unbiased Risk Estimation), with orthogonaly invariant estimators, sample eigenvalues are modified (for the best) yet sample eigenvectors untouched.
All block orthogonal transforms are quite easy to deal with.

a) Overlap and windows are more involved. An early attempt is the LOT (Lapped Orthogonal Transform), and their generalizations (GenLOT - Generalized Lapped Orthogonal Transforms , GULLOT - Generalized Unequal Length Lapped Orthogonal Transforms , etc.). They can be seen as DCT (or Hadamard) matrices with additional degree-of-freedom. By decomposing them with rotations, it is possible to optimize their design by imposing constrains, like "Coding gain", "DC leakage" (Linear phase paraunitary filter bank with filters of different lengths and its application in image compression, 1999), that can be related to their performance in data compression. Trivia: read Properties of the Eigenvectors of Persymmetric Matrices with Applications to Communication Theory to see how correlation structures impose properties of eigenvectors
b) Some noise properties can be estimated in a data sparsifying domain. A classical example in the wavelet domain (and filter-banks generally) is the MAD (median absolute deviation) estimator of the Gaussian noise in a combination of subbands.

However, orthogonality has limits, imposes additional constraints that sometimes spoil other desirable properties. Plus, the human sensory systems care little about  strict orthogonality, and you can gain degrees of freedom (even in computational efficiency) by discarding strict orthogonality. The filterbanks used in the JPEG2000 image compression standard are not orthogonal, but biorthogonal. [UPDATE] In some cases, you may want integer or dyadic rational operations (for memory or speedup), and this is rarely compatible with orthogonality. In video coding, a lot of decorrelating operations are not strictly orthogonal.

On most practical cases though, a lot of systems remain at least "close to orthogonal" in some sense.


Separate Simple Signal into Exponential Components Everyone here should be familiar with the trignonometric identity

$$2 \cos(x)=e^{jx}+e^{-jx}$$

Let's say I have a real world analog signal given by the following function

$$x(t)=2\cos (\omega_1 t)$$

and I only want a single exponential component, specifically

$$y(t) = e^{j\omega_1 t}$$

where $\omega_1$ is a real constant. How to do this? Is it possible?

Your problem seems to be finding the analytic signal to a given real signal.

What you apparently need is the Hilbert transform. This will allow you to "recreate" the imaginary part of $y(t)$ from the real part alone, which is essentially your measurement $x(t)$.

Note that the transform can only be approximated as it is a acausal filter, and only works best when the signal $x(t)$ is available in the computer for all $t$.

A simple implementation is to Fourier transform $x(t)$ and set the negative frequency components to zero. After transforming back you obtain (the complex signal) $y(t)$.


inverse z transform performed on 6th order IIR filter we are told to find coefficients and impulse response of IIR filter of order of 6. There are 6 zeros and 6 poles in the design. Pole and zero pairs are conjugate and poles are within the unit circle whereas zeros are on the unit circle and. I can write transfer function with no problem, but I cannot expand it as its order is too big to work with. I must work my way up to finding impulse response and filter input/output coefficients, but I am stuck terribly with the transfer function. What can I do? long divisions are used in second order polynomials, no example of it for larger orders. Is there a easy way? I am up for hours trying to figure it out.

You can write your transfer function as the sum of three 2-nd order transfer functions (partial fractions), and then the impulse response will be the sum of the impulse response of each.
To help you with the partial fractions you can use this calculator


Signal compare - cross corelation in frequency domain matlab I need compare two signals. Signals have simliar shapes but completely different values (as on picture):

So I think cross corelation in frequency domain should gave me a number that shows me how the signals are similar. But I have no idea how to do that in matlab.

Where I can find code to do it?

Yes, cross-correlation will provide some information on the similarity of the two waveforms. See the function reference for xcorr. Look into proper normalization so that the result is between something like 0 and 1 or -1 and 1. As a matter of fact, this is pretty well discussed in another stackexchange question.

As for frequency domain cross-correlation, I don't think that's what you want to do. Doing time-domain cross-correlation results in multiplication in the frequency domain, so you're looking to maximize the total energy in F(x)*F(y), where F(x) is the FFT of x.


evaluating individual terms of 2D DFT (Goertzel?) I'm looking for a fast way to do a 2D discrete Fourier transform of an image at many arbitrary frequencies. I know the Goertzel algorithm works for 1D, but is it possible to generalize it in 2D? Or any other useful algorithms? Two solutions that I know, but which I'm afraid might be too slow, are directly using the DFT equation, and zero padding and doing an FFT. Alternately, the frequencies might not need to be completely arbitrary but could be sub-arrays of the array of frequencies given by an FFT with high frequency resolution.

The 2D DFT is separable, so you can conceptually compute it by applying 1D DFTs in both directions. For the Goertzel algorithm you can first apply it for a specific column for each row, and then use that column and apply it again to obtain the individual term. Most likely, you can rewrite the Goertzel algorithm to do this in one go, but the computational savings are quite small (for an $N \times N$ image, you would need $N^2 + N$ multiplications using the separable approach, while a 2D version would give $N^2$ multiplications).


Show that any continuous-time signal $x(t)$ can be represented as $x(t)= x_e(t) + x_o(t)$ Show that any continuous-time signal $x(t)$ can be represented as $x(t)= x_e(t) + x_o(t)$

where $x_e(t) =\frac{1}{2}[x(t) + x(-t)]$ and $x_o(t) = \frac{1}{2}
[x(t) − x(-t)]$
are even and odd functions, respectively.

Substitute the last two equations into the first and it should be clear.

iI.e., if you have $x_e(t), x_o(t)$  given to you. Substitute these into the equation for $x(t)$ and you should see that both sides are equal to each other.

(Note that this proves the equality, but doesn't prove that $x_e(t)$ is even or that $x_o(t)$ is odd.)


Is Cross correlation between stationary random signal and noise zero for all lags? If $x(n)$ is a stationary random process (let say small 20ms window of speech which can be assumed stationary) and $v(n)$ is said to be uncorrelated noise, can we say that cross correlation of $x(n)$ and $v(n)$ is always zero, for all lags? E.g.

$$r_{xv}(\tau) = E\{x(n) v(n+\tau)\} = 0, \forall \tau$$

Or,

$$\sum_{n=0}^{N}x(n) v(n+\tau) = 0, \forall \tau$$

where N is the signal length?

As Marcus alludes to in the comments on the question, you are using two different definitions of cross-correlation.

The first definition of cross-correlation, the statistical one, is:
$$r_{xv}(\tau) = E\{x(n) v(n+\tau)\} = 0, \forall \tau\tag{1}$$
where, the expectation operator means:
$$
E\{ z \} = \int_{-\infty}^{+\infty} z p_z(z) dz \tag{2}
$$
where $p_z(z)$ is the probability density function of the random variable argument of $E\{ \cdot \}$, z in equation (2).  Equation (2) is sometimes called the ensemble average of $z$.

It is true that if $x$ and $v$ are uncorrelated, then (1) is true.

The second definition in equation (3) (or the normalized version of it) is what we usually use because we don't know what $p_z(z)$ is.  In order to use the second definition, we have to make various assumptions about the relationship between $x$ and $v$, one of which is ergodicity (that time-averages of one realization of each can be substituted for the ensemble averages).

$$\sum_{n=0}^{N}x(n) v(n+\tau) = 0, \forall \tau \tag{3}$$

The normalized version is:
$$\frac{1}{N}\sum_{n=0}^{N}x(n) v(n+\tau) = 0, \forall \tau \tag{4}$$

Now we are using time averages to find the cross-correlation.  For any given realization of $x$ and $v$, we are not guaranteed that equation (3) ( or (4)) is true.

However, what generally happens is that the mathematics surrounding any theory using cross-correlation uses (1) while the implementation uses (3).

That is why implementations of cross-correlation, such as xcorr in Matlab, do not output a vector of zeros with uncorrelated input time series.

The confusion about cross-correlation seems to be widespread. As @msm notes in the comments, the Wikipedia page for cross-correlation makes no mention of its statistical origins.

The whole point of correlation (either auto or cross) is to see how similar one signal is to another, at various time-lags.  For the auto-correlation case, the aim is to see how "predictable" the signal is --- knowing the values up until sample $n$, how much can I say about the sample values after time $n$ before they occur?

Let's start with a Wikipedia page that hasn't forgotten these origins.

This page starts with the original version of normalized correlation:
$$
\rho_{XY} = \frac{\textrm{cov}(X,Y)}{\sigma_X\sigma_Y} = \frac{E[(X - \mu_X)(Y - \mu_Y)]}{\sigma_X\sigma_Y}
$$

This is another point of confusion for many: staticians generally require that the correlation is normalized to being between -1 and +1. Signal processing engineers tend to do away with this requirement.

The trouble with $\rho_{XY}$ as defined above is that, as in the earlier part of my answer, taking the expectation requires knowledge of the statistical properties of the signals, which we often don't have or have to guess.

That is why the actual cross-correlation is usually substituted for the sample cross-correlation:
$$
r_{XY} = \frac{\displaystyle\sum_{n=1}^{N} (x_n - \bar{x}) (y_n - \bar{y})}{n s_X s_Y}
$$
where $s_X$ and $s_Y$ are the sample standard deviations and $\bar{x}$ and $\bar{y}$ are the sample means.


Amplify a signal and phase shift it by multiplying by a complex number Still need help 25th October 2019

I have a real time domain signal contained in an array, it's just two different frequencies summed together. If I want to amplify this signal by 10 and phase modulate it by 0.5 degrees, will this work?

Change the real signal to a complex signal. $ (e^{\omega_1t} + e^{\omega_2t})$
Multiply it by the phasor $Ae^{i \theta}$
So the calculation is $Ae^{i \theta} (e^{\omega_1t} + e^{\omega_2t})$
Where $A=10$ is an amplitude change and $\theta=0.5 \pi/180$ is the phase change caused by the equipment
Take the real signal of this calculation to get the result

The numbers I have used, 10 and 0.5, are just an example.

Please can someone share knowledge please?

Multiplying by a complex number causes a rescaling by the magnitude and a rotation by the angle.  This is why when you multiply two complex numbers, you multiply the magnitudes and add the angles.

The only signal for which a complex multiplication (rescale and rotation) is the equivalent of a phase shift is a single complex pure tone (looks like a spring, or threads on a bolt).  The frequency is the translation value between the rotation amount and the shift size in the time domain.  Thus if you have a signal that is a sum of complex pure tones with different frequencies, each will be shifted by a different amount by the same non-real multiplication factor, and the sum of the shifted tones will not match a shift of the sum.

Since a real pure tone is the sum of two complex pure tones of different frequencies, it cannot be shifted by a complex scalar multiplication, nor can any sum of them.

That's all there is to it.


Classifying sleep stages from only EEG I'm currently working on a project that requires me to classify sleep stages (Awake W, N1, N2, N3 and REM) based on only an EEG. Various algorithms and classifying standards (such as Rechtschaffen & Kales or AASM) actually require EMG and EOG in order to distinguish from REM sleep and wakefulness, which I don't have. Additionally, most algorithms use techniques such as Wavelet analysis, ICA, Neural Networks and such, but I have only learned up to Fourier Transform. Which means that I have to detect and classify these stages and their components (detect sleep spindles in order to classify N2, for example) based on only frequency analysis.

Power spectra analysis is poor since low frequency components tend to have high powers so the program ends up over-classifying 30s epochs as deep sleep stages.

I have a .txt file telling me exactly what's going on in every 30s epoch, so I have a priori knowledge regarding when is the subject awake or asleep (and which sleep stages the subject is in).

I appreciate any insight into this problem.

The standard approach would be to

Choose your feature vectors.
Partition your data set into training and testing.
Calculate the feature vectors for all your data (both training and test).
Use the training data set to train your classification algorithm.
Apply the trained classification algorithm to the test data set and calculate your correct classifications and incorrect classifications to get your precision / recall.

The trick will be step 1: choosing your feature vectors, given that you've only got Fourier coefficients (and values derived from them).


Minimum Phase & maximum phase i have 2 questions:
1-If zeros on the unite circle that will be min or max phase?and this system has low-pass frequency response or not ?
2-Is the min phase has the shortest impulse response?

Let me first state the commonly used definitions of minimum-phase and maximum-phase for discrete-time systems:

  A causal and stable discrete-time system is said to be strictly minimum-phase when all its zeros are inside the unit circle.

Of course, all its poles are inside the circle as well, otherwise it wouldn't be causal and stable.

  A causal and stable discrete-time system is said to be strictly maximum-phase when all its zeros are outside the unit circle.

Again, causality and stability requires that all the poles are inside the unit circle as well.

If you have zeros on the unit circle (and all others either inside or outside the circle), the system is usually referred to as loosely minimum-phase or loosely maximum-phase. A system with all its zeros on the unit circle is a degenerate case of that definition, in the sense that it is loosely minimum-phase and loosely maximum-phase at the same time. But in that case these definitions do not make much sense. These definitions are meant for systems that have some zeros inside (or outside) the circle and some on the circle.

Your question about the low pass frequency response can't be answered, because the concepts of minimum-phase and maximum-phase have absolutely nothing to do with the magnitude response of the filter. You can have a minimum-phase low pass filter, and a maximum-phase low pass filter with exactly the same magnitude response. The same is true for all other frequency selective filter characteristics (high pass, band pass, etc.).

Concerning your last question, it depends what you mean by shortest. If you mean the number of samples of the impulse response, then the answer is no. As an example, take a length $N$ FIR minimum-phase impulse response. If you flip the coefficients (i.e., time reverse them) you obtain a maximum-phase filter of the same length with the same magnitude response.

The minimum-phase property means that of all filters with the same magnitude response, the impulse response of the minimum-phase filter is maximally concentrated near time index $n=0$ in a mean square sense. I.e., if $h[n]$ is a minimum-phase impulse response, then

$$\sum_{n=0}^{N}|h[n]|^2\ge \sum_{n=0}^{N}|g[n]|^2\tag{1}$$

for all $N\ge 0$, and for any $g[n]$ with the same magnitude response as $h[n]$.


Pole Magnitude and Damping Ratio relationship I know that the damping ratio of a system is defined by the angle of the pole, calculated with respect to the left hand side $x$-axis. Could one infer though, that if the magnitude of the poles is small (i.e. the conjugate poles lie closer to the origin), then the damping ratio of the system would be small ,and result in a slower system response?

The magnitude of the poles alone cannot determine the damping ratio, which in turn is related to the rate at which the oscillations about the steady state value will decay.
In fact it is the ratio of the absolute value of the real part of the poles and the magnitude of the poles that gives the value of $\zeta$.
If you think about a second order system its characteristic polynomial would be:
$$
s^2+2\,\zeta\,\omega_n\cdot s \,+\, \omega_n^2 = \\\left(s\,+\omega_n\cdot \zeta-j\,\omega_n\,\sqrt{1-\zeta^2}\right)\cdot \left(s\,+\omega_n\cdot \zeta+j\,\omega_n\,\sqrt{1-\zeta^2}\right)
$$
with poles: $p_{1,2}\,=\,-\omega_n\cdot \zeta \pm j\,\omega_n\,\sqrt{1-\zeta^2}$
However, the above pole forms are valid for $0 \le \zeta \lt 1$. These systems are called underdamped.
Whereas the damping ratio is:
$$
\zeta = \dfrac{|\Re(p)|}{\omega_n}
$$
As a result of the above the magnitude and the damping ratio relate to each other through:
$$
\zeta = \dfrac{|\Re(p)|}{||p||}
$$
You are correct in that the angle of the pole, with respect to the left hand side real axis relates to $\zeta$, and this is the formula:
$$
\cos{\theta} = \dfrac{|\Re{p}|}{\omega_n} = \dfrac{\zeta\cdot\omega_n}{\omega_n} = \zeta
$$
The conclusion that the magnitude alone, whether small or big, cannot determine the decay rate of the oscillations is also apparent from the fact that if we normalize the poles with the natural frequency $\omega_n$ their magnitude will be always equal to 1, this means that if the angle of the pole somehow varies, then they will always be on the unit circle. In such case as the angle varies the damping ratio will also vary but the magnitude will always be the same no matter what the angle.
So its the real part of the normalized pole that determines the damping ratio
rather than its magnitude, alone.
Note: A slow, from a steady state point of view, pole is in fact one that has a small real part, while a fast is one that has larger real part, in absolute value.
Systems with $\zeta = 1$ are critically damped, having conjugate imaginary poles with magnitudes at the system's natural frequency:
$$
p_{1,2} = \pm\, j\, \omega_n
$$
So this means that the hypothetical system has a pair of conjugate poles with equal magnitude:
$$
||p_{1,2}|| = \omega_n
$$
Those systems may exhibit oscillatory behavior sustained longer than a underdamped system, for example, a sinusoidal response. The poles angle is $90^o$ degrees because their real part is zero. Hence $\zeta = \cos(90^o) = 0$ implying an undamped system but the magnitude of the poles is $\omega_n$, thus it depends on the natural frequency of the system which could be huge in value.
In many engineering applications, the measure that is used extensively for pole quality characterization is the $Q$-factor.
In fact, the quality factor and the damping ratio relate to each other through the following formula
$$
Q = \dfrac{1}{2\zeta}
$$
In addition, you may come across transfer functions with characteristic polynomials parameterized with $Q$ instead of $\zeta$.
$$
s^2+2\,\zeta\,\omega_n\cdot s \,+\, \omega_n^2 =
s^2+\dfrac{\omega_n}{Q}\cdot s \,+\, \omega_n^2
$$


Maximum Likelihood Estimator (MLE), MMSE and LS - Are All of Them Regressor, Estimator and Predictor? Can all three criteria ML, MMSE, and LS be called regressor, estimator, and predictor ?

If not, an intuitive explanation of why they can't be, would be good.

Let's see:

Maximum Likelihood (ML) - A Method to estimate parameters (MLE - Maximum Likelihood Estimator). Given a PDF of samples, what is the parameter value which maximizes the probability of the observations.
Least Squares (LS) - Least Squares model for estimating parameters. Sometime coincides with the MLE (Usually in the case of Gaussian Noise).
Minimum Means Square Error (MMSE) - Bayesian Estimator based on Loss Function which is the second moment of the error.

All of them are estimators in a model (PDF of Parameters in the ML, Numeric Model in the LS or Bayesian Model for the MMSE).

Predictor is a term for predicting future data based on given data. Any of the method above can do that with different assumption of the model.
Regressor is a term for predicting model parameters from a given data. Again any of the above can be used to do so.

So you need to separate the means you use (Type of estimator and model) for the object you try to achieve (Predict, Regress, Smooth, etc...).