Downsampling a PSD that was generated by zero-padding - worth it? Let us say we have a signal $y(k)$ that is of length $K$. I can compute its PSD, also of length $K$. Let us say that this length was not good for finding peaks, so I computed a PSD that was of length $10 \ K$ instead, by doing a zero-padding in the time-domain. I understand that this is equivalent to interpolating the frequency domain result.

So my question is, at this point, if I resample the actual $10 \ K$ length PSD back to be $K$ samples long, would it give me the same result as if I took a $K$ length PSD to begin with? If not, is there any advantage to doing it this way, or would I just be better off sticking to taking a length $K$ PSD to begin with?

Consider these two cases where you operate on a signal $y[k], K = 0, \ldots K-1$:

You take the DFT of the signal $y[k]$ to obtain a length-K frequency-domain signal.
You first zero-pad $y[k]$ to obtain a longer signal $y_p[k], k = 0, \ldots, 10K-1$. You then take its DFT. Finally, you downsample the DFT by discarding all but the first of every 10 samples (i.e. you keep the first, tenth, twentieth, and so on).

The two processes will yield the exact1 same result.

1 In practice, that is, when using finite-precision arithmetic, it's likely that you won't arrive at results that are identical bit for bit. However, in theory, and with the use of infinite arithmetic precision, the two approaches are equivalent.


What characterizies 'causality' for a finite FFT? Causality of a LTI transfer function $G(\tau)$ in the continuous time domain, i.e. for
$$y(t)=\int G(\tau)x(t-\tau)d\tau$$
is characterized by
$$G(\tau < 0) = 0$$
By the way, in the frequency domain, causality is characterized by the Kramers-Kronig relations, which is why I am particularly interested in the Fourier transform $G(\omega)$.

Now I am analyzing the functions $x(t)$ and $y(t)$ in discrete time and for a finite set of equidistant samples (measurements). In order to work in the Fourier domain, I will have to assume that they are periodic functions (or rather series):
$$x(t+T)=x(t) \qquad y(t+T)=y(t) \qquad G(t+T)=G(t)$$
Even with such techniques as zero-padding, they are still periodic as seen from the view of FFT.

But what means causality for periodic functions? If $G(t)$ is the response of a periodic Dirac comb, it is not clear if it is the causal response of the Dirac impulse at $t=0$ or the non-causal 'response' at some of the later Dirac impulses at $t=nT$ where $T$ is the measurement duration. Strictly speaking, it is the response to all the Dirac impulses in the comb and it makes no sense to define the notion of 'before or after' the impulse.

Am I missing something? Is there a proper notion of causality for periodic functions and hence on finite domains?

I am interested in this question and since no-one has answered yet I will offer my further thoughts (up for debate or confirmation):

I agree with the OP's perspective that the Discrete Fourier Transform (DFT) fundamentally represents when cast to domains that extend to $\pm \infty$ periodic sequences in time (and in frequency) and therefore is fundamentally non-causal due to the time periodicity specifically.

However, to utilize these transforms with context of causality, we can consider the cases that best match an approximation of the Continuous Time Fourier Transform for the same causal representation of our sampled system. Such cases are the ones with no time aliasing, which would occur when $x[n]=0$ for $n\ge N/2$ given $N$ samples with $n \in [0, N-1]$. When we can restrict ourselves to causal sequences specifically (meaning when we know the underlying process is causal) then this restriction can be extended to be $x[n]=0$ for $n\ge N$ given $N$ samples with $n \in [0, N-1]$. This means basically that we know the time duration of the data capture exceeds the expected response time of the underlying continuous time process.

Further in frequency for such an equivalently causal time domain waveform the phase would be going increasingly negative versus frequency indicating the delay that is indeed properly being modeled.

It is the counter case, when for a known causal sequence $x[n]$ has significant non-zero values in the upper-half of the time domain samples right up to the $N-1$ boundary that we cannot be reasonably assured we are not observing the effects of time domain aliasing. We are no longer able to uniquely distinguish between causal and non-causal sequences, or more specifically and practically when we know the sequence is causal; we are no longer able to distinguish from underlying continuous-time responses that went beyond the sample duration of the sequence due to time domain aliasing (so specifically in the OP's case there is no way to know if the solution is representing longer response times than the samples are providing).

This is not very different from what I would typically do when evaluating results in the frequency domain: if I observe strong spectral content right up to the Nyquist boundary I am not assured that my sampling rate is high enough (or filtering is not tight enough) as frequency domain aliasing is likely occurring. While if the spectrum rolls off to sufficiently low levels prior to the Nyquist boundary (and if I can assume the system in question would not have higher spectral content at further offsets) then I am reasonably assured that the digital representation accurately represents the continuous time spectrum. This obviously requires knowledge of the system and what the likely spectral content is beyond the sampling rate, but is very analogous to the time domain challenge presented here.

For further details on matching an underlying CTFT with the DFT scroll down to the "Exact CTFT Result from DFT" section of my answer at this post: Why do we have to rearrange a vector and shift the zero point to the first index, in preparation for an FFT?


How to convert a filtered time-domain signal into frequency one? I have a signal in matlab like this:

I have designed a band-pass filter in 16-32 MHz range:

and applied the filter to the signal,

the result is like this :

My Question:

what is the easiest way to convert this time-domain signal into a frequency based signal?

can I export the signal into workspace and use the data variable?

Taking the FFT will give you a frequency representation of the data:

y=fft(sig3.data);


Options for DSP hardware for Under Grad Lab? We have TI 6713 dsp kits in our Lab but since those kits are very old(They were started/released near 2000) we are now considering to get new dsp kits ,so our under grad students may clear their hardware concepts related to dsp in a better way
I have seen a question on DSP SE that is some how relevant and  especially the answer of Hilmar there.
But am i bit confused?
Is it okay to use kits like Rasberry PI or Aurdino,for teaching dsp Lab to under grad students?especially as alternative option to dsp kits of TI or Analog devices
If possible,Please also recommend other hardware options
Link of relevant question
Hardware kits for DSP Lab?

I think using a Raspberry PI or Arduino (without ADCs or DACs) has no advantage over using just a PC with Matlab or Python. The thing that separates DSP hardware from "normal" processors is that they can do certain operations very fast and/or in parallel. If you want to do something cool, I would use a FPGA with ADCs and DACs like the Red Pitaya for lab experiments.


Does a power spectral density that follow a power law imply that the time hisory follows a power law? This is a question related to this previous question "I cannot find reference (paper) of this relation $u(t)t^{α} ↔^{FT} f^{-(α+1)}$".
The signal that I work with g(t), is the time history of a velocity component flucutuations. This signal comply with the rule:
$\bar{g}=lim_{T->+\infty}\int_{0}^{T}g(t)dt=0$
The Power spectrum of g(t) follows in a frequency range a power law "-5/3", in other frequency ranges the  slope can change. (SEE for example, figure 1 [here])2.
The power law frequency range, following the rule $u(t)t^{α} ↔^{FT} f^{-(α+1)}$ should imply that also the time history has a power law. The signals I work with does not seems to have a power law.
Is this can be due to a masking performed by the behaviour of the other frequency range? The rule
$\bar{g}=lim_{T->+\infty}\int_{0}^{T}g(t)dt=0$
can violate some hypotheses at the base of $u(t)t^{α} ↔^{FT} f^{-(α+1)}$" ?

should imply that also the time history has a power law

No. The PSD is an incomplete description of the signal since it has lost the phase information. There are infinitely many time sequences that have the same PSD. For example a unit impulse, white noise and a linear sweep have the same PSD (roughly speaking) but are completely different time-domain wave forms.
In your case $u(t) \cdot t^{\alpha}$ is one signal that fits your PSD but there are many others that will fit too but have a totally different shape.
This being said, your observation of a conflict has some merit. Your first equation simply states that $g(t)$ is mean free which in turns implies that the PSD at $f=0$ must be zero, but $0^{-5/3}=1$. I'm guessing that your power law applies only to a limited frequency range which doesn't include $f=0$ and hence there is no conflict.


Adjustment of window length with sliding/moving RMS method for frequency drifts I'm using a sliding RMS mechanism to compute RMS of a signal, i.e. with a window size of N, I add the squared value of the new sample $x[n]$ to the running total while deducting the squared value of the $x[n - N]$ sample and then performing a square root to get the RMS value.

This method works fine when a cycle's worth data in the signal under question lines up exactly with the length of the window. For ex,

Let the window length $N = 512$, frequency of the signal $f=100 \ \mathrm{Hz}$ and sampling frequency $f_s= 51200  \ \mathrm{Hz}$.

With the above settings, every set of 512 samples would have a cycle's worth data of the signal and the computation works fine.

Now, the issue begins when the frequency of the signal deviates from $100  \ \mathrm{Hz}$ giving more or less number of samples in the same 512 samples window which throws off the accuracy of the sliding RMS algorithm.

This happens because the division by 512 to get the RMS, is no longer valid due to the fact that this particular cycle of the signal is getting complete before the end of the 512 samples window. However, the running squared total would already have the contributions from the previous samples which were divided by 512 while the samples of this cycle would require a different window length.

Is there a way to compensate for the addition/lack of samples due to the frequency drifts while using this technique?

You're right that when you're using this method to calculate the RMS value of a sinusoid, it's sensitive to getting an integer number of periods of that sinusoid in the observation window. I can see a couple ways to tweak your approach:

Since it's important to get full periods into your estimate, you could add some tracking of the sinusoid's frequency and phase. In the most complicated implementation, this could be a full-blown phase-locked loop. I've seen embedded RMS measurement systems before where this approach is used. Basically, you create  a PLL that is locked to your input signal, then use the locked phase and frequency to ensure you take samples across as close to an integer number of periods as you can.

A somewhat simpler way of doing this is as follows: consider processing datasets like those you plotted above. If the data should have approximately one period of the signal, then you can get a finer estimate of the period by locating the sinusoid's zero crossings. Once you've estimated the period, then you can use that subset of the data for your RMS calculation. This might be sufficient, depending upon the accuracy that your application requires.
Thinking a bit outside the box: is the input signal a very pure sinusoid (with negligible noise), and is it always significantly oversampled? You could dispense with the RMS approach altogether and simply look for peak amplitude, then use this to calculate the sinusoid's RMS level.


Downsampling in Matlab simple example I'm trying to visualise downsampling in the frequency domain in matlab

I take a simple sinusoid, perform an fft and plot a two sided spectrum. Then I downsample the time domain signal (downsamplefactor D=2) and perform the same fft and two sided spectrum plot.

What i would expect to see: (This uses M, i use D as downsample factor)

(from:
Frequency Representation of Downsampled Signal)

So i'm expecting my frequencies to get stretched by a factor D and the amplitude should get scaled by 1/D. The stretching works as expected, however my amplitude doesnt scale. What am I missing/not seeing? Matlab code:

% normal signal -> fft -> plot
Fs = 100; %sampling frequency
T = 1/Fs; %sampling time
L = 500;  %signal length
f = 5;    %sine frequency

dF = Fs/L;      %frequency bin step
t = (0:L-1)*T;  %time vector
y = 2*sin(f*2*pi*t);    %sine

Y = fft(y);     %fft
P = abs(Y/L);   %extract magnitude plot

P2 = fftshift(P); %rearranges P to 2sided spectrum:negfreq,0component,posfreq

F2 = dF*(-L/2:1:L/2-1); %for when L is even
%F2 = dF*((-L+1)/2:1:(L+1)/2-1); %for when L is odd
F2rad = F2*pi/(Fs/2); % *pi /Nyquist to convert to [-pi pi]

figure, hold on
%plot(F2,P2) %frequency plot
plot(F2rad,P2) %radian plot

% downsampled signal -> fft -> plot

D = 2;
yd = y(1:D:L); %can also use downsample command
yd = downsample(y,D); %equal to above

Ld = length(yd); %downsampled length

Yd = fft(yd);

Pd = abs(Yd/Ld); %0component,posfreq,nyquist,negfreq
Pd2 = fftshift(Pd); %rearranges P to nyquist,negfreq,0component,posfreq

dFd = (Fs/D)/Ld; %frequency bin step

Fd2 = dFd*(-Ld/2:1:Ld/2-1); %for when L is even
%Fd2 = dF*((-L+1)/2:1:(L+1)/2-1); %for when L is odd
Fd2rad = Fd2*pi/(Fs/(2*D)); % *pi /Nyquist to convert to [-pi pi]

%plot(Fd2,Pd2) %frequency plot
plot(Fd2rad,Pd2) %radian plot

You are removing the scaling effect yourself when you perform the following operations:

Pd = abs(Yd/Ld)
P = abs(Y/L)

If you don't divide by the length, then you'll see that one has half the amplitude of the other.


What would be the variance for complex number? When $x$ is a zero mean random variable then

$$\sum_{n=1}^N x_n x_n^T = N \sigma^2_x\,\text,$$

where the variance is $\sigma^2_x$.

I'm considering Complex Normal Distributions where the real and imaginary part are uncorrelated.

https://en.wikipedia.org/wiki/Complex_normal_distribution explains about the form of the distribution for complex normal distribution.

I have a confusion because the denominator for the real case in the distribution has a 2 but for the complex that is not there.

Does this mean that if $x$ is a complex valued random variable, then the variance becomes half i.e.,
$$\sum_{n=1}^N x_n x_n^H = N \sigma^2_x/2\,\text,$$
where the variance is $\sigma^2_x/2$ because the variance gets equally distributed in the real and imaginary component?

I have this doubt because when implementing, if I need to generate a complex noise of variance 1, I would be doing (in Matlab)

noise = sqrt(1/2) * (randn(N,1) + 1j*randn(N,1))

Since each component (real and imaginary) needs to have variance 1/2, such that their sum becomes 1.

So, the variance $\sigma^2$ is half mathematically. Is my understanding correct?

UPDATE based on valuable information provided under the comments : I am considering at the circularly symmetric complex normal distribution, where real and imaginary part are completely uncorrelated.

I will focus on the reason of the factor $1/2$ and leave aside the estimation things.

The exact understanding should be : if a scalar Gaussian random variable (rv) is circular symmetric, its real and imaginary parts must be uncorrelated (this is equivalent to independence if they are assumed jointly Gaussian) and identically distributed with zero mean. Thus, your Matlab code is correct for a rv $\sim \mathcal{CN}(0,1)$.

The story behind is that a complex random variable (rv) is simply a vector of two real random variables. A vector of $n$ complex rv is indeed a vector of $2n$ real rv.

You are talking about the case $n=1$, 1 complex Gaussian rv $Z = Z_r + jZ_i$ or a vector of 2 real Gaussian rv $[Z_r,Z_i]^T$. As with real Gaussian rv which is described by its variance, a real vector $[Z_r,Z_i]^T$ must be described by its covariance matrix
$$\mathbb{E}\left\lbrace [Z_r,Z_i]^T \times [Z_r,Z_i] \right\rbrace = \mathbb{E}\left\lbrace\begin{bmatrix}
        Z_r^2 & Z_rZ_i \\
        Z_iZ_r & Z_i^2 \\
        \end{bmatrix}\right\rbrace$$

Take a look at the variance $\mathbb{E}\left\lbrace ZZ^H \right\rbrace = \mathbb{E}\left\lbrace Z_r^2+Z_i^2 \right\rbrace$ and pseudo-variance $\mathbb{E}\left\lbrace ZZ^T \right\rbrace = \mathbb{E}\left\lbrace Z_r^2-Z_i^2+j2Z_rZ_i \right\rbrace$, the covariance matrix of real vector (or complex scalar) is fully described by the variance and pseudo-variance. Threfore, you need both variance and pseudo-variance to characterize a complex Gaussian rv (with the prior condition that its real and imaginary parts are jointly Gaussian).

Now we use the circular symmetry property : $e^{j\phi}Z$ has the same probability distribution as $Z$ for all real $\phi$. This leads to $\mathbb{E}\left\lbrace e^{j\phi}Z(e^{j\phi}Z)^T \right\rbrace = \mathbb{E}\left\lbrace e^{j2\phi}ZZ^T \right\rbrace = \mathbb{E}\left\lbrace ZZ^T \right\rbrace$ for all $\phi$ thus $E[ZZ^T] = 0$ and variance $\mathbb{E}\left\lbrace ZZ^H \right\rbrace$ is sufficient statistic for $Z$. Note that $\mathbb{E}\left\lbrace ZZ^T \right\rbrace = \mathbb{E}\left\lbrace Z_r^2-Z_i^2+j2Z_rZ_i \right\rbrace \implies \mathbb{E}\left\lbrace Z_r^2 \right\rbrace = \mathbb{E}\left\lbrace Z_i^2 \right\rbrace$ and $\mathbb{E}\left\lbrace Z_rZ_i \right\rbrace=0$, the real and imaginary parts are uncorrelated then independent (because they are Gaussian), with the same variance, this is the reason of the factor $1/2$.

To sum up, your code is correct because you are estimating circular symmetric complex Gaussian rv. The jointly Gaussian assumption between real and imaginary parts must be used. If this is not about circular symmetric rv (or real random vector with two elements), you must calculate the pseudo-variance.

For more details and to understand the formula of the wikipedia article, you can read this article Circular Symmetric Gaussian R.Gallager.


Detection of sine signals with random amplitudes Problem statement:

I am designing a NP detector for the following detection problem:

$\mathcal H_{0}: x[n] = A_0\cos(2\pi f_0n) + w[n]$
$\mathcal H_{1}: x[n] = A_1\cos(2\pi f_0n) + w[n]$

where:

$f_{0}: \text{known parameter}$
$w[n] \sim \mathcal N\left(0, \sigma^2\right), \text {WGN}$
$A_{0}\ \ \sim \mathcal N\left(a_0, \sigma_0^2\right)$
$A_{1}\ \ \sim \mathcal N\left(a_{1},\sigma_1^2\right)$

What didn't work for me:

o I have tried to  directly apply Neyman-Pearson theorem in the way it's done in "Fundamentals of Statistical Signal Processing, Volume II, Detection Theory" by Steven M. Kay. But I have failed to derive test statistics for optimal detector.
o I have tried to use the amplitude of the spectral component for $f_{0}$ as a test statistics and estimate threshold of making decision from given $P_{FA}$ (probability of false alarm). It worked fine for me, but I'm still concerned that this detector is not optimal.

Question:

Is there some way to design optimal NP detector for such problem?
Any ideas, advices, links will be much appreciated.

Since both the null hypothesis $\mathcal H_0$ and the alternative hypothesis $\mathcal H_1$  are signal + noise, this is not only the detection of a random signal in WGN, but also a discrimination problem (i.e. "classification") of one signal from the other both embedded in WGN. The Neyman-Pearson detector decides $\mathcal H_1$ if:

$$\frac{p(\mathbf x; \mathcal H_1)}{p(\mathbf x; \mathcal H_0)}> \gamma$$

Which is:

$$
\frac{\displaystyle\frac{1}{\left[2\pi\left(\sigma_1^2 + \sigma^2\right)\right]^\frac N2}\exp\left[-\frac{1}{2\left(\sigma_1^2 + \sigma^2\right)}\sum_{n=0}^{N-1}\left(x[n]-a_1\cos\left(2\pi f_0 n\right)\right)^2\right]}{\displaystyle\frac{1}{\left[2\pi\left(\sigma_0^2 + \sigma^2\right)\right]^\frac N2}\exp\left[-\frac{1}{2\left(\sigma_0^2 + \sigma^2\right)}\sum_{n=0}^{N-1}\left(x[n]-a_0\cos\left(2\pi f_0 n\right)\right)^2\right]}>\gamma
$$

Taking the logarithm on the LHS you get:

\begin{align}
T(\mathbf x)&=\frac N2 \ln\left(\frac{\sigma_0^2 + \sigma^2}{\sigma_1^2 + \sigma^2}\right) +\frac{\sum_{n=0}^{N-1}\left(x[n]-a_0\cos\left(2\pi f_0 n\right)\right)^2}{2\left(\sigma_0^2 + \sigma^2\right)}-\frac{\sum_{n=0}^{N-1}\left(x[n]-a_1\cos\left(2\pi f_0 n\right)\right)^2}{2\left(\sigma_1^2 + \sigma^2\right)}\\
&=\frac N2 \ln\left(\frac{\sigma_0^2 + \sigma^2}{\sigma_1^2 + \sigma^2}\right) +\left[\frac{\left(\sigma_1^2-\sigma_0^2\right)}{2\left(\sigma_0^2+\sigma^2\right)\left(\sigma_1^2+\sigma^2\right)}\right]\cdot\sum_{n=0}^{N-1}x^2[n] \\&\quad+ \left[\frac{\left(a_1-a_0\right)\sigma^2-a_0\sigma_1^2+a_1\sigma_0^2}{\left(\sigma_0^2+\sigma^2\right)\left(\sigma_1^2+\sigma^2\right)}\right]\cdot\sum_{n=0}^{N-1}x[n]\cos\left(2\pi f_0 n\right)\\&\quad+\left[\frac{a_0^2}{2\left(\sigma_0^2 + \sigma^2\right)}-\frac{a_1^2}{2\left(\sigma_1^2 +\sigma^2\right)}\right]\cdot\sum_{n=0}^{N-1} \cos^2\left(2\pi f_0 n\right)
\end{align}

The inequality becomes $T(\mathbf x) > \ln\gamma$. Keeping data-dependent terms on the LHS and defining $\gamma'$ as:

$$
\gamma'=\ln\gamma-\frac N2 \ln\left(\frac{\sigma_0^2 + \sigma^2}{\sigma_1^2 + \sigma^2}\right)-\left[\frac{a_0^2}{2\left(\sigma_0^2 + \sigma^2\right)}-\frac{a_1^2}{2\left(\sigma_1^2 +\sigma^2\right)}\right]\cdot\sum_{n=0}^{N-1}{\cos^2\left(2\pi f_0 n\right)}
$$

We get:

$$
\underbrace{\overbrace{\frac{\left(\sigma_1^2-\sigma_0^2\right)}{2\left(\sigma_0^2+\sigma^2\right)\left(\sigma_1^2+\sigma^2\right)}\sum_{n=0}^{N-1}x^2[n]}^{\rm Quadratic} + \overbrace{\frac{\left(a_1-a_0\right)\sigma^2-a_0\sigma_1^2+a_1\sigma_0^2}{\left(\sigma_0^2+\sigma^2\right)\left(\sigma_1^2+\sigma^2\right)}\sum_{n=0}^{N-1}{x[n]\cos\left(2\pi f_0 n\right)}}^{\rm Linear}}_{T'(\mathbf x)}> \gamma'
$$

The test statistic $T'(\mathbf x)$ contains both a linear term and a quadratic term, the detection should then be based on both the variance characterized by the quadratic term and the "mean" characterized by the linear term. More like a composite energy detector. That's my two cents for now.

By the way, this problem looks like a special case the general Gaussian detection problem where the signal has a deterministic part and a random part. And this under $\mathcal H_1$ is modeled as the Bayesian general linear model as follows:

$$
\mathbf{ x =  H\theta + w}\tag{1}
$$

Where $\mathbf H$ is a $n\times p$ deterministic matrix , $\mathbf \theta$ is a $p\times 1$ random vector with PDF $\mathcal N\left(\mu_\theta, \sigma_\theta^2\right)$ and independent of  $\mathbf w$.

If the problem was to detect the presence of a signal (any) versus noise, then one could transform the initial definition to for instance $\mathcal H_0$ noise only versus $\mathcal H_1$ signal 1 + signal 2. And this could be brought to the case in $(1)$.

UPDATE:

Note two special cases:

When the random amplitudes have equal variances but different means, that is if $\sigma_0^2 = \sigma_1^2$ but $a_0 \neq a_1$, the detection problem reduces to the one of simply tracking the "mean" as the quadratic term (i.e energy part) cancels out and the test statistic becomes

$$
T'(\mathbf x)=\left(\frac{a_1-a_0}{\sigma_0^2+\sigma^2}\right)\sum_{n=0}^{N-1}{x[n]\cos\left(2\pi f_0 n\right)}
$$

And for the case where $a_0 = a_1$ but $\sigma_0^2 \neq \sigma_1^2$ you get both the quadratic and linear term with the same coefficient by a factor of $-2a_0$.

$$
T'(\mathbf x)=\frac{\left(\sigma_1^2-\sigma_0^2\right)}{2\left(\sigma_0^2+\sigma^2\right)\left(\sigma_1^2+\sigma^2\right)}\sum_{n=0}^{N-1}x^2[n] + \frac{\left(\sigma_0^2-\sigma_1^2\right)a_0}{\left(\sigma_0^2+\sigma^2\right)\left(\sigma_1^2+\sigma^2\right)}\sum_{n=0}^{N-1}{x[n]\cos\left(2\pi f_0 n\right)}
$$

So, the random amplitudes $A_0$ and $A_1$ play part in the detection not only through their respective means $a_0$ and $a_1$ but also through their variances $\sigma_0^2$ and $\sigma_1^2$.


Impulse Response to LTI I am new to DSP, and I am self-studying using mostly Proakis.  I have a question.  There are some examples in the text where you will be given the impulse response of an LTI system, and then asked to solve something/prove something, so forth.  My confusion involves the following: How do I know that there exists an LTI system with the given impulse response?  For instance, if a problem is posed like one of the examples in Chapter 2, "given a relaxed LTI system with impulse response:

$$h[n] = a^n u[n]$$

where $-1 < a < 1$," (and $u[n]$ is the discrete-time unit step function) how do I know that this impulse response actually comes from an LTI system?  Sorry if this question is simple.  I'm probably not seeing something.

There is a general question and a specific question.

To answer the specific question, you can guess at an LTI and check if it satisfies your spec:

$$ h[n] = a^n u[n]  \qquad \text{where} \ |a| < 1 $$

Here is my guess:

$$ y[n] = a \cdot y[n-1] + x[n]  \qquad \text{where} \ y[-1]=0 $$

you can set $x[n] = \delta[n]$ and you will find that $y[n] = h[n]$ for all integer $n$ and you are done with the specific question.

To see if a conceptual discrete-time LTI can be made to satisfy a general impulse response spec, $h[n]$, there are a few requisite properties to satisfy:

Is the system causal?  This means the output $y[n]$ does not react to future samples of the input $x[n]$.  And if it's both LTI and causal, that means the impulse response does not react to the input impulse before that input impulse takes on non-zero values.  And that means $$ h[n] = 0 \qquad \forall \ n<0$$  You can see that this is satisfied by your example spec.
Is the LTI system FIR or IIR?  If it's FIR, that means that the non-zero length of the impulse response, $h[n]$ is finite in length (let's call that FIR length "$M$") and that means $$ h[n] = 0 \qquad \forall \ n<0, \ n \ge M$$ Now it turns out that your example is not FIR (which means it is IIR, and Infinite Impulse Response LTI system), but if it were FIR, then your general system can always be realized (or "actualized") as $$ y[n]=\sum\limits_{m=0}^{M-1} h[m] \, x[n-m] $$  That can always be done for an FIR.  It might not be the most efficient way to make an FIR for extremely large $M$, but it is a method that, if you can afford it, will always work.
If it's an IIR (which means the condition in paragraph 2 above is not satisfied), then you have some restrictions about what can be realized and what cannot.  And you have some other conditions to nail down.  One is "BIBO stability" which means, for any bounded input $x[n]$, the output $y[n]$ is also bounded.  Now, I don't know all of the implications, but that does mean that systems that blow up, like $$ h[n] = a^n u[n]  \qquad \text{where} \ +1 < |a| $$ can be realized (same input/output equation as in 2), but they will blow up.  I don't think you want that.  BIBO stability also means that impulse responses like $$ h[n] = \tfrac{1}{n+1} u[n] $$ will also blow up because $\tfrac{1}{n+1}$ does not converge fast enough.  If, say, $x[n]=u[n]$ (which is bounded), the output $y[n]$ will blow up.
So it turns out that BIBO stable LTI will need to either be FIR (they're always stable from an input/output POV) or, if the LTI is IIR, then the impulse response must decay exponentially.  That means, for large $n$, $$  \Big|h[n] \Big| \le A \ a^n \ u[n] \qquad \text{for some} \ |a| < 1  $$  If that's true, you can conceptually make an IIR filter that approximates any $h[n]$ that satisfies the causality condition (paragraph 1) and the BIBO stability conditions (paragraph 3).  The way we learn to do that is by taking a DSP course after learning LTI system theory.


Is Differentiation as a system, is an invertible system? is the following system invertible?

as I understand it, invertible means finding an inverse function which should return back the original input from an output of the given system.

if so I understand correct. Then one could say, yes the system is invertible and the inverse function is

but for example, If our original input is x + 3 then its derivative is 1,  and if I try to use the inverse function it returns x + c.   is the nature of not knowing what exactly c equals makes the system non-invertible; meaning the original input cannot be restored. Does that make the system non-invertible

Yes, your understanding is correct. The differentiation is not an invertible system. And, the reason is precisely, because of the possible unknown constant which gets nulled by differentiation.

Basically, the differentiation of $f(t)$ and $(f(t) + c)$ is indistinguishable $f^{'}(t)$, which makes the system, many-to-one, and hence non-invertible.


Difference between cross correlation definitions (signal processing / random processes) I'm trying to understand the similarities/differences between the following definitions for cross correlation)
Signal Processing: (or it's discrete equivalent)
$$R_{XY}(\tau) = \int_{-\infty}^{+\infty} x^*(t) y(t+ \tau) dt$$
Random Processeses$$R_{XY}(t_1, t_2) = E\begin{Bmatrix}
X_{t_1} Y^*_{t_2}
\end{Bmatrix}$$
Are these related?

If the process is stationary it is ergodic, then time averages are equal to probabilistic averages.


Frequency response of a long signal I am trying to get the frequency response of a system, so, in theory I would need to do:
H = fft(y)./fft(x);

I have measured the system with a long stepped sine log sweep four times as I am planning to perform a spectral averaging to improve my SNR. This gives me a total of almost 5.000.000 samples per file. I have calculated the STFT of smaller chunks and averaged them as it looked like a too long signal to do an FFT. I used the cross spectrum and auto spectrum for that (I did the same procedure 4 times, one per audio file/repetition).
However I think that's not the correct way to go as I am averaging a non stationary signal. Should I just do an FFT of the whole signal? Isn't it too long? Do you have any idea of how I should proceed?
Thanks a lot in advance!

Assuming that the system response is smooth, without narrow peaks or narrow valeys. And that the spectrum of the input signal does not vary so much from one chunk to another.
I think this approach is better than doing a single FFT, as you know the FFT of a random signal will have lots of randomness. Averaging may also remove effects due to quantization errors.
Since the FFT bins are equally distributed in the spectrum, each bin corresponds to a very narrow band of your spectrum, and the energy in this band is highly dependent on the signal you used, or noise, if the FFT is big. In order to have a smooth curve you can compute the energy in a wider band by applying a moving average, for instance.


Does $\cos(bt)\cdot u(t)$ have a Fourier Transform? If it does,
$$\int_{-\infty}^{\infty} \cos(bt)\,u(t)e^{-j\omega t} dt = \int_{0}^{\infty} \cos(bt)\,e^{-j\omega t} dt = \int_{0}^{\infty} \frac{e^{jbt} + e^{-jbt}}{2}\,e^{-j\omega t} dt$$
Then how do we carry on from this point?
$$ \frac{1}{2} \int_{0}^{\infty} \left[e^{j(b-\omega )t} + e^{-j(b+\omega )t}\right]dt$$
Thank you.

The integral doesn't converge in the conventional sense, so you can't solve it with standard methods. Assuming that you know (or can look up) the Fourier transform of the unit step function $u(t)$, it is straightforward to compute the Fourier transform of $\cos(\omega_0t)\,u(t)$ using the modulation property:
$$\mathcal{F}\big\{u(t)\big\}=U(\omega)=\pi\delta(\omega)+\frac{1}{j\omega}\tag{1}$$
$$\begin{align}\mathcal{F}\big\{\cos(\omega_0t)\,u(t)\big\}&=\frac12\big[U(\omega-\omega_0)+U(\omega+\omega_0)\big]\\&=\frac{\pi}{2}\big[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)\big]+\frac{1}{j}\frac{\omega}{\omega^2-\omega_0^2}\tag{2}\end{align}$$


Cross-correlation peaks in acoustic multi-path conditions Consider a speaker and a microphone are placed in a multi-path environment and the audio signal to be transmitted is LFM(chirp) signal, the time of arrival of the audio signal is computed by cross-correlation at the receiver side. My question is:

If there are severe reflection, reverberation phenomenons but the direct-path propagation acoustic signal is also received. When doing cross-correlation, will the first correlation peak (which is supposed to belong to direct path propagation signal) shift due to reflections & reverberations?

thank you very much

Depending on the circumstances, you can also have other than air borne propagation like through a concrete floor and the support structures of your projector and sensors.  the fastest path isn’t necessarily going to be the direct air path.


I have a question on SIFT I am reading a paper on SIFT by Lowe.

In it, gradient magnitudes and orientations are calculated using these formulae:

$$$$
I don't know why m(x,y) is calculated above. I heard that gradient magnitude is usually calculated like the following.

$$\sqrt{f_x ^2 + f_y^2}
$$

Horizontal and vertical gradients are computed by taking neighbor pixel differences:

$$g_{x}=L(x+1,y)-L(x-1,y)\\g_{y}=L(x,y+1)-L(x,y-1)$$

Gradient magnitude is computed the same way as in your formula:

$$m(x,y)=\sqrt{g_x^2+g_y^{2}}$$

Replacing $g_{x}$ and $g_{y}$ with above will give you the original formula.

Gradients are usually computed by forward or central difference formula, but I think in this case the above is still accurate enough since the image $L$ is smoothed.


How to get fundamental frequency of a signal using autocorrelation? I am trying to get the fundamental frequency of a signal that only has a single pitch. I coded out the autocorrelation function using FFT and already got the autocorrelation result. Unfortunately, I don't know how to get the fundamental frequency from the autocorrelation result. Can someone help me? My JAVA code for getting the fundamental frequency using autocorrelation is below:

public double getPitch(double[] buffer, int firstSample, int lastSample, double sampleRate)
{
    int lengthOfFFTWindow=lastSample-firstSample;
    double[] input_buffer=new double[lengthOfFFTWindow];
    DoubleFFT_1D fft = new DoubleFFT_1D(lengthOfFFTWindow);
    double[] autocorrelation_values=new double[lengthOfFFTWindow];
    double[] fftData = new double[lengthOfFFTWindow * 2];
    double max=-1;
    double max_i=-1;
    //FFT on each sample in each window
    for (int i = 0; i < lengthOfFFTWindow; i++) {
        // copying audio data to the fft data buffer, imaginary part is 0
        fftData[2 * i] = buffer[i+firstSample];
        fftData[2 * i + 1] = 0;
    }
    fft.complexForward(fftData);
    for (int i = 0,j=0; i < fftData.length; i += 2,j++) {
        // complex numbers -> vectors, so we compute the length of the vector, which is sqrt(realpart^2+imaginarypart^2)
        autocorrelation_values[j] = Math.sqrt((fftData[i] * fftData[i]) + (fftData[i + 1] * fftData[i + 1]));
    }
    fft.complexInverse(fftData, false);
    for(int i=0;i<autocorrelation_values.length;i++)
    {
        if(max<autocorrelation_values[i])
        {
            max=autocorrelation_values[i];
            max_i=i;
        }
    }
    return max/2;
}

Is it correct to return the maximum autocorrelation value divide by 2 as the fundamental frequency? I keep getting wrong answers when I do that.

The value of the peak of the autocorrelation function (maximum) is not really important - though it can be used as a basic feature for voiced/unvoiced classification or more generally for "pitchedness".

What matters is the index at which the function takes its maximum (arg-maximum). This is very intuitive to understand: if the signal has a period of $\tau$ samples, the original signal and the signal shifted by $\tau$ samples should be very similar, so the autocorrelation function will have a peak when evaluated at $\tau$.

Thus, the interesting variable is max_i in your code. One thing I notice is that you search through the whole autocorrelation_values vector. But the autocorrelation function is maximum at 0, so your maximum will always be 0! I advise you to think about a meaningful range for the fundamental frequency you want to detect, and restrict the search in this range to get a meaningful result.

I hadn't checked the rest of your code in details, I don't think it is the purpose of this site to dive into coding/implementation problems.


What is the relation between windowing and hopping in audio DSP While using short-time Fourier transform (STFT), it is common to use window function, usually tapered.

Will this not alter the original signal $x[n]$ and the STFT may no more be the true STFT of $x[n]$?
Also, what is the concept for sliding the frames not by $100\%$ of its width but a fraction of that? i.e, why should frames hop in small steps (for $75\%$ overlap, frames hop by $25\%$). Is there a conceptual explanation with as little math as possible?

What many people do not seen to realize is that using a finite length transform (such as one iteration of an STFT FFT) on a longer signal already alters the longer signal by windowing (e.g. ignoring the portions of the signal that do not fit in the current finite length FFT).  Not using a tapered window means that one is using a rectangular window.  Using a rectangular window may alter the signal spectrum even more than using a tapered window (from, say the infinite length FT of an infinitely long stationary signal, especially if the signal is not perfectly integer periodic in the STFT FFT aperture width).

The FFT spectrum is altered the least near the center of the window (where a tapered window has a gain of unity or a maximum gain, and non-tapered signal is farthest from both potential edge transients).  Therefore the closer together the window centers of individual STFT FFTs (the shorter the hops), the less the alteration or distortion of the total STFT spectrum.

Lots of overlap also provides greater time localization or locality (e.g. the STFT won't miss or slice up a short event that is right between two windows).
 However the greater the overlap, the greater the computation required.  So a compromise overlap of 50% or 75% window overlap is common.


Chirp after T seconds Suppose we have this signal from this wiki page

$a(t)=\sin(\varphi_o + f_0t + (f_1-f_0)t^2/T)$

The article tells that T is a period of chirp modulation in which frequency changes. So my question is how to create signal with constant $f_1$ frequency after T?

how to create signal with constant $f_1$ frequency after T?

Find the instantaneous phase at the end of the $T$ interval and keep adding it to the phase that makes the $\sin()$ "go round".

A plain simple oscillator at some frequency $f$ has its phase increasing at a constant rate.

Here is an example (in GNU Octave but easily transferable to other platforms as well):

Fs = 44100; % Sampling frequency in Hz
T = 1; % Total duration of the signal in seconds
f = 120; % Frequency of the oscillator
t = 0:(1./Fs):(T - (1./Fs)); % Time vector
p = 2.0 .* pi .* t; %Phase vector
y = sin(f.*p);

Now y contains our 120 Hz sinusoid sampled at 44.1kHz.

Notice here that phase (p) is a "straight line" that grows at a constant rate of $\frac{2 \pi f}{Fs}$.

When the oscillator is "chirping", the rate of change of the phase is variable. For example:

Fs = 44100; % Sampling frequency in Hz
T = 1; % Total duration of the signal in seconds
f0 = 1; % Start chirp at
f1 = 120; % end chirp at (Better keep f1>f0)
t = 0:(1./Fs):(T - (1./Fs)); % Time vector
p = 2.0 .* pi .* t; %Phase vector
y = sin(f0.*p + 2.*pi.*(((f1-f0)/(2.*T)).*t.^2));

Now, this looks like:

And the rate of change of the phase of the $\sin$ is:

Now, from the first example, the phase "step" is f.*p(2)-f.*p(1) == 0.017097 and from the second example, after doing:

s = diff(f0.*p + 2.*pi.*(((f1-f0)/(2.*T)).*t.^2)); % Get the first derivative of phase to find the **rate of change**.

The "last" rate of change is s(end) == 0.017097.

So, to have the oscillator stay at the chirp's concluding frequency, keep accumulating the phase at the constant rate you are after from the last known value of the chirp's phase (to avoid phase jumps). For example:

r = f0.*p + 2.*pi.*(((f1-f0)/(2.*T)).*t.^2);  % The phase as above.
z = r(end) + cumsum(ones(1,T.*Fs).*0.017097); % A running sum of duration T.*Fs (samples).

Now, if you feed the combined result of r, z into the oscillator y with something like y = sin([r,z]); it would look like:

...with a rate of change of:

Hope this helps.


Detection of regions of voice in music I want to detect the regions of voice in a song. Currently, I am using the fact that the vocals in the song are recorded as center panned. So, I am removing all center panned audio from the song. However, this removes other instruments such as bass also. So, it sometimes gives false indication of voice being present.

I want a way ( better? than this naive approach?) to be able to detect portions of the song where a singer is singing. There is no specific type of music, I'd like an implementation for any type of song. A sample might be http://www.youtube.com/watch?v=w_Rut4qm33g but I up for any song for which this can be done, even if it has some specific properties.

My current implementation is -

[y, fs] = wavread('song.wav');

subplot(2, 1, 1);
plot(y)

left = y(:,1);
right = y(:,2);
yOut(:,1) = (left - right)/2;
yOut(:,2) = (left - right)/2;
wavwrite(yOut, fs, 'tmp2.wav');

fftL = fft(left);
fftR = fft(right);

N = size(fftL)
N = N(1);

M = size(fftR)
M = M(1);

N = min(N, M);

for i = 1:N
  dif = fftL(i,1) / fftR(i,1);
  dif = abs(dif);
  if (dif > 0.7 && dif < 1.5)
    fftL(i,1) = 0;
    fftR(i,1) = 0;
  end;

end;

leftOut = ifft(fftL);
rightOut = ifft(fftR);
yOut(:,1) = leftOut;
yOut(:,2) = rightOut;

subplot(2, 1, 2);
plot(yOut)

N = size(yOut);
N = N(1);
for i = 6:N
    yOut(i) = 1/5* (yOut(i-1) + yOut(i-2) + yOut(i-3) + yOut(i-4) + yOut(i-5));
end
wavwrite(yOut, fs, 'tmp.wav');

Warning: this is still an active research topic so there is no "canned" solution.

As you have stated, middle channel suppression is not a very robust method. It can help in getting a very rough "Karaoke" version but you can't infer anything reliable from the resulting audio.

A simple approach would consist in extracting a bunch of audio features on each signal frame, and use machine learning techniques (for example support vector machines) to discriminate the "voice present"/"no voice" classes. This works relatively well (reported recognition rates in the literature are in the 70-80% range), but the downside is that you have to manually annotate material to get training data. You might try to use karaoke files (with precise timing of the lyrics) and align them with the audio to get a kind of automatic transcription, but remember that your system will only be as good as your training data. Your training database will have to cover all the possible genres of music on which your classification system will be unleashed. See Rocamora's Master's thesis or Ramona et al's paper. There are dozens of other papers presenting variations of this scheme - with different feature sets, post-processing/temporal-smoothing operations, etc. Follow the references in the papers.

If you can't afford to spend countless hours in Transcriber or Wavesurfer annotating data, there are still unsupervised methods that would reasonably work. I recommend you to look into lead melody separation methods such as:

"REpeating Pattern Extraction Technique (REPET): A
Simple Method for Music/Voice Separation" (Rafii / Pardo). This is a very simple method which works very well on "commercial" pop music (with a very stereotypical structure and modern computer-based production).
"Source/Filter Model for Main Melody Extraction From Polyphonic Audio Signals" (Durrieu et al). Complex algorithm, but the code is available.
If you want something really simple, just extract and track the most salient pitch and assume it is the lead melody. This works to some extent in recordings in which the singer voice is well isolated from the mix (other instruments EQ'ed so that they do not overlap much with it).

You don't need to bother here with resynthesis - just use the energy of the extracted spectrum / f0 track.

A last angle of attack is to use a multi-f0 estimator or a sinusoids tracker, and check if one of the tracks has frequency fluctuations typical of vibrato or glissando. In commercial pop music, most instruments have a very stable pitch contour (keyboards, fretted bass and guitars), and the voice is the only instrument with a very wobbly pitch contour - unless it has been doctored by autotune.


Compute the time derivative of a noisy digital signal? The issue is that my signal is very noisy. I need extract its time derivative as accurate as possible.  P.S. I do not have any prior knowledge on the signal (black box).
On forums some suggested Savitzky-Golay filter.
Any idea please? If so, is there any c++ library for the purpose?
In fact, for my application I need to compute optical-flow like information for control purpose. I compute an estimate using image information. Then I need to compute the time derivative of this estimate.
4th order Savitzky Golay filter introduces delay, yet I need the output in real-time (real time control). For info:

The signal is regularly sampled;
The noise is not defined but bounded;
The output needs to be real-time: delay-minimal;
The signal is black box: I only get a measure each iteration.

If computational efficiency is important to you, and your signal's bandwidth is not to large relative to the $f_s$ sample rate, three simple causal differentiators that tend to attenuate high-frequency noise are described at DSPRelated.

A plot of the results of that article is below.

Applying these three filters to your example signal (or an attempt to mimic it) yields the following plots.

R Code Below

# Q27420
h_ref <- c(-1/16, 0, 1, 0, -1, 0, 1/16)
d2 <- c(0.5,0,-0.5)
d3 <- c(-3/16, 31/32,0, -31/32, 3/16)

H_ref <- abs(fft(c(h_ref,rep(0,1000))))
f2 <- abs(fft(c(d2,rep(0,1000))))
f3 <- abs(fft(c(d3,rep(0,1000))))

w1 <- seq(0,length(H_ref)-1,1)/(length(H_ref)-1)
w2 <- seq(0,length(f2)-1,1)/(length(f2)-1)
w3 <- seq(0,length(f3)-1,1)/(length(f3)-1)

n1 <- length(w1)/2
n2 <- length(w2)/2
n3 <- length(w3)/2

plot(w3[1:n3],f3[1:n3]/1.2, col="blue", lwd=5,type="l", lty=2)
lines(w1[1:n1],H_ref[1:n1]/1.63,col="black", lwd=5)
lines(w2[1:n2],f2[1:n2], col="green", lwd=5)

legend(0.0, 1.5, c("Proposed","Reference", "Difference"),
       col=c("blue","black","green"), lty=c(2,1,1), lwd=c(5,5,5));

x <- c(seq(-1.6,-0.4,0.01), seq(-0.4,-0.2,0.001))
xn <- x + rnorm(length(x))/10

f <- seq(0,1,0.1)
m <- f
d4 <- signal::fir2(100,f,m)

y1 <- filter(xn,h_ref, circular=TRUE)
y2 <- filter(xn,d2, circular=TRUE)
y3 <- filter(xn,d3, circular=TRUE)
y4 <- filter(xn,d4, circular=TRUE)

par(mfrow=c(4,1))
plot(xn)
title("Original signal")
plot(y4,lwd=4,col="grey")
lines(y1,col="blue")
title("Deriviative calcualted using the reference technique")
plot(y4,lwd=4,col="grey")
lines(y2,col="black", lty=2, lwd=3)
title("Deriviative calcualted using the proposed technique")
plot(y4,lwd=4,col="grey")
lines(y3,col="green", lwd=5)
title("Deriviative calcualted using simple difference")


Spectral subtraction - Python Is there a good, documented open source implementation of spectral subtraction in Python?

I've searched a lot but I've only come across pieces of code where I can't really understand what's going on as they are undocumented. For example. Can someone either point me to some implementation or explain step by step what is happening in the above link or here for example.

Of the two Python code examples you found, I think the second one is clearer and likely to serve you better in your efforts to understand spectral subtraction and develop a Python script-file for your work.

The "noise_subtraction.py" script-file implements a basic power spectral subtraction method.  A good paper that describes this method is:  "Enhancement of speech corrupted by acoustic noise", Proceedings of the International IEEE Conference on Speech, Acoustics and Signal Processing, 208-211, 1979.  Since you are student at Imperial College, you should have no problem finding this paper.

A related resource that you'll want to find is the Speech Enhancement: Theory and Practice textbook.  It describes several spectral subtraction methods and includes functional Matlab script-files that demonstrate the methods.  In particular is a Matlab implementation of the method described in the Berouti et al. paper.

I know you are looking specifically for Python code, but if you review the material I've suggested and are able to obtain the basic power spectral subtraction Matlab code in the Loizou textbook you'll be in a position to use the "noise_subtraction.py" code you found.

As a final note, see this site for additional papers and Matlab code on speech related processing.

I hope this helps.


How to periodically estimate states of a LTI if the output is measured irregularly? How can I periodically estimate the states of a discrete linear time-invariant system in the form $$\dot{\vec{x}}=\textbf{A}\vec{x}+\textbf{B}\vec{u}$$
$$\vec{y}=\textbf{C}\vec{x}+\textbf{D}\vec{u} $$if the measurements of its output $y$ are performed in irregular intervals? (suppose the input can always be measured).

My initial approach was to design a Luenberger observer using estimates $\hat{\textbf{A}}$, $\hat{\textbf{B}}$, $\hat{\textbf{C}}$ and $\hat{\textbf{D}}$ of the abovementioned matrices, and then update it periodically every $T_s$ seconds according the following rule:

  If there has been a measurement of $y$ since the last update: $$\dot{\hat{x}}=\hat{\textbf{A}}\hat{x}+\hat{\textbf{B}}\hat{u}+\textbf{L}(y_{measured}-\hat{\textbf{C}}\hat{x})$$
  If not:
  $$\dot{x}=\hat{\textbf{A}}\hat{x}+\hat{\textbf{B}}\hat{u}$$

(I have omitted the superscript arrows for clarity)

I believe that there may be a better way to do this, since I'm updating the observer using an outdated measurement of $y$ (which is outdated by $T_s$ seconds in the worst case).

Thank's in advance.

My intuition on this is that you would actually perform a single instantaneous update in the moment that you receive a measurement, where this update is dependent upon the time since the last update.

The reason for this intuition goes as follows: Consider finely discretizing the system. Let's assume tht the matrices $A,B,C,D$ are known. In fact, let's just ignore the matrix $D$, since we can equivalently define observation $\tilde{y} = \mathbf{C}\vec{x} = \vec{y} - \mathbf{D}\vec{u}$. Finally, let's assume some distribution over the initial value $\vec{x}(0)$. In fact, call it gaussian for the niceties of it. This gives us a hidden markov model which can be solved by gaussian message passing on the resultant probabilistic graphical model. By "solve" I mean solving $E[\vec{x}[t]]$, where $\vec{x}[t]$ is the discretized version of the original $\vec{x}$. The factor graph (see paper Factor Graphs and the Sum-Product Algorithm) corresponding to the (hidden, since $\vec{x}$ is not observed) markov system with simply be a long chain of $\vec{x}$ variables and the deterministic relationship between them, every once in a while having an observation coming in. Realtime inference in this graph would consist of essentially a kalman filter, except that most of the time is just the "prediction" step. Only in the instances that there is an observation would you have an "update" step.

Basically, what I'm speculating is that you should actually do something like the following:

  If there is no measurement at time $t$, then simply carry on as
  $$\dot{\hat{x}} = \mathbf{A}\hat{x}+\mathbf{B}\vec{u}$$
  However, when an observation is received at time $t$, Update $\hat{x}$ as
  $$\hat{x}(t) = \hat{x}(t^-) + \mathbf{L}(\vec{y} - \mathbf{C}\vec{x}(t^-))$$

where $t^-$ is the simulated state up until right before the observation was received. This $\mathbf{L}$ would be designed by you, and should depend on the time since the last observation. Exactly the right way to do this escapes me so far unfortunately, but the finely discretized kalman filter analogy may shed some light on what is reasonable here.

Note that this is only considering obtaining an estimate for the hidden state $\vec{x}$. If the state estimate is then to be used for control purposes, one needs to take care about the fact that there are discontinuities in the derivative $\dot{\hat{x}}(t)$ for each $t$ when an observation is received. This can likely be mitigated by simply assigning to those values of $t$ the left or right limits of $\dot{\hat{x}}$ around that observation time.

If I obtain more information on this, I will update.


Correlated phase noise and delay 2 identical signals generated by the same PLL/VCO should have the same phase noise profile and their respective phase noises should be correlated.

Suppose one of them is delayed and I start sampling both simultaneously , if I now compare the phase noise of their respective samples taken at a certain time instance would the phase noise still be correlated? and if not how long would the delay have to be in order for their phase noise to be uncorrelated (if decorrelation is involved)?

My guess is that since phase noise has a trapezoidal shape and has a low pass filter frequency response, it is roughly the reciprocal of the loop BW that determines the length of the delay that would cause phase noise to become uncorrelated , the wider the loop BW i.e. the shorter the delay.

Does my reasoning make sense?

No, it has nothing to do with loop BW (as you stated that you are observing the final output, so it is the same signal at all offset frequencies;  as a very simplified explanation the loop bw indicates the demarcation in phase noise components between the VCO and the reference frequency).

If the output of the PLL/VCO is the same, meaning we split it coherently into two outputs and observe both outputs with no delay between the two, then they will be completely correlated at all offset frequencies. However, if there is a delay between the observations of the two signals, the components of the phase noise at very small frequency offsets from the carrier will still be correlated, while offsets much further away from the carrier will be decorrelated. As the delay between the two observations increases, the "frequency corner" where we start to decorrelate will go down in frequency as given by the inverse of the delay.  The corner frequency is the inverse of the delay (in radian frequency for delay in seconds), and the decorrelation for the phase noise component goes down as $1/\sqrt{f}$, where $f$ is the offset frequency from the carrier.

To determine the amount of correlation assuming the two waveforms have the same scaling (as long as there is no amplitude difference, which would apply if we simply delayed one of them as a delay does not change the amplitude), we can simply subtract the two and measure the resulting output power of the waveform (averaged detector).  The amount of cancellation indicates the degree of correlation; if the waveforms are completely correlated at all frequencies we would see complete cancellation, and if they were completely uncorrelated at all frequencies we would see a summation in power (+3 dB). In the continuous time domain such a delay and subtract operation has a high pass response, as suggested by the equivalent filter of that process shown in the graphic below, where the frequency shown represents offset frequency from the carrier:

The frequency response to a noise process of the delay and subtraction operation, is given by:

$$H(\omega)= \sqrt{1 - e^{-s\tau}}$$

(I believe the square root is included as I have shown since this is the result of summing random variables; given that independent random variables would add in variance, while the response is a magnitude quantity or the resulting standard deviation).

This shows the decorrelation intuitively, as components fully correlated would subtract, creating a cancellation, while the components fully uncorrelated, being random noise, would add in power (thus +3 dB higher). In terms of phase noise specifically this is clear when you consider a "slice" of phase noise at a a 3 Hz offset from the carrier (for example, as shown in the graphic below), represents a noise term that is sinsusoidally varying at a 3 Hz rate (at that given power level within 1 Hz of BW at that offset as given by the typical dBc/Hz scale of the vertical axis in phase noise spectrum plots). In comparison, consider a phase noise component from an offset of 0.03 Hz; that represents phase varying much more slowly (it takes 33.3 seconds to go through one cycle); surely you can see how even with a 1 second delay would still be quite correlated to this term, with a significant cancellation ratio (from such a subtraction).

The graphic below is specific to a slide I have on explaining such considerations with regard to carrier tracking loops, however directly applies the the explanation above, showing the cancellation of the lower frequency phase noise offsets.

Relationship to Cycle-to-Cycle Jitter
It is worth noting at this point how this explanation unifies the relationship between phase noise and cycle-to-cycle jitter; as the measurement of cycle-to-cycle jitter involves subtracting the rising edge of a clock signal from the previous rising edge; which is the block diagram in the first figure above where $\tau$ is the time duration of one clock cycle. Thus if you multiply the frequency response of that filter with the phase noise spectrum (to filter it), and then integrate to get the total rms phase noise, the rms jitter is found by converting phase to time based on the clock frequency used.


How do I apply a Chebishev filter? I read a paper about a brain-computer interface. In this paper the authors reported "each signal has been filtered with an 8-order band-pass Chebishev Type I filter which cut-off frequencies are 0.1 and 10 Hz and has been decimated according to the high cut-off frequency". I tried to design this filter with scipy:

import scipy.signal as signal
signal.cheby1(8,0.05,[0.1,10.0],btype='band',analog=0,output='ba')

The result was:

Warning: invalid value encountered in sqrt
(array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
    nan,  nan,  nan,  nan,  nan,  nan]), array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
    nan,  nan,  nan,  nan,  nan,  nan]))

I have no background in signal processing, so I actually don't know what I am doing. I don't know whether they used a IIR or FIR filter or whether I have to scale the cut-off frequencies or whether I'm using the wrong ripple. I hope you can help me.

The main issue with the example you gave is that the filter design function cheby1 is returning all NaNs, which isn't going to be a very good filter. The problem is how you're specifying the passband/stopband edge frequencies. This particular function is meant to emulate MATLAB's cheby1 function; the frequencies that you give it should be normalized, such that a value of 1 corresponds to half of the sample rate.

import scipy.signal as signal
fs = whatever_the_sample_rate_of_the_filter_input_is_going_to_be
signal.cheby1(8,0.05,[0.1/(fs/2),10.0/(fs/2)],btype='band',analog=0,output='ba')

I don't have SciPy handy, but that should at least correctly design the filter you want.


Deriving Time-Series Data from PSD when the Inverse Fourier Transform doesn't converge So I'm trying to get some sample time-series data from pink noise (PSD: $S(f) = 1/f$). However, when I try to evaluate the Inverse Fourier Transform integral

$\int \frac{1}{\omega} e^{j \omega t} dt$

from negative infinity to positive infinity, I get a non-converging integral. I tried a Taylor Expansion, but I'm running into the same issue. Is there some other way around this?

Even if you had succeeded in deriving the inverse Fourier transform of $1/|f|$ (note that you need the absolute value) then you still wouldn't have a time-series with a $1/f$ behavior. What you would get is the auto-correlation sequence.

A good way to generate pink ($1/f$) noise is to generate (approximately) white noise with a random number generator, and filter that noise with a filter whose frequency response approximates $1/\sqrt{f}$ over a large frequency range. Of course the frequency response needs to have a finite value at $f=0$, but it is possible to design a filter approximating $1/\sqrt{f}$ for $f>f_{min}$ with a relatively small value for $f_{min}$. If you have a question as to how to design such a filter you should probably formulate this as a new question. You can also have a look at this site.

And by the way, the inverse Fourier transform of $1/|f|$ does exist in terms of distributions: check this.


Do you really have to always satisfy Nyquist with all signals if your purpose is to reconstruct a signal? If not, when can you not comply with it?

No, because this is a sufficient condition (for regularly sampled signals), and not a necessary one. This condition restricts the space of all possible continuous signals to a subspace of discrete sequences that contain the same information.

Suppose that you can constrain the signal space, eg limited band-width, positivity, parametric models, sparsity, etc. Then many theoretical works ensure perfect reconstruction (in theory) under Nyquist-like conditions. Literature abounds, under the names of compressive or compressed sensing, sparse sampling, finite-rate of innovation, etc. Some are "exactly exact", some are exact within a high probability, under quantization bounds, etc.

The topic is very wide and active, I cannot provide more details unless the question is more precise. A couple of papers however:

A Systematic Review of Compressive Sensing: Concepts, Implementations and Applications, 2018
Sparse Sampling of Signal Innovations, 2008
A Review of Sparse Recovery Algorithms, 2018


Ideas for a college assignment I have to do a project for college. The project consists in writing a paper about some application of digital signal processing and then implement the solution in MATLAB. It doesn't have to be the original.

We're supposed to check papers to come with ideas but I thought that perhaps someone here could give some ideas.

Let me suggest a project for you that I think will fit with what you've learned and with the scope of the project, as far as I can judge it: Denoising of audio signals by spectral subtraction techniques.

The basic method is relatively simple. You need to cut your signal into overlappling blocks, you window each block, compute an FFT, and modify the magnitude of the FFT coefficients according to the estimated SNR in that frequency bin. High SNR means no modification, low SNR means attenuation by some gain factor $g$ satisfying $0<g<1$. For reconstructing the signal you combine the overlapping blocks (after IFFT) and you use the phase of the noisy signal, so only the magnitude is modified.

The nice thing about it is the simplicity of the basic method, but also that you can tweak a lot of things (gain computation, SNR estimation, etc.) to improve the perceptive quality of the denoised signal. If you do it in a very simplistic way you'll hear many artefacts, such as musical noise.

Also have a look at this and this answer for some more information on spectral subtraction.


Why can FFT only operate on images with specific properties? Can FFT only operate on Grayscale images? If Yes, why?
Can FFT only operate on images with dimensions of power of two? If Yes, why?

The FFT can only operate on grayscale images, yes. Each color images contains multiple channels and only if they are displayed together (depending on which color model you use or rather what channels you use), the image will become a full color image. A good example is an RGB image. It contains three grayscale channels with varying intensity. Each channel is colored and R+G+B added together. So if you want to use the FFT on a color image, you have to take it on each color separately or convert the image to grayscale.

The Radix-2 FFT can only be used with signals with $2^n$ samples, yes. The FFT takes advantages of the number of samples in the signal by reordering the signals samples and then taking the fourier transform: Even numbered samples and odd numbered samples. By cutting the size of the transform in half, you only need a quarter of multiplications as opposed to using the DFT on the full signal. But what stops us from reordering the samples again? - Nothing. So we reorder again and again and again until we end up with a "one bit DFT", do the transform and reorder all of the samples recursively. This allows us to keep the transform simple, but it only works with powers of 2.

Using a different Radix, for example 4, the same principle applies (Radix-4 further reduces the number of multiplications compared to Radix-2), but now the signal has to be $4^n$, a power of four: 4, 16, 64, 256, 1024, ...


Why do convolution kernels such as Gaussian, Laplacian, LoG almost always seem to be expressed in integers? I'm a total newb in search of some deeper understanding, but I'm not able to read the maths behind these on Wikipedia.

If I understand correctly, you get the new value for each pixel by multiplying each of its neighbours' values by the corresponding member of the convolution kernel/matrix. For Laplacian there seem to be two commonly used such kernels:

 0 -1  0
-1  4 -1
 0 -1  0

and:

-1 -1 -1
-1  8 -1
-1 -1 -1

But why do they use only integers? I believe modern devices have floating point as fast as faster integer math these days.

Is this just a leftover from the past when floating point was slow and now it's become a kind of tradition? Perhaps people just re-using the old info they find from those days without deeper understanding? Or is there an actual reason why integers are used?

In fact, won't the use of integers here result in artefacting? What am I missing?

The Laplacian kernel with the 4 in the middle results from summing second derivatives along the two axes ([1,-2,1]). Those are the right values to use, you can show this by writing out the math for the second derivative and set the distance h to 1 (or search for discrete approximation to derivative). This kernel hasn’t been rounded, the values just happen to be integer.

The other kernel, with 8 in the middle, does not provide correct values for the Laplacian, but uses more input values and hence is a bit more robust against noise (but not by a whole lot).


Transposed convolutional layer Can someone define the weights in a transposed conv2D kernel used to perform interpolation (NN or bilinear or whatever)?

The idea is to get "convinced" that one can perform upsampling (interpolation) by the process of transposed convolution (aka. "deconv layers" although this name is not entierly correct)

e.g. take a small image of size 2x2 e.g with grayscale values. (lets forget depth)

How would you defined your "transposed conv2D" kernel (shape, weights) in such as way as to perform this upsampling/interpolation?

1D nearest neighbor interpolation is $K = \begin{bmatrix}1 & 1 & 0\end{bmatrix}^T$ or $K = \begin{bmatrix}0 & 1 & 1\end{bmatrix}^T$. This depends on whether your library flips the kernel (i.e. cross correlation vs convolution). To increase the sampling frequency just increase the kernel size e.g. $\begin{bmatrix}0 & 1 & 1 & 1 & 0\end{bmatrix}$. See also my other answer.
For 2D NN interpolation compute the outer product $K^TK = \begin{bmatrix}0 & 0 & 0\\0 & 1 & 1\\0 & 1 & 1\end{bmatrix}$.


How to remove the boundary effects arising due to zero padding in scipy/numpy fft? I have made a python code to smoothen a given signal using the Weierstrass transform, which is basically the convolution of a normalised gaussian with a signal.

The code is as follows:

#Importing relevant libraries
from __future__ import division
from scipy.signal import fftconvolve
import numpy as np

def smooth_func(sig, x, t= 0.002):
    '''
    x is an 1-D array,
    sig is the input signal and a function of x.
    '''
    N = len(x)
    x1 = x[-1]
    x0 = x[0]

# defining a new array y which is symmetric around zero, to make the gaussian symmetric.
    y = np.linspace(-(x1-x0)/2, (x1-x0)/2, N)
    #gaussian centered around zero.
    gaus = np.exp(-y**(2)/t)

#using fftconvolve to speed up the convolution; gaus.sum() is the normalization constant.
    return fftconvolve(sig, gaus/gaus.sum(), mode='same')

If I run this code for say a step function, it smoothens the corner, but at the boundary it interprets another corner and smoothens that too, as a result giving unnecessary behaviour at the boundary. I explain this with a figure shown in the link below.
Boundary effects

This problem does not arise if we directly integrate to find convolution. Hence the problem is not in Weierstrass transform, and hence the problem is in the fftconvolve function of scipy.

To understand why this problem arises we first need to understand the working of fftconvolve in scipy.
The fftconvolve function basically uses the convolution theorem to speed up the computation.
In short it says:
convolution(int1,int2)=ifft(fft(int1)*fft(int2))
If we directly apply this theorem we dont get the desired result. To get the desired result we need to take the fft on a array double the size of max(int1,int2). But this leads to the undesired boundary effects. This is because in the fft code, if size(int) is greater than the size(over which to take fft) it zero pads the input and then takes the fft. This zero padding is exactly what is responsible for the undesired boundary effects.

Can you suggest a way to remove this boundary effects?

I have tried to remove it by a simple trick. After smoothening the function I am compairing the value of the smoothened signal with the original signal near the boundaries and if they dont match I replace the value of the smoothened func with the input signal at that point.
It is as follows:

i = 0
eps=1e-3
while abs(smooth[i]-sig[i])> eps: #compairing the signals on the left boundary
    smooth[i] = sig[i]
    i = i + 1
j = -1

while abs(smooth[j]-sig[j])> eps: # compairing on the right boundary.
    smooth[j] = sig[j]
    j = j - 1

There is a problem with this method, because of using an epsilon there are small jumps in the smoothened function, as shown below:
jumps in the smooth func

Can there be any changes made in the above method to solve this boundary problem?

Also I tried removing the zero padding in the fft source code and replaced it with a constant value, but it gave undesired results.
Can you suggest a way of removing this zero padding in the scipy fft source code?

Can there be any changes made in the above method to solve this
  boundary problem?

Yes there can be, but it all depends what you want to occur at the boundary. If you don't want it to be smoothed then simply replace the edge with the original signal for a sensible number of samples, depending on the width of your gaussian. You gaussian might be > 0 for say 10 samples, so replace the 5 samples at the start and end of the signal with the original samples.

  Can you suggest a way of removing this zero padding in the scipy fft
  source code?

Yes. Pad it yourself. If you want to consider the wrap around, as if your signal is truly periodic, then fabricate the wraparound effect for a power of 2 in size by carefully padding your own signal. Just as a simple idea, multiply your window by 2 and find the nearest power of 2, then fill the middle of the copy with zeros to pad it out

[1234567890] [12345] - zero pad here - [67890]

Now when you apply the convolution, you should get the same effect for the first window's worth of the result as if the convolution kernel had wrapped around the original signal. To illustrate, take your gaussian and zero pad it at the end, and slide through the convolution to see for yourself that it yields the same thing

 original   left    padding    right
[1234567890 12345 000000000000 67890]
[xgaussianx 00000 000000000000 00000]

so the value for the right most sample comes from this:

 original   left    padding    right
[1234567890 12345 000000000000 67890]
[sianx00000 00000 000000000000 xgaus]

Which is equivalent to what we actually want

[1234567890]
[xsianxgaus] (gaussian shifted by 5 samples and wrapped around)

and the value for the left most sample comes from

 original   left    padding    right
[1234567890 12345 000000000000 67890]
[00000xgaus sianx 000000000000 00000]

Which is also equivalent to what we actually want

[1234567890]
[sianxxgaus] (gaussian shifted by -5 samples and wrapped around)

All this relies on deciding what it is you actually want to happen. If you understand the mathematical process of convolution, and you understand exactly what it is you want to calculate, you can fabricate the conditions under power of 2 scenarios to ensure you get fast computation. Generally speaking most people just opt for the replace because their signal windows are not periodic and have sharp discontinuities at the edges.


DFT: Basis functions and Significance of dividing frequency by Sample length A time domain signal can be decomposed into sinusoids which are based on basis functions.

For a N sampled input, a cosine basis function is defined as:

$$C_k[i] = \cos\left[\dfrac{2\pi k i}{N}\right]$$

Range of Frequency $k = 0 \ldots \frac{N}{2}$, Range of $i = 0 \ldots N-1$

Why is the sinusoid angle divided by $N$?

This normalises the frequency such that the period of the first cosine is N samples. Then all subsequent cosines are N/2 samples, N/3 samples.... all the way to N/(N/2) = 2 samples, at the Nyqust limit where k = N/2. If you want the full DFT, including negative frequencies, k should range from 0 to N-1.


How to calibrate a signal with respect to frequency I have recordings from accelerometers (non-flat response with respect to frequency) and I'd like to convert my voltage time series into acceleration values.  Is there a standard method for doing this?

My proposal to achieve this is:

Break the signal into equally spaced segments overlapping with 50%.
Apply a fft (with appropriate window) to each segment and apply the calibration values to each frequency.
Convert the segment back into a time series by applying an ifft and reversing the window transform.  We not have calibrated segments of the time series, although they are likely disjoint between window segments.
Apply a triangle window to all overlapping segments and add them into a single time series.  The triangle window should result in smooth transitions between window segments, resulting in a continuous, calibrated time series.

Are there any problems with this approach?  The window choice is dependent on use case.  Are there other approaches to achieve the same goal?

Potentially a lot of problems.

You really want to apply a filter to your data that the calibration curve represents, which seems to be your conceptual approach, only frequency domain filtering doesn't involve overlapping and windows, which one does in spectral estimation.

The hard part is to convert your calibration curve into a filter.


Negative exponential signal's energy and power Earlier I have dealt with exponential functions multiplied with unit step function. But, energy and power of exponential function alone comes out to be infinite when I put limits of the integrals to from $-\infty$ to $+\infty$. How can I find its energy or power if the signal is not multiplied with unit step function?

An ideal exponential signal $x(t)=e^{at}$ , which extends from $-\infty$ to $\infty$ has infinite energy and infinite power, as for real $a >0$ ( and similarly for real $a < 0$) you have

$$
\mathcal{E}_x = \int_{-\infty}^{\infty} e^{2at} dt = \lim_{t \to \infty} \frac{1}{2a} e^{2a t} \to \infty
$$

and

$$
\mathcal{P}_x = \lim_{T \to \infty} \frac{1}{T} \int_{-T/2}^{T/2} e^{2at} dt = \lim_{T \to \infty} \frac{e^{a T}}{2aT}  \to \infty
$$

Therefore It's neither energy nor power signal. Note that it's not a practical signal and exists only for mathematical purposes.


In MATLAB ,do the input signals need to be of the same sampling rate for the conv function to work correctly? eg. z=conv(u,v) where u is the input signal encoded in 24KHz sampling rate and v ,the impulse response encoded in 44kHz .
z is played at 24KHz rate .
In the convolution function ,we basically perform a convolution of input and impulse response samples.But reading across the net, I found that for convolution the input and impulse signals must be signals encoded at the same rate.
How true is this?

What an interesting question! An impulse response sequence, taken on its own, does not have a specific numerical sample rate (measured in Hz) associated with it. If I applied a [1, 0, 0, 0, ...] unit impulse sequence to a 3-point moving average lowpass filter then the filter's impulse response would be the sequence [1/3, 1/3, 1/3, 0, 0, 0, ...].  There's no specific numerical sample rate (measured in Hz) associated with either the input unit sample sequence or the filter's impulse response sequence. If you insist on talking about the somewhat unusual notion of the sample rate of a filter's impulse response I suppose you can say, "The sample rate of a filter's impulse response is always equal to the filter's input data sample rate."

For example, if to my 3-point lowpass moving average filter I apply an input sequence whose sample rate is 1000 samples/second (1 kHz), then the filter's output sequence has a sample rate of 1000 samples/second. On the other hand, if I apply an input sequence whose sample rate is 24000 samples/second (24 kHz) to the same lowpass filter, then the filter's output sequence will have a sample rate of 24000 samples/second. (WARNING: I suggest you never discuss the sample rate of a filter's impulse response during a job interview.)


Quadrature modulation I have a question:
What is the difference between
$ I \cos(t) - Q \sin(t) $    and   $ I \cos(t) + Q \sin(t) $

except that they are phase shifted 90 degrees ?
On the internet I can find both modulators with adding quadrature data and substracting quadrature data.

http://www.ni.com/tutorial/4805/en/  at the end
http://zone.ni.com/reference/en-XX/help/372058G-01/nirfsa/iq_modulation/

thank you

The idea is that you can transmit two signals, $I(t)$ and $Q(t)$, over the same bandwidth and at the same time, and still recover each independently of the other. The math is pretty simple. If the transmitted signal is $s(t)=I(t)\cos(2\pi f_ct)-Q(t)\sin(2\pi f_ct)$, then (ignorning noise, ignoring a factor of 1/2 and assuming coherent reception) the receiver can calculate
\begin{align*}
r(t)\cos(2\pi f_ct) &= I(t)\cos^2(2\pi f_ct)-Q(t)\sin(2\pi f_ct)\cos(2\pi f_ct)\\
&=I(t)+I(t)\cos(2\pi 2f_ct)-Q(t)\sin(2\pi 2f_ct).
\end{align*}
After a low-pass filter, the receiver has recovered $I(t)$. Likewise,
\begin{align*}
r(t)\sin(2\pi f_ct) &= I(t)\cos(2\pi f_ct)\sin(2\pi f_ct)-Q(t)\sin^2(2\pi f_ct)\\
&=I(t)\sin(2\pi 2f_ct)-Q(t)+Q(t)\sin(2\pi 2f_ct).
\end{align*}
from which $Q(t)$ can be recovered by an low-pass filtering and multiplication by -1.

The identities used here are: $\cos(A)\cos(B)=0.5(\cos(A+B)+\cos(A-B))$, $\cos(A)\sin(B)=0.5(\sin(A+B)-\sin(A-B))$ and $\sin(A)\sin(B)=0.5(\cos(A-B)-\cos(A+B))$.

Of course, you can also define the transmitted signal as $v(t)=I(t)\cos(2\pi f_ct)+Q(t)\sin(2\pi f_ct)$ and avoid the sign inversion when estimating $Q(t)$. The main reason not to, though, is that the complex envelope of $s(t)$ is $I(t)+jQ(t)$. Since most of the time it is preferrable to work with the complex envelope, $s(t)$ is chosen to produce the most convenient expression for it.


Consider an ideal low pass filter $H(\omega)$, and the input to this filter is the periodic square wave $x(t)$. Find the output $y(t)$ The solution to the problem is $$y(t) = 5 + \frac{20}{\pi} \sin(\pi t) + \frac{20}{3\pi} \sin(3 \pi t) $$ and to get that the solution says to find the Fourier series expansion of $x(t)$ and I am having trouble doing that. I find that $C_0 = 5$ and that $$C_k = \frac{5}{j k \pi} \left( 1 - e^{-j k \pi} \right)$$. I don't know how to compute $$C_0 + \sum_{k=1}^{\infty} C_k e^{j     k w t} $$. Also since the cut off is $4 \pi$ so when the Fourier series expansion greater $4 \pi$ those values will not be included as seen in the solution...

So you actually solve(d) this problem by computing the CTFS of the input $x(t)$:
$$ C_k = 5 \frac{ 1 - e^{ - j k \pi } }{ j k \pi } = \begin{cases} 0 &, \text{ for k even , } k \neq 0 \\ 10/jk\pi &, \text{ for k odd} \end{cases} $$ which gives the coefficient associated with the complex exponential of frequency $\omega_k = k \omega_0$ for the fundamental frequency $\omega_0 = \frac{ 2 \pi }{T_0} = \frac{ 2 \pi }{2}  = \pi$ , as the period of the input $x(t)$ is $2$ seconds.

Since the ideal lowpass filter passes those frequencies below its cutoff frequency $\omega_c = 4 \pi$, then we would only take those coefficiencs whose frequencies fall below $\omega_c$ yielding those values of $k = 0,\pm 1, \pm 3$ . Those indices higher than $\pm 5$ are filtered out by the ideal lowpass filter.

Therefore the output $y(t)$ is:
$$ y(t) = 5 + \frac{10}{j\pi} \left( e^{j\pi t} - e^{-j\pi t} \right) + \frac{10}{j3\pi} \left( e^{j3\pi t} - e^{-j3\pi t} \right) $$ which simplifies to:
$$ y(t) = 5 + \frac{20}{\pi} \sin(\pi t) + \frac{20}{3\pi} \sin(3\pi t) $$


Reducing a signal to Acos(omega * t+ phi) with negative omega*t I am trying to solve a phasor addition problem, reducing from its original form to
$$X(t) = A \cos(\omega_0 * t + \phi)$$

The original equation is :
 $$X(t) = 2 \sin(\omega_0 * t + 45) + \cos(\omega_0*t)$$

I converted the sine term into cosine. The result that i got was
$$X(t) = 2 \cos(\pi/2 - \omega_0*t - \pi/4) + \cos(\omega_0*t) $$
Which reduces to :
$$X(t) = 2 \cos(\pi/4 - \omega_0*t) + \cos(\omega_0*t) $$
Up to this point I verify in MATLAB that the original function and the function in cosines are equal by graphing them and then calling the $norm()$ function and verify that it is something very close to 0.
http://www.mathworks.com/help/matlab/ref/norm.html

I then perform the phasor addition process.
$X(t) = Re(2e^{j(\pi/4)}e^{j\omega_0 t})+Re(e^{j0}e^{j\omega_0 t})$

$$X_1= 2e^{j(\pi/4)} = 2(cos(\pi/4) + jsin(\pi/4)) = 2(\frac{\sqrt2}2+j\frac{\sqrt2}2) = \sqrt2 + j\sqrt2$$
$$X_2 = e^{j0} = cos0 + jsin0 = 1$$
$$X_3 = X_1 + X_2 = \sqrt2 + j\sqrt2 + 1 = 2.41421 + j1.41421$$
 To turn it into back into the stated form of $X(t) = A \cos(\omega_0 * t + \phi)$ i find the $R$
$$R = \sqrt{2.41421^2 + 1.41421^2} = 2.79793$$
$$\theta = arctan(1.41421/2.41421) = .529902$$

The final result being $X(t) = 2.79793\cos(\omega_0 t + .529902) $

I then go to MATLAB to plot the original signal and the reduced signal that I found but when I call the $norm(originalSignal - reducedSignal)$ I get a number much greater than 0.

I am trying to figure out why this is returning such a high number. Some of my guesses are the negative $-\omega_0 t$ in the first term of the function when I convert sine to cosine $X(t) = 2 \cos(\pi/2 - \omega_0 t - \pi/4) + \cos(\omega_0 t) $ I don't know if the process changes for a situation like that. Problems I have worked through and have gotten the correct answer always seem to have the $\omega_0 t$ as positive part inside the cosine term.

Below is the Matlab code.

t = [0 : .01 : 1];
originalWave = 2*sin(50 *t + pi/4) + cos(50*t);
reducedWave = 2.79793*cos(50*t+.529902);

subplot(2,1,1)
plot(t, originalWave)
subplot(2,1,2)
plot(t, reducedWave)

norm(originalWave - reducedWave)

Any ideas on what could be wrong with my process or answer?

I appreciate any advice or help. Thanks.

You flipped a sign on the phase, the answer is -.529....
It's not
$$X(t) = Re(2e^{j(\pi/4)}e^{j\omega_0 t})+Re(e^{j0}e^{j\omega_0 t})$$

It's
$$X(t) = Re(2e^{-j(\pi/4)}e^{j\omega_0 t})+Re(e^{j0}e^{j\omega_0 t})$$


Pattern recognition in time series 4x3000 vector I have a vector, here is a sample of some data from some heat flow data:

I would like to identify features in this image. In the example above I have identified one feature I would like to automatically recognize. There are others I would like to identify. In this example, the data in all four vectors moves up at the same time, however this can happen at different rates and sometimes it can be delayed.

Like this one:

My regular set of tools is inadequate for this type of problem. I have no problem detecting peaks or using thesholding. This is the first time I've had to do this with multi dimensional data. The other problem is the data is very noisy so peak detection and thresholding yeilds false positives.

I would like to plug this into a neural network and classify it, but the problem is in this 4x3000 vector I might only have 4 or 5 examples of features that I could train with. So it might be hard to train a CNN if I don't have many examples.

There are many tools that work with 1x(n col) vectors, but I need to detect patterns that coincide in all four vectors. Another problem I have is

How can I compare a 4x5 section of the 4x3000 matrix and find similarities?

The think I would really like to do is take a 4x5 matrix and compare it with the rest of the vector and see if and where it matches (so really there would be 2995 matrices of size 4x5). Are there neural networks that can make multidimensional comparisons like this?

Assuming your time series are the same length, take your data and produce spectrograms, $S_k$, for each row of the $4 \times 3000$ data matrix $D$ that you have. Since these individual time series should be similar, this means you can try to extract features by stacking these spectrograms next to each other horizontally into a large matrix. So if you have $S_k$ as the spectrogram for row $k$, you can create the overall spectrogram as $S = \left[ S_1, \cdots, S_4\right]$.

From here, find the modulus of the elements (since they will be complex), and define the result as $\bar{S} = |S|$. Then I recommend performing a Non-negative Matrix factorization on $\bar{S}$ such that you get the low-rank approximation $\bar{S} \approx FW$ with some small specified number of features. $F$ in this case will be your features (they will be column vectors in $F$) and $W$ will be your new data in feature-space. You can look at $W$ to see where a given feature from $F$ appears. Given you've generated $F$, all you need to do to do to find $W$ for individual/new time series results is compute the spectrogram for that time series, call it $\hat{S}$, and then do $\hat{W} = F^{+} |\hat{S}|$ where $F^{+}$ is the pseudo inverse of $F$. Now you can look at $\hat{W}$ to get an idea of where certain features are being used and also build models using this lower dimensional representation.

Note that the number of features you use will have to be played with. I would look at the modulus of a sample spectrogram to get a feel for how much the frequencies are changing over time and consider roughly how many dominant features you think could be used to approximate it and then set the number of features as that to start. As an example, I have done this stuff above to identify different sounds in an audio clip and I only had to reduce it to the number of dominant sounds (there were 3 of them) so that hugely minimized my dataset size. Images from this example are shown below:

One thing you might consider as well is taking the separate data streams and trying to fuse them using some sort of filtering algorithm or using Hidden Markov Models. Then just perform the above spectrogram computations on the filtered result.


When should we using moving average in algorithm design? I'm new to signal processing. And I just read G.720.1. I find there is a lot of moving average used.
And I can also recall there is a lot of moving average been used in other audio processing algorithm.

For example in G.720.1, the ZCR (The zero crossing rate of the input 10ms frame) are input of moving average algorithm and the result are used for sound detect. But I assume ZCR should have big change in continues frames.

So when will you think you need moving average in an algorithm design?
And if we need one, how we should select the weight?

So when will you think you need moving average in an algorithm design?

If you mean moving average filters; moving averages as the name suggest are computed as averages of samples say,  $M-1$ previous samples (+ current sample) from input $x[n]$ to get an average output $y[n]$, repeating the process to get all $y$ samples. Computing their mean to get an averaged $y[n]$ signal. What you get is a smoothing effect.

As explained by Sanjit Mitra in the chapter on Discrete-Time Systems in his book. When your data is corrupted with noise and you happen to have multiple measurements of the same data samples, a good estimation of your signal is doing the ensemble average. But if you only have one set of measurement, then moving averages is a way of doing it. The result is a smoothing of your noisy signal. An M-point moving average is defined as a FIR system as:

$$
y[n] = \frac{1}{M}\sum_{k=0}^{M-1} x[n-k],\quad 0\leq n-k \leq M-1
$$

  And if we need one, how we should select the weight?

As seen in the equation, all the weights are equal and dependent on the length of the chosen filter:

$$
h[n] = \frac{1}{M}
$$

The choice of the length is critical. An efficient implementation however is done by introducing a recursive form and requires only 2 additions and 1 division.

Again, depending on your application, there are other variants  of the moving average that may be more suited to your signals.


filterbank: understand the different responses at the center frequency of each filter I'm computing filterbank by applying 26 triangular filters on a Mel-scale to the power spectrum of an audio frame to extract frequency bands, and I found that some references use filterbank with equal responses at the center frequency of each filter, see the following figure:

A formula for calculating these is as follows:

Then, they take the log of the resulting energies to compute the log filter bank energies.

But, in some other references, they use filters with decrease responses at the center frequency of each filter, i.e. something like this:

I want to understand the difference between these two methods, and which one is more accurate for usage in a speech recognition application based on deep learning?

And what are the formula and different steps to calculate the second one?

References:

http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/

http://www.cs.cmu.edu/afs/cs/user/bhiksha/WWW/courses/yahoo2009/01-02.featurecomputation.ppt

Normalizing the filterbanks by their widths is optional and totally up to you (similarly to the warping scale Mel/Bark). Depending on your application, you can start without normalization and see what results you are getting. Personally I prefer to keep it fixed and have one knob less for turning. There are more important parameters to tune, such as warping scale, number of banks and coefficients.

For example in the case of DNN's normalization is redundant, since you will be normalizing the input features anyway (at least that's what you should do), so any way you will end up with the same distribution.


Problem reading and rendering Spectrum I'm now processing PCM codes extracted from .wav file and want render the spectrum myself. First I try to read a spectrum on my foobar2000 player as below.

My problem with this spectrum is: what do x-axis and y-axis represent? IIRC, the x-axis represents different frequency(or angular speed) and the scale should be whole-number times of fundamental, i.e. 50. However, when I divide next by its preceding, e.g. $94\div 69$, the outcome is roughly always 1.37. As to the y-axis, why the scale is negative? What is the calibrated value, i.e. $P_2$ in $10\log_{10}{\frac{P_1}{P_2}}$.

If I want to render the amplitude spectrum with PCM code, should some normalization be first taken(I know the original 8-bit PCM code ranges from -128 to 127, even normalized range is -1.0 to 1.0. It seems I can't use these values directly to represent amplitude of the signals)?

Thanks and Best Regards!

The x-axis is indeed frequency, but the numbers are quite strange to me. One would expect a scale in powers of 2 (octaves) or of powers of 10 (decades).

The y-axis represents the relative magnitude of the signal at that particular frequency. The 0dB level is what you would get from a sine wave with the maximum amplitude (-1..+1). So your P2 is 1.0.


What's the relation between the sampling rate and FFT size in OFDM In OFDM system,

Is there relationship between Sampling rate and FFT size? For example, if we downsample or oversample the received signal, Should we decrease or increase the FFT size accordingly?

thank you

Yes, you can use different sample rates and corresponding FFT sizes to accomplish the same thing for OFDM (mod and demod). For example, 10 MHz LTE can be sampled with 15.36 Msps and modulated/demodulated with a 1024-point FFT for 15 kHz subcarrier spacing and  the 600 subcarriers around DC give the 9 MHz of utilized BW. Alternatively, you can use 30.72 Msps and 2048-point FFT to mod/demod and get the same 600 subcarriers.

What is effectively happening is that the basis vectors are resampled (by 2x in this case), but they remain orthogonal and at the same subcarrier spacing.


RF demodulation for dummies Given a signal that's been FM encoded on a carrier wave, what are the individual steps I need to take to demodulate it, and recover the original signal?

I don't have a background in signal processing, I have a background in programming. I'm trying to solve a real-world problem: I'm trying to demodulate a raw analog RF signal captured from a Laserdisc at 40 MHz. This RF signal contains a composite video signal centered at 8.1 MHz, two analog audio channels at 2.81 MHz and 2.3 MHz, and one digital data channel encoded at 1.75 MHz and below.For a handy reference image, here's a visual representation:

More specifics on the frequency ranges and other related info can be found in a concise summary here, and in the official documentation describing the final encoding format here.

Now, I've captured a signal like this straight from a Laserdisc. Under a spectral view, it looks like the following:

Here I can clearly see the composite video signal encoded, in this case showing the end of a field, the vblank period, and the start of a new field. I can see the two analog audio channels encoded below, and at the bottom, a band of low-frequency noise which contains the digital channel data. I can also see the harmonic sidebands coming through, giving lower intensity "echoes" of the data at lower and higher frequencies than the encoded frequency.

When the composite video is decoded to a waveform, I expect to see a signal like this:

The problem I'm having right now is, I don't really understand the math involved in demodulating the FM signals from the combined RF signal here.

I understand I can use some basic band-pass filtering to eliminate frequencies outside the bands I'm interested in, which would allow me to separate the audio, video, and digital signals from each other, but after that it's the actual process of de-modulating the FM-encoded data back into their data streams (i.e., turning FM encoded data back to a basic waveform), that I'm not getting right now.

For an answer here, I'm after something that explains this in layman's terms, so that if I was a green student walking into a classroom on signals processing for the first time (which I pretty much am here) I can understand what you're talking about. I still need it to be specific though, because ultimately, I need to understand this well enough so that I can write real code that performs these operations. I'm after a pure digital solution to this problem.

A good answer might constitute a few chapters in a good DSP text, making this question overly broad.  So I'll sketch an over-simplified skeleton, where you will need to seek a decent rationale elsewhere.

Pretend your captured signal is just a vector (array of samples).

Move your signal of interest "down" to "baseband" (resulting in either 2 vectors (X,Y) or 1 complex vector (if you know complex arithmetic).
Do this by multiplying by a vector of cosines, oscillating at -8.1 MHz, to get vector X, and by a vector of sines, oscillating at the same (negative!) frequency, to get vector Y.  Low pass filter (to remove audio etc.) these 2 vectors with a cut-off somewhere around 2 or 3 MHz (see specs for your video).

The angle of a 2D vector of X,Y pairs (use atan2() to compute) will rotate with a speed proportional to the modulating signal (e.g. the video signal you want).  The speed of rotation is the first derivative of the angle with respect to time.  Or roughly, subtract successive angles, to get a time vector proportional to dA/dt.  Unwrap/rewrap phase discontinuities as needed.  Low pass filter (maybe somewhere around 2 or 3 MHz) to remove noise.  Scale as needed.  Result should be a vector containing your raw video signal (perhaps at the wrong scale, perhaps noisy, perhaps inverted, perhaps with/without gamma: fix if possible).


MATLAB: filter nonperiodic noise in the signal I have gone through the below related posts of filtering noise from the signals

MATLAB, How to filter a discrete signal?

MATLAB: Filter a noisy EKG signal

But the signal I want to filter contains noise which is not periodic as shown in the images below:

Can anyone suggest a filtering technique for the noise labeled with red colour in the image?

You can do a non-decimated one dimensional wavelet transform on the image. Since the data is more or less homogeneous within a certain range, the outliers will be denoted by spikes in the wavelet coefficients. You can then use a threshold and remove those coefficients and perform inverse non-decimated wavelet transform.

You can check this paper also.

Do you mind sharing a copy of the data in .mat ? I would like to test it out in order to be sure.