How to effectively teach topic of frequency domain? This is a soft question. The frequency domain has made revolution in the field of Mathematics,  Digital Signal and Image Processing etc. Some of concepts which are very difficult to analyse in spatial or time domain can be very easily understood in the frequency domain using theorems like Fourier transform & series, wavelets, Laplace transform, etc.

Well, I am doing under graduation in 4 years Electronics and Telecommunication Engineering. We have semester pattern system (2 semesters in 1 year). In the 3rd semester, we have subject named "Signals and Systems" under which the above topic is covered.

But my experience is that when it is being taught in the academic curriculum or course, it is not given as much importance as it has to be. Most of the teachers or professors don't know how teach this topic. They teach with the help of ready-made PowerPoint presentations (e.g. for Laplace transform, single PPT which is mostly downloaded from internet) and sometimes they take chalk in their hand. So the students like us get hardly any knowledge about the frequency domain in academic course compared to discussions in the forum like the very useful http://www.DSP.stackexchange.com and http://www.math.stackexchange.com.

So can anybody give explanation about how to effectively teach topic of frequency domain? You can also explain with the help of any subtopic like Fourier domain , Laplace or wavelets etc.
So that all of us would study this difficult topic with interest.

In my personal opinion, it's best to study how the master teachers of the field have presented this information, and to the extent possible, copy their approaches and their style.

There are two teachers in the field who have more or less taught a generation of DSP engineers their craft, either by direct teaching, or by teaching teachers.  They are Alan V. Oppenheim and Richard G. Lyons.  In Oppenheim's case, the Signals and Systems courses that he taught are still available on MIT's web site for free viewing, and his Discrete-Time Signal Processing is still the bible for many engineers on the topic.  Lyons still teaches and writes from Besser Associates, and his "Understanding Digital Signal Processing" book is simple and lucid.

You can't go too far wrong by copying the masters.


What is meant by LFM (Linear Frequency Modulation) waveform? I m a newbie in radar experiment on which currently i m working, and I m trying to understand the transmitted waveform which is LFM (Linear Frequency Modulation), but unable to interlink between the following parameters of transmitted signal waveform ,can someone explain to me how do these are connected to LFM?

Center frequency  : 150kHz
Sampling Frequency: 4000kHz
Bandwidth         : 50kHz
Pulse length      : 500microseonds

An LFM pulse is one in which the "instantaneous frequency" changes linearly over the duration of the pulse.  By "instantaneous frequency" I mean the rate of change of phase.  Over its $\tau=500\mu s$ duration, a pulse with "bandwidth" $B=50 kHz$ and center frequency $f_c=150 kHz$ would be defined by

$f(t)=A\cos(\theta(t))=A\cos(2\pi(f_c-B/2)t+\pi(B/\tau)t^2+\theta_0)$

Taking the derivative of phase with respect to time, we have

$\partial\theta(t)/\partial{t}=2\pi(f_c-B/2)+2\pi(B/\tau)t $

We can see that this sweeps from $2\pi(f_c-B/2)$ to $2\pi(f_c+B/2)$ over the duration of the pulse.

An LFM pulse provides better range resolution than an unmodulated pulse of the same duration.


Complex conjugate and IFFT I asked a question over on stack overflow.

I'm having a slight problem however.  As suggested by Paul R I am mirroring my lower $n/2$ bins into the upper $n/2$ bins.

I have a few questions however.

I'm assuming $n$ is $2^{fftOrder}$?  Thus for an order 10 IFFT I would be mirroring the bottom 512 frequency bins in to the top 512?  This would also produce me 1024 real samples in the end?
When mirroring does this mean that bin 0 would go to bin $n - 1$, bin 1 to bin $n - 2$ and so on?
I appear to be getting some form of reflection at a quarter of my sampling frequency?  This means that the top half of my frequency range is actually appearing over the bottom half.  Has anyone any ideas why?
Does anyone have any idea why I am getting a striping effect where I get 512 samples of one magnitude and the next 512 of a different? (This "could" be where my mirroring is coming from)

I think thats it for now.  Any help would be hugely helpful!!

Edit: Giving you an image of the reflection will be near enough impossible as its very hard to see!  Basically I'm sampling at 22 kHz and getting nothing but "noise" above 11 kHz.  What "ought" (though obviously I'm doing something wrong) to be above 11 kHz seems to get reflected down from 11 kHz to 0.  Its very odd.

As for the image of the stripes they can be seen here

Each stripe, I believe, is 512 samples wide.

I'm not sure what you're getting at with "reflection" or "striping," but you want to duplicate your data points such that the frequency-domain signal (the input to the IFFT) is conjugate symmetric. That means, for $N$ even:

$$
X[k] = X^*[N-k], k = 1, 2, \ldots, \frac{N}{2}-1
$$

Note that $k$ is zero-based; the first element in the vector is $X[0]$. The above indicates that you do not duplicate the zeroth and $\frac{N}{2}$-th elements in the input vector. This relationship is based on the periodicity of the discrete Fourier transform; the first half of the frequency-domain vector corresponds to "positive" frequencies covering the angular frequency range from $0$ to $\frac{(N-1)\pi}{N}$. The second half covers from $\pi$ to $\frac{(N-1)2\pi}{N}$, which, taking advantage of the $2\pi$ periodicity in the frequency domain, is equivalent to $-\pi$ to $\frac{-2\pi}{N}$. In order for the frequency-domain vector to have conjugate symmetry, then it can be easily shown that the above relationship must be true.

Speaking to your original goal of converting a spectrogram into an audio signal that sounds meaningful, there's no guarantee that what you get out will be in any way pleasing. As was pointed out on Stack Overflow, the spectrogram has no phase information; the method suggested there assumes the phase of each frequency-domain bin is zero (so the conjugate operation shown above is superfluous). While the human ear isn't terribly sensitive to phase distortion in audio, you might have to do something else depending upon exactly what you're hoping to achieve.


How Is the Formula for the Wiener Deconvolution Derived? Wikipedia - Wiener Deconvolution shows this formula:
$$
\ G(f) = \frac{H^*(f)S(f)}{ |H(f)|^2 S(f) + N(f) }
$$
Could someone explain the derivation?
Specifically, where does the squaring ($|H(f)|^2$) come from?

The Wiener Filter can also be derived by another (Easier) way.

Let's assume the following model:

$$ y = h \ast x + n $$

Namely the data is a result of a linear combination (Convolution) of $ x $ with Additive Noise.

If we assume the noise model is Gaussian and our data is also formed by a Gaussian distribution then we should try to minimize (MAP Estimator):

$$ \hat{x} = \frac{1}{2 {\sigma}_{n}^{2}} \left\| h \ast x - y \right\|_{2}^{2} + \frac{1}{ 2 {\sigma}_{x}^{2} } \left\| x \right\|_{2}^{2} $$

Since the convolution is Linear Operation and assuming we in finite dimension then it can be rewritten as:

$$ \hat{x} = \frac{1}{2} \left\| H x - y \right\|_{2}^{2} + \frac{{\sigma}_{n}^{2}}{ 2 {\sigma}_{x}^{2} } \left\| x \right\|_{2}^{2} $$

This is just a Tikhonov Regulaized Least Squares with the solution given by:

$$ \hat{x} = \left( {H}^{T} H + \frac{{\sigma}_{n}^{2}}{{\sigma}_{x}^{2}} I \right)^{-1} {H}^{T} y $$

Now, if we use the discrete form of the circular convolution then $ H $ is a circulant matrix which means all can be solved in the Frequency Domain which results in the same equation as above.


Why the requirement of the GCD of the lengths of all circuits in the graph being one? I am reading A Mathematical Theory of Communication. The second requirement of an ergodic process confuses me (emphasis mine):

  All the examples of artificial languages given above are ergodic. This
  property is related to the structure of the corresponding graph. If
  the graph has the following two properties the corresponding process
  will be ergodic:

  The graph does not consist of two isolated parts $A$ and $B$ such that it is impossible to go from junction points in part $A$ to
  junction points in part $B$ along lines of the graph in the direction
  of arrows and also impossible to go from junctions in part $B$ to
  junctions in part $A$.
  A closed series of lines in the graph with all arrows on the lines pointing in the same orientation will be called a “circuit.” The
  “length” of a circuit is the number of lines in it. Thus in Fig. 5
  series BEBES is a circuit of length 5. The second property required
  is that the greatest common divisor of the lengths of all circuits in
  the graph be one.

The author later elaborates

  If the first condition is satisfied but the second one violated by
  having the greatest common divisor equal to $d > 1$, the sequences
  have a certain type of periodic structure. The various sequences fall
  into $d$ different classes which are statistically the same apart from a
  shift of the origin (i.e., which letter in the sequence is called
  letter 1). By a shift of from $0$ up to $d - 1$ any sequence can be
  made statistically equivalent to any other. A simple example with $d =2$
  is the following: There are three possible letters $a, b, c$. Letter
  $a$ is followed with either $b$ or $c$ with probabilities
  $\frac{1}{3}$ and $\frac{2}{3}$ respectively. Either $b$ or $c$ is
  always followed by letter $a$. Thus a typical sequence is

  $$a\,b\,a\,c\,a\,c\,a\,c\,a\,b\,a\,c\,a\,b\,a\,b\,a\,c\,a\,c.$$

  This type of situation is not of much importance for our work.

But unfortunately, I still didn't get his point. Can you show me how to perform "a shift of from $0$ up to $d - 1$" to make two sequences statistically equivalent? In fact, the process already seems ergodic to me, but that's just my intuition, and I cannot justify my argument with rigorous reasoning.

In addition, what makes the above-mentioned process not of much importance in communication?

This is my understanding: the statistics of the source described in the paper depend on which character is produced first. If the first character is $a$, then one of the source properties is that letters in odd positions are always $a$. However, if the first letter is $b$ (in other words, a shift of $1$ in the circuit), then letters in odd positions can be either $b$ or $c$. This is what makes the source non-ergodic.

Note that the statistics when choosing $b$ as first letter and $c$ as first letter are exactly the same. In other words, there are $d=2$ different classes of sequences, and those in each class have the same statistics.

Regarding the importance of these sources: for one, Shannon's theorems depend on the source being ergodic. From a more practical standpoint (and this is purely my interpretation), non-ergodic sources are not good models for real-world processes; referring back to the example, no practical language produces repeated letters in a fixed periodic pattern.

I recommend reading "An Introduction to Information Theory: Symbols, Signals and Noise" by John R. Pierce, along with Shannon's paper. It explains many of the concepts in the paper in a very didactic manner. For example, it takes four pages to very clearly explain ergoidc sources, with nice examples.


Confused about DFT symmetric behavior for real-valued signals? I have a simple and might be a preliminary question about symmetric property of DFT of real valued signals. From definition we know, DFT coefficients of a real-valued signal is conjugate symmetric, so $X_{N-k}=X_{k}^*$ . However, I cant see this in calculated DFT, for example consider the below:

N=1000;
fs=10^6;
t=(0:N-1)./fs;
x=2*sin(2*pi*t*10000);
plot(x);
figure;
stem((abs((fft(x)))));

The non-zero bins are 11 and 991, where I expected them to be 11 and 989, even considering that indexes starts at 1 in MATLAB,and subtracting bin numbers by 1,  the non-zero bins must be 10th and 990th bins.

Would any one please point out, where is my misunderstanding?

$N=1000$ and your signal frequency is $10000$. So there should be a peak at bin $k=10$. And another one at $N-k=1000-10=990$. However since MATLAB starts indexing at $1$ instead of zero, you see peaks at $11$ and $991$.


Frequency warping when integrators are replaced with backward-euler and forward-euler integration Resonant controllers are used in the power industry. The transfer function
is
$G_{res}(s) = \frac{K_i s}{s^2 +2\omega_c s+ \omega_o^2}$
The "textbook" discrete implementation is depicted in the image The first integrator is a forward Euler and the second one a backward integrator. It is not possible to use Tustin integrators because this would create algebraic loops. Note that $\omega_o$ is not a constant, it will vary slowly and will be fed by a PLL or FLL.
Here's my question. How can one characterize the frequency warping?
In one paper, they evaluate frequency warping by applying the bilinear transform to $G_{res}(s)$, I feel like this is wrong.

To the extent this helps, here are some interesting magnitude and phase plots of the backward Euler (here using the method of impulse invariance mapping which for 1/s has the same result as backward Euler but as Matt L points out in the comments is not generally the case) vs Tustin (method of Bilinear Transform mapping) integrators. Notice the Euler on the left has dominant phase distortion while the Tustin on the right has more dominant amplitude distortion (no phase distortion from an analog integrator at all!).

The analog frequency response for a pure integrator as $1/s$ is -20 dB/decade in magnitude ($1/f$) and -90° for phase.


Odd problem with really long sinc filters I'm trying to build a "perfect" low-pass filter as an exercise to test that I've understood some basics in signal processing.

I've computed the required length of a sinc function such that at either end the value is 0 for a given integer bit depth of output. For example, 16-bits would require a filter that is 41722 taps.

I'm then applying this filter to a frequency sweep file sampled at 48kHz, with Fc set at 0.5, thus it should act as a low-pass filter with the cutoff frequency at 12kHz.

When I feed this filter some data I get unexpected results. For the 16-bit filter, the low-pass effect is clear, but there is data beyond 12kHz as well as other corruptions. For a 20-bit filter, 667543 taps, the output is totally wrong.

In the 16-bit case, applying a blackman window to the sinc function improves the output and makes it look correct, but this doesn't help in the 20-bit case.

I've attached a couple of images to illustrate the problem. I'm fairly certain that my sinc function is generating correct data. I'm using vDSP_convD to apply the filter (I've also tested with a handwritten function, with no change.)

Am I misunderstanding how a sinc function can be used as a low-pass filter, or am I running into other problems?

This is the output of the 16-bit, 41722 tap filter:

This is the output of the 20-bit, 667543 tap filter:

This a graph of -200..200 on the 41722 tap filter:

This is the frequency part of the FFT of the centre of the 41722 tap filter.

A frequency sweep has non-zero bandwidth over the length of any FIR filter.  Use a long enough filter (compared to the sweep rate) and nearly half the spectrum of the sweep  could end up on each side of the cutoff.

Also, use a narrow enough transition band, and you might end up with Gibbs horn ringing right near the edge of the transition frequency.


Discrepancy between Gaussian FFT and its Fourier Transform I am trying to do the FFT of a Gaussian signal and comparing it to the theoretical Fourier transform. For infinitely small time step $dt$ and infinitely long signal length $T$, the 2 should become identical right?

First consider the theoretical DTFT of $x[n]$
$$ X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j \omega n} = \sum_{n=0}^{\infty} a^n e^{-j \omega n} = \frac{1}{1 - a e^{-j \omega} }$$

Then consider the N-point DFT of the truncated signal:
$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} n k} = \sum_{n=0}^{N-1} a^n e^{-j \frac{2\pi}{N} n k} = \frac{1 - a^N}{1 - a e^{-j \omega} }$$

A simple inspection of those DTFT abd DFT yields the following conclusion:
$$ \boxed{ X[k] = (1-a^N) X(e^{j\omega_k}) } $$

Where the DTFT is evaluated at the DFT frequency of $w_k = \frac{2\pi}{N}k$.

And the error between $X[k]$ and $X(e^{j \omega_k})$ is:
$$ E[k] = X(e^{j \omega_k}) - X[k] = a^N X(e^{j \omega_k}) $$

As can be seen, the error due to the finite length DFT truncation is multiplied by the constant  $a^N$ and the frequency dependent variable $X(e^{j \omega_k})$ for each DFT frequency of $\omega_k$. It's obvious from this that the error plot will have a shape of the actual DTFT $X(w)$ scaled by $a^N$.

In order to make this error smaller than the floating point round-off errors (such as $10^{-16}$ for 64 bit IEEE floating point numbers) you have to make sure that for all frequencies $\omega_k$ for $0 \leq k \leq N-1$ you should satisfy $E[k] < 10^{-16}$ that is:$$ E[k]=a^N X(e^{j \omega_k}) < 10^{-16}$$

Note that this error is frequency dependent due to $X(e^{j \omega_k})$ term: the higher the DTFT at that frequency, the higher will be the error and vice versa. So for a randomly fluctuating DTFT $X(e^{j \omega_k})$ the error will also randomly fluctuate.

In the analysis of the Gaussian signal $x[n] = e^{-\alpha n^2} = a^{n^2}$ where $ 0 < a = e^{-\alpha} < 1$ , an analytic expression for the error will be difficult to simplify as you have to evaluate the following sum:

$$ X_g(e^{j\omega}) = \sum_{n=-\infty}^{\infty} a^{n^2} e^{-j \omega n} = \beta a^{\gamma w^2}$$ where $\beta$ and $\gamma$ should be determined.

And the N-point DFT of the truncated (and shifted by $N/2$) Gaussian signal:
$$ X_g[k] = \sum_{n=0}^{N-1} a^{(n-N/2)^2} e^{-j \frac{2\pi}{N} n k} = ? $$

Unfortunately since the series sum may not have a closed form expression (or it will be quite complicated) therefore a simple formulation of the error does not exist as were the case for the simpler signal $x[n] = a^n u[n]$.

EDIT: The following exposition shows the same results as attained above however from a different perspective which may be less intuitive to follow. I suggest the above perspective for a clear understanding and the below perspective a complementary case.

Consider the following infinite length signal $x[n] = a^n u[n]$, with $|a|<1$, where the discrete-time Fourier tranform (DTFT) $X(e^{j\omega})$ of $x[n]$ is:

$$X(e^{j\omega}) = \text{DTFT} \{x[n] \} = \sum_{n=0}^{\infty} a^n e^{-j \omega n} = \frac{1}{1-a e^{-j\omega} } $$

In order to practically compute the DTFT $X(e^{j\omega})$ using a digital computer, the relative transform, DFT (discrete Fourier transform) is used. For finite length signals it's well known that DFT $X[k]$ will be samples of DTFT $X(e^{j \omega})$. But for infinite length signals (such as $x[n]=a^n u[n]$) the computed N-point DFT (which should practically use a finite length truncated version of the infinite length signal) will be an approximation to the true theoretical DTFT. The question is ; what's the relation between the true DTFT spectrum $X(e^{j\omega})$ and the computed DFT samples  $X[k]$ ?

Now it's obvious that for infinite length signals the values of DFT $X[k]$ will not be the samples of the function $X(e^{j\omega})$. For example given $x[n]=a^nu[n]$ and $X(e^{j \omega}) = \frac{1}{1-a e^{-j\omega}}$ then we see that  $$X[k] \neq X(e^{j\frac{2 \pi}{N}k}) = \frac{1}{1-a e^{-j\frac{2 \pi}{N}k}}$$

First lets consider the following, assume you have artificially obtained $N$ uniform samples from the function $X(e^{j\omega})=\frac{1}{1-a e^{-j\frac{2 \pi}{N}k}}$ as
$$ X_f[k] = \frac{1}{1-a e^{-j\frac{2 \pi}{N}k}} $$ where $X_f[k]$ stands for the  DFT samples of a finite length signal $\tilde{y}[n]$ whose N-point DFT would be $X_f[k]$. Therefore we can assert that the inverse DFT of the $X_f[k]$ would yield the N-point discrete-time sequence $\tilde{y}[n]$ as:

$$ \tilde{y}[n] = IDFT\{ X_f[k] \} = \frac{1}{N} \sum_{k=0}^{N-1}  \frac{1}{1-a e^{-j\frac{2 \pi}{N}k}} e^{ j \frac{2\pi}{N} k n} $$

It can be shown from the periodicity assumption of the IDFT that the result is also equivalent to:
$$ \tilde{y}[n] = IDFT\{ X_f[k] \} = \sum_{r=-\infty}^{\infty}  x[n-rN] $$

Where the signal $x[n]$ is nothing but the theoretical infinite length inverse of the DTFT $X(e^{j\omega}) =  \frac{1}{1-ae^{-j\omega}}$ which is as we know
$$ x[n] = IDTFT \{\frac{1}{1-ae^{-j\omega}} \} = \frac{1}{2\pi} \int_{-\pi}^{\pi} \frac{1}{1-a e^{-j\omega}} e^{j\omega n} d\omega = a^n u[n] $$

Therefore we have the following:

$$\tilde{y}[n] = \sum_{r=-\infty}^{\infty}  x[n-rN] = \sum_{r=-\infty}^{\infty}  a^{n-rN}u[n-rN] $$

It will be interesting to consider the value of $\tilde{y}[n]$ for the range $0 \leq n \leq N-1$ for which we have:

$$\tilde{y}[n] = a^n \sum_{r=-\infty}^{\infty}  a^{-rN}u[n-rN] = a^n \sum_{r=-\infty}^{0}  a^{-rN} = a^n \sum_{r=0}^{\infty}  (a^N)^r = a^n \frac{1}{1-a^N} $$

Which shows the relation between finite length signal $\tilde{y }[n]$ and the infinite length signal $x[n] = a^n u[n]$ as:
$$\tilde{y}[n] = \frac{x[n]}{1-a^N}$$ where $x[n] = a^n u[n]$.

One can see that the modification involved on the computed value of finite length approximation to the infinite length signal $x[n]$ is a constant weight of $K = \frac{1}{1-a^N}$ for the signal $x[n]= a^n u[n]$

Now using the linearity properties of DFT one can deduce the following:
Given the pair of N-point sequence $\tilde{y}[n]$ and its N-point DFT $X_f[k]$ which is the uniform samples of the true theoretical DTFT $X(e^{j\omega})$ associated with infinite length signal $x[n]=a^n u[n]$ we have the following

$$\tilde{y}[n] \leftrightarrow X_f[k] $$
$$\tilde{y}[n] = K x[n] \leftrightarrow X_f[k] \Rightarrow x[n] \leftrightarrow \frac{1}{K} X_f[k] = (1-a^N) X_f[k]$$

This last statment is concluded as follows: given an N-point finite length truncation $x_N[n]$ of the infinite length signal $x[n] = a^n u[n]$ for the range $0 \leq n< N$, then the N-point DFT of the N-point sequence $x_N[n]$ will be $(1-a^N)X[k]$ where $X[k]$ are the samples of the true DTFT spectrum $X(e^{j\omega}) = \frac{1}{1 -a e^{-j \omega}}$ associated with the infinite length signal $x[n]=a^n u[n]$. So the error associated with the DFT computation is the fixed weight $1-a^N$ per sample of the DFT.

You should perform a similar analysis for the Gaussian signal $x[n] = \beta e^{-\alpha n^2}$ to find the error associated per computed DFT sample.


Calculate signal to noise ratio from the Hillbert envelope of a stacked signal I am working with seismic waveforms and more specifically with a stack of years of daily seismograms, that I would call the reference waveform. I am trying to calculate the SNR for each daily waveforms by using the envelope of the reference, but I have doubts on how I should calculate the noise.
Concretely, I am trying to reproduce what was done in this paper (page 8).

I believe It should be simple, but I have poor basics in signal processing and struggle with notations.

Based on this section of the paper:

Here is my attempt at doing a calculation like that.

The plot above shows three plots:

The average of the dailies (in green) and all the daily data (in black).
The function $\sigma$ calculated from the dailies.
The SNR calculated from $\sigma$ and the dailies.

Scilab code for what I did is below.

// 25624
Tdata = 33;
stack_average = wfir("hp",Tdata,[.3 0],"hm",[0 0]);

Ndays = 30;

ccs = zeros(Ndays,Tdata);

for k=1:Ndays,
    ccs(k,:) =  stack_average + 0.05*rand(1,Tdata,'norm').*window('hm', Tdata);
end;

clf
subplot(311)
plot(ccs','k')
plot(stack_average,'g')
xtitle('$\mbox{All } cc \mbox{ in black and } <cc> \mbox{ in green}$')

sigma_N_t = sqrt((mean(ccs.^2,'r') - mean(ccs,'r').^2)/(Ndays-1));
s_N_t = abs(hilbert(mean(ccs,'r')));
SNR_N_t = s_N_t ./ sigma_N_t;

subplot(312)
plot(sigma_N_t)
xtitle('$\mbox{\rm Level of noise,}\ \sigma$')

subplot(313)
plot(s_N_t, 'r')
xtitle('$\mbox{\rm Signal to noise ratio,}\ SNR$')


Does Gamma Correction Interfere with Sub Pixel Motion Estimation? Advanced Video Coding (AVC) and High Efficiency Video Coding (HEVC) work on gamma corrected, non linear, colors.

When interpolating the value at subpixels, a multi tap algorithm processes the gamma corrected values of the original pixels directly, without converting into linear space and back, like a GPU does when doing texture mapping.

Isn't this wrong?

When one does an operation Wrong / Right aren't strictly defined.
In most cases the questions are:

What's the model?
Is the model reasonable?

When you do sub pixel motion estimation on frames which are highly correlated you can assumes the change in values, even after non linear operation like Gamma Function, is small.

So basically the model assumes that even after Gamma Correction the spatial correlation between 2 frames is kept.

You can do a small test by generating image supposedly in Linear Color Space.
Generate another image as a result of a small movement of the previous image.
Have a look on the spatial correlation between the 2.
Now transform each using Gamma Function and measure the correlation as well.

If their model is reasonable, the correlation in both cases will be similar.


MVDR Beamformer, minimize but do not change pattern in main beam Is there a way to minimize all signals that are not coming from the mainbeam?

I would like to keep the main beam, or at least +-3° around the main peak stable.

When I estimate the covariance from the data itself, it always also minimizes my wanted signal in the main beam. I can not estimate the covariance in the absence of the signal.

Any advice?

Use a white noise gain constraint on the filter coefficients, which is a quadratic inequality constraint. When you have sources in the main beam that are not exactly on the unity gain constraint, the filter coefficients get large.  You can control the magnitude by diagonal loading the covariance matrix. There is a full derivation in "Detection, Estimation, and Modulation Theory, Volume 4, Array Processing" by Harry VanTrees.  It solves the problem but I'm not aware of a nice online algorithm. There is a paper by Cox and Zeskind that has an online algorithm but it isn't optimal no matter what Harry Cox thinks.

Another approach people use is to use a direction finding algorithm prior (before) to set the unity gain constraints (actually LCMV ) and then beam form. I prefer the white noise gain constraint.

Bob Zeskind suggested using an unmodified MVDR but forming about 3 times as many beams as you would in a conventional beamformer in a conversation I had with him. It's not elegant but does seem to have uses.


Frequency spectrum of a set of images in Matlab I wish the compute the frequency distribution of a set of images. I basically am keen on finding out how the the frequency varies in the images.

I know I can use the 'fft2' function to compute the 2-dimensional Fourier Transform. I can do this for each image. But I am not sure how to obtain the isotropic frequency distribution of all the images. I want to find a method to compare the frequency distribution.

What is the best way to go about this problem? Thanks for your help.

If I'm interpreting the question correctly, you're looking for a way to find a distribution of the frequencies in each images. It sounds like you might need some sort of 3D scatterplot.

If you have 20 images that are 256x256, how I imagine this working is,

Take your 2D FFT or 2D DCT (or blockwise FFT/DCT) for each image
Normalize all the values so you can have a good comparison between each image
Concatenate all the frequency domain images into a 256 x 256 x 20 matrix
Filter out the values that are below a certain threshold because they do not represent a significant component of your image. For example, if the magnitude of the values in the corners that represent high frequency information are very low, then removing them won't make a huge difference to your image because most of the information is elsewhere. But removing those data points create holes in your 3D plot so you can see the information that you're interested in
Plot the data in 3D. If you're using Matlab, you can use scatter3(). The plot will be easier to read if you apply a heatmap.

Otherwise, if you're looking for the frequency information of a single image, there's several ways to represent the data depending on your problem. If you clarify the question you have and it's application, I can edit this answer.

Histogram

If you take the histogram of a grayscale image, you'll see a distribution of the intensity values of each pixel.

This easiest way to see this is to open up the image in a photo-editing software, like Photoshop, GIMP, Picasa and select histogram from the menu.

If you want to use Matlab, the hist() function plots a histogram similar to the one below. You can directly pass your image variable into the hist() function. But it's important to pick a good bin size, because different bin sizes can dramatically change the representation of the data. There are ways to calculate optimal bin sizes, but you can always guess and check.

2D FFT

Follow the Wikipedia article for the theoretical explanation, but basically the 2D FFT gives you an image of the intensity of low frequency to high frequency information. The center of the image is low frequency information (such as flat surface, sky, skin, walls) and the corners are high frequency information (edges, grainy noise, intricate patterns).

If you're using Matlab, you need to shift the 2D FFT to put the low frequency information in the center of the image. It's also a good idea to normalize the 2D FFT so you can see the details better; it's similar to contrast enhancing a photograph to see the details.

2D DCT

Follow the Wikipedia article for the theoretical explanation. The 2D DCT is very similar to the 2D FFT; It also gives you an image of the intensity of low frequency to high frequency information, but the low frequency information is in the top-left corner and the high frequency information is in the bottom-right corner. This makes the information easier to manipulate (at least for me).

Generally, the 2D DCT is used on blocks of an image (8x8 pixels). This gives you frequency information that is still related to its location in the spatial image.

The 2D DCT is useful for encoding images. A variant of JPEG encoding uses the 2D DCT transform.


Texture detection and region properties on an image I have an interesting problem that I'm trying to solve. So for example, if we have a black and white image containing two textures (A and B), I'm interested in knowing the pixel values of the boundaries that would completely cover the individual textures.

I thought of using cross correlation, which would give me a cluster of the positions, of the template (legends) on the image but is there a way to use that to get the pixel values of the boundaries (these are typically irregular)? Also, is there any other better way of doing this?

A real world example of this would be in detecting the regions of rainfall on a map showing three levels of rainfall in a year. Each level is assigned a texture thats present in the legend, which is used for texture matching on the same image.

Here's a simple process:

Assign texture measurements to each region in the image.
Use a region labeling / connectivity algorithm (or region growing algorithm) to join adjacent regions having the same texture measure.
Implement a simple edge-following algorithm to trace the contour of each region.

Law's Texture Measures are an older but still useful technique for determining the texture in an image, and they may be sufficient for you to distinguish texture A from texture B throughout the image. See the section "Laws Texture Energy Measures" in the Wikipedia article:

http://en.wikipedia.org/wiki/Image_texture

As a first step, calculate all texture measures and determine which particular measure (e.g. Edge or Spot) allows you to distinguish one texture from the other most easily. (If you post some pictures I could help you identify a texture measure.)

If you have only two textures, A and B, then you can treat them as foreground and background, and a standard region labeling algorithm will work. To make it easier to see what's happening in the processing, you could generate a new image by assigning A texels (texture elements, little chunks of texture) to the color white, and B texels to the color black. Region labeling and/or contour-following algorithms would then find the connected white and black regions. The findContours() function in OpenCV will work well.

http://en.wikipedia.org/wiki/Connected-component_labeling

The same Wikipedia article includes both the traditional multipass algorithm as well as a single-pass algorithm. I haven't implemented the single-pass algorithm described there, but I have worked with the single-pass algorithm described in the paper "A Component-Labeling Algorithm Using Contour Tracing Technique" by Chen and Chang. The paper by Chen and Chang also describes a standard contour-following algorithm that can be implemented quickly.

If you have more than two textures, then you might use a watershed or mean shift algorithm to cluster regions together after you've remapped textures to colors. Although this remapping from texture to color isn't necessary, it does make the process easier to debug and understand.


Super basic questions on statistical process Before starting: I am really a beginner in statistical process in time. I mainly do quantum information and while learning aspect of quantum noise I realized that I am actually too weak on basics of statistical process.

Thus, I would like to first tell you what I basically understand and that you confirm or correct me. Then I have some more specific questions.

Allright, here is what (I think) I understood, and I would like to check

Let's consider a random process. I consider a random variable $X$.
First, I can fix a given time $t$, and ask questions about the statistical properties of $X(t)$. For example its mean value $<X(t)>$, its standard deviation $<(X(t)-<X(t)>)^2>$, and I could imagine many other quantities. In principle without further assumptions I need a lot of information (an infinite one ?) to specify all the characteristics of the statistics of $X(t)$.

Furthermore: here it is what is happening at a given time $t$. Even if I fully described the statistics of $X(t)$ I would lack information. For example there might have some correlations between what happens at $t$ and what happens at $t'$. Thus I also need to specify some properties linking those different times. They are not given by what happens at a given time $t$.

An example I have in mind is a variable that verifies: $X(t) \in [0,1]$, uniformously distributed. But the correlations are such that $X(t'>t)=X(t)$. Then if I do a first experiment I will have $X(t)=0.2$ for example. And for all further times $X(t')=0.2$ as well because of correlations. If I re-run the experiment I can have $X(t)=0.87$ for example but again the values for further times will be the same as $X(t)$. Do you confirm my understanding ?

To specify what happens in time we thus need other information like autocorrelation functions: $C(t,t')=<X(t)X(t')>$. And probably many other quantities that I am not aware of.

More specific questions

What is a rigorous definition of stationnary statistical process ? What I basically understand is that the statistical properties of $X(t)$ actually do not depend on $t$. And the statistical properties of $X(t)X(t')$ only depend on $t-t'$. But I guess there is a more precise definition than this.

In a paper I am reading, they assume a gaussian, stationnary process. They say that knowing the average value $<X>$ and the autocorrelation function allows to fully specify the process. I do not understand this. I agree that it will fully describe the statistical behavior at a given time (because from autocorrelation function we have access to the variance in addition to the mean we already knew). But how is the autocorrelation function the only information we need to specify the statistical relationship between different times ?

Finally, basic question but... what is the name of the precise field that deals with all this ? It is not statistics (this is way too broad), I tried some keyword with "statistical noise" but I didn't find something good. Related: is there a good reference to really give the basics ? I don't think I need much.

Answer to your question related to Gaussian Random Process :

  For a Gaussian Random Process, Wide-Sense Stationarity implies Strict-Sense Stationarity as well.

That means, if mean $\mu_X(t)$ is independent of $t$ and is a constant and Auto-correlation function $R_{X}(t, t+\tau)$ is dependent only on time-shift $\tau$, then the Gaussian Process is actually Strict-Sense Stationary.

Why?
Because a Gaussian Random Process, $X(t)$, is defined as that random process for which if we freeze time at $n$ instances to get random variables $\{X(t_1), X(t_2), ..., X(t_n)\}$, then these random variables are jointly Gaussian Distributed. So, we can say that if we have a Gaussian Random Process $X(t)$ and we freeze time at $n$ instances, $t_1, t_2, ..., t_{n}$, then we get a Jointly Gaussian Random vector $\mathbf X$ of $n$ dimension. And for a jointly Gaussian Random Vector, its joint Probability Density is fully determined by its Mean vector $\mathbf \mu_{\mathbf X}$ and Covariance Matrix $\mathbf C$, given by:

$$f_{\mathbf X}(\mathbf x) = \frac{1}{(2\pi)^{\frac{n}{2}}|C|^{\frac{1}{2}}}\cdot \exp(-\frac{1}{2}(\mathbf x - \mu_{\mathbf X})^TC^{-1}(\mathbf x - \mu_{\mathbf X}))$$

Also, note that the Auto-Correlation and Covariance Matrix $\mathbf C$ are related for a Wide-Sense Stationary Process, that, all the elements $\mathbf C_{i,j}$ can be computed from the Auto-Correlation function and constant mean of the WSS Random Process. $$\mathbf C_{i,j} = \mathbf E[(X(t_i) - \mu_X)(X(t_j)-\mu_X)] = \mathbf R_X(t_i-t_j) - \mu_X^2$$
Thus, if the mean of a Gaussian Random process is constant $\mu_{X}$ and Auto-Correlation is only dependent on time-shift $\tau$, then the Gaussian Random Process is Strict Sense Stationary, because we know all the joint Distributions completely for any $n$.This is not true for any other Random Process in general.

For detailed definition of Strict-Sense Stationary Random processes and Wide-Sense Stationary Random Process, see the rest of the explanation.

The way I visualize Random processes $X(t, \zeta)$ is: for every outcome $\zeta$ in a random experiment, the outcome $\zeta$ is mapped to a function $x(t, \zeta)$. This is like Random Variable concept. A random variable $X$ is a rule for assigning a value or set of values $x(\zeta)$ to each outcome $\zeta$ of a random experiment. Similarly, a Random process $X(t)$ is a rule of assigning a function $x(t,\zeta)$ to each outcome $\zeta$ of a random experiment. Three important intuitions are as follows:

If we freeze the outcome $\zeta$ of that random experiment, we will get a deterministic function of $t$, $x(t)$.
If we freeze $t$, then we get a Random Variable $X(\zeta)$.
If we freeze both $t$ and $\zeta$, we get a number $x$.
$$X(t, \zeta) : ensemble \ of \ function \ of \ t \ which \ are \ uncertain$$

Consider a Random process $X(t)$ (dropped $\zeta$ to de-clutter) $\forall \ t \in \mathbb R$.

$X(t)$ will be called Strict-Sense Stationary(SSS) Process if its statistical properties are time-shift invariant. This means that the random processes $X(t)$ and $X(t+\tau)$ have the same statistics for all $\tau \in \mathbb R$.

What I mean by statistical properties is: If we freeze time at $t$, and get random variable $x(t)$ from the random process $X(t)$, then we can define $f(x(t))$ as the probability density function of $x(t)$. Similarly, we can define joint probability density function of $n$ such random variables $x(t_1), x(t_2),..., x(t_n)$ by fixing time at $t_1, t_2,...,t_n$, as $f(x(t_1), x(t_2),...,x(t_n))$. Then, Strict-Sense Stationary says that:
$$f(x(t_1), x(t_2),...,x(t_n)) = f(x(t_1+\tau), x(t_2+\tau),...,x(t_n+\tau)), \ \forall \ \tau \in \mathbb R$$

We can put $n=1$ in the above expression to figure that, $1^{st}$ order probability density of SSS Random processes are independent of time $t$:
$$f(x(t)) = f(x(t+\tau)) = f(x)$$

Similarly, if we put $n=2$, and look at the $2^{nd}$ order probability density $f(x(t_1), x(t_2))$ will also depend only on time-shift $\tau$ between $t_2$ and $t_1$. Since we can choose any value of $\tau$, consider the case that we choose $\tau = -t_2$ always, then:
$$f(x(t_1), x(t_2)) = f(x(t_1-t_2), x(0))$$
The above expression says that joint density function of $2$ random variables $x(t_1)$ and $x(t_2)$ which we got by freezing time at $t_1$ and $t_2$, is always equal to the joint density function of $2$ random variables $x(0)$ and $x(t_1-t_2)$ which we get when we freeze time at $0$ and $(t_1-t_2)$.
This means that $2^{nd}$ order density function is independent of time and only dependent on the time shift $(t_1-t_2)$.

Wide-Sense Stationary(WSS) Random Processes are those Random Processes which satify following $2$ requirements:

Mean of the Random Process, $\mathbf E[X(t)]$, is independent of time and is a constant.
$$\mathbf E[X(t)] = \mu_X = c$$
Auto-correlation function of the Random Process, $\mathbf R_X$, is only dependent on time-difference $\tau = t_1 - t_2$
$$\mathbf E[x(t_1)x^*(t_2)] = \mathbf R_X(\tau)$$

For a WSS Process $X(t)$, the above requirement implies that the variance of the Random Process, $\mathbf E[(X(t)-c)^2]$, is also a constant.
$$\mathbf E[(X(t)-c)^2] = \mathbf E[(X(t))^2] - c^2 = \mathbf R_X(0) - c^2$$

Note that there is no requirement on statistical properties of Wide Sense Stationary Random Processes.

The field of Stochastic Processes deals with these topics and a good starting reference would be Chapter: 9 of "Probability, Random Variables and Stochastic Processes, by A. Papoulis & S. Unnikrishna Pillai"


problem with Fourier transform I am reading some paper, and I am having some trouble with some Fourier transform,

Suppose that $F(\omega)$ is the Fourier transform of $f(x)$, i.e. where

$$F(\omega)=\int_{-\infty}^{\infty}f(x)e^{-i\omega x}\,dx.$$

What is $$\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{F(\omega-a)}{i\omega}e^{i\omega x}\,d\omega ?$$
I was thinking of using the time integration property, since $i\omega$ appears at the denominator, but I am not sure how to start it.

Since

$$f(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)e^{i\omega x}d\omega$$

you are looking for the inverse transform of

$$G(\omega)=\frac{F(\omega-a)}{i\omega}$$

If you look at a table of Fourier transform correspondences you will find

$$\frac{F(\omega)}{i\omega}\Longleftrightarrow \int_{-\infty}^x f(y)dy$$

and

$$F(\omega-a)\Longleftrightarrow e^{iax}f(x)$$

Consequently, the function whose Fourier transform is $G(\omega)$ is the integral of a modulated version of $f(x)$:

$$G(\omega)\Longleftrightarrow \int_{-\infty}^x e^{iay}f(y)dy$$

EDIT:

The correspondence

$$\frac{F(\omega)}{i\omega}\Longleftrightarrow \int_{-\infty}^x f(y)dy\tag{1}$$

assumes that

$$F(0)=\int_{-\infty}^{\infty} f(x)dx=0$$

This assumption is often made implicitly because it is a standard requirement (obviously, a sufficient one) for the Fourier transform to exist. If

$$g(x)=\int_{-\infty}^x f(y)dy$$

it simply means that $g(x)$ is a well-behaved function satisfying $\lim_{x\rightarrow\infty}g(x)=0$. However, it is also possible to compute the Fourier transform of functions that do not satisfy this condition. In this case, the correspondence (1) has to be modified:

$$\frac{F(\omega)}{i\omega}+\pi F(0)\delta(\omega)\Longleftrightarrow \int_{-\infty}^x f(y)dy\tag{2}$$

which obviously simplifies to (1) if $F(0)=0$. If this condition is not satisfied, the result for the inverse transform of $G(\omega)$ given above must be modified by an additive constant:

$$G(\omega)\Longleftrightarrow \int_{-\infty}^x e^{iay}f(y)dy-\frac{1}{2}\int_{-\infty}^{\infty}e^{iax}f(x)dx$$


DFT as an Orthogonal Basis Change In one of the homeworks that I am dealing with for Linear Systems course, I have encountered with such a statement: Consider $\mathbb{C}^N$ the vector space of N dimensional complex vectors. We can define a basis
$F=\{f_1,\ldots,f_N\}$ where
\begin{align}
f_k = \begin{bmatrix}
    f_{k,1}       \\
    f_{k,2}        \\
    .\\
    .\\
    f_{k,N}
\end{bmatrix}
\end{align}

and $f_{k,l} = \frac{1}{N}e^{\frac{j2\pi(k-1)(l-1)}{N}}$. It is straightforward to show that those vectors are orthogonal, but I have no idea about how to show that those vectors are linearly independent, and their span is $\mathbb{C}^N$. Could you please give me a clue?

Hint 1: you can rewrite products of the exponential argument:

$$ f_{k,l}=\frac{1}{N}\left(e^{\frac{j2\pi(k-1)}{N}}\right)^{l-1}$$

and recognize a Vandermonde matrix.

Hint 2: consider the dimension of the subspace spanned, and the linear independence.


Can an infinite bandwidth signal be an infinite-time signal? While decomposing the impulse train using the Fourier series the number of the samples are infinite hence representing a infinite length of train.

But while studying the Nyquist Criteria in the frequency domain it is shown
the expression to have $m$ samples where  $-\infty<m<\infty$

While studying the concept of the Nyquist Criteria in the ISI model where for the system to be interference free the value of the sampling instant must be given as following:

$$h\left(nT_s\right)=\begin{cases}1;&n=0\\0;&n\neq0\end{cases}$$

For more reference here. I cannot understand that the signal is having infinite duration and infnite bandwidth accordingly is it even possible.

Let's consider the continuous time case.

Take a Gaussian $f(t) = \frac{1}{\sqrt{2 \pi}} e^{-t^2/2}$, which is positive for all time.

The Fourier transform of a Gaussian is another Gaussian: $F(\omega)= c e^{- d \omega^2 }$ for some constants $c,d$  that you can work out on your own, so a Gaussian occupies all time and all frequency. That is, you can be both infinite in time and infinite in frequency.

Your trade-off between resolution in time and frequency is given by the uncertainty principle. The Gaussian attains equality in the uncertainty principle.

It seems like you're somewhat confused as to what Nyquist's sampling theorem says -- see this page for details (or your favorite signal processing text). You're assuming the signal is band-limited (i.e. is non-zero in some finite interval of frequencies) in order to be able to sample the signal at some prescribed rate in order to be able to perfectly reconstruct the signal.


Light Detection and Ranging (LIDAR) Range Bias Compensation with Light Dependent Resistor (LDR) I have a distance sensor (VL53L0X), which measures the distance to an object. It is placed outside. I currently do 100 measurements in the best available range mode, then I take the median and insert it to a database.

The current distance is about 1m.

This is how it looks like during a period af about 12 hours. As can be seen there is a light dependent error in the data.

I also have a LDR (light dependent resistor) output (log scale on y-axis), which shows this more. At "peak" time the error is about 15cm.

Green: distance in meters, yellow, Voltage from LDR (light dependent resistor). The darker it gets, the higher the voltage of the LDR. At night it's about 2.5V.

Now, the distance hasn't changed in those 12 hours. My question is: how can I compensate this error? I am looking for a term on how this procedure is called.

My first guess was to calculate the difference to 1m for all distance points, then maybe fit the LDR output to that error points.

Data: https://pastebin.com/2fXcSUQB

Column 1: distance in meters, Column 2: voltage from LDR

Ideally, you would like to explain the LIDAR signal as a combination of one or more other variables, here LDR and/or Temperature and then remove the effect of those variables from the LIDAR signal.

The first part can be achieved with linear (hopefully) regression and specifically Multiple Linear Regression where you would be fitting a model like LIDAR = Constant + LDR_Coef * LDR + Temp_Coef * Temp.

With the availability of just the LDR data (code in GNU Octave but easily transferable to other platforms too):

%Load the data
D = dlmread("myfile.txt",";"); %myfile.txt is the file you provided which is delimited by semicolon.
% Adjust the LDR signal
D(:,3) = log(D(:,3));
% Fit the model using regress.
% For more information see: https://octave.sourceforge.io/statistics/function/regress.html
[B, BINT, R, RINT, STATS] = regress(D(:,2), [ones(size(D(:,3))), D(:,3)]);

After you do this, you get your model coefficients in B and a few other things (more about this later on).

It is now possible to express LIDAR as a function of LDR:

% Plot LIDAR and LIDAR as a function of LDR
r=0:(length(D)-1);
plot(r,D(:,2),B(1)+B(2).*D(:,3));
xlabel('Time');ylabel('LIDAR / LDR');legend('LIDAR','LDR');grid on;

This now looks like:

This concludes with the fitting part. Now we want to remove the variation due to LDR, which is the "easy" part:

LIDAR_adjusted = D(:,2) - B(2).*D(:,3);
plot(r,D(:,2),LIDAR_adjusted);
xlabel('Time');ylabel('LIDAR');legend('LIDAR','LIDAR_adjusted');grid on;

Which results in:

The STATS vector contains a number of coefficients that inform you about how good is this fit. In this case 7.6591e-01   1.4952e+03   0.0000e+00   7.8229e-04. So, the $R^2$ is not ideal but not too bad too (that is the first number, it can be between $0$ and $1$ and it is ideal at $1$), the rest of the numbers are not bad too (the error is very small), but if you look at the residuals (that is the difference between the LIDAR and the LIDAR as a function of LDR) it is not exactly noise. Here is the residual (R):

plot(R);grid on;

Which results in:

There are clearly trends in the residual which means that LDR on its own cannot explain the variation in the LIDAR signal. You can add 'r' as a regressor in there, which would return three coefficients for the model LIDAR = B(1) + B(2)*D(:,3) + B(3)*r and that would "take away" that linear upwards trend as a function of time and improve the error but you still get those two humps because the sensor responses are non-linear.

To conclude with this, if you wanted to add the temperature to the model too, just change the regress line to:

[B, BINT, R, RINT, STATS] = regress(D(:,2), [ones(size(D(:,3)), D(:,3), D(:,4)]);

Assuming that you will add another column for the temperature. The rest of the process remains the same (with reasonable adjustments for removing the components further below).

The other thing that you can do of course is a Discrete Fourier Transform (just on the LIDAR) and retain the first few harmonics from it for detrending. The latter assumes that there is a periodic well defined effect of light (and / or temperature and other parameters) on LIDAR which can be approximated with a set of periodic components. And in that case, the adjustment comes as a function of the time of day "for free".

Hope this helps.


Using the flip-and-drag method to draw a convolution So, we are asked to graphically convolve (using the flip-and-drag method) the following pair of signals: $$x_1(t) = \delta(t+1)-2\delta(t-1) \ \ \ \ \ \ \ \ \ \ \ \ \ x_2(t) =  \begin{cases} -1-2t, & -\frac{1}{2}\leq t\leq 0 \\ -1+2t, & 0\leq t\leq \frac{1}{2} \\ 0 & |t|>\frac{1}{2}\end{cases}$$

where $\delta(t)$ is the dirac delta function.

I don't know how to draw using Latex, so I am going to try and describe you what I've drawn so far using words (as best as I can). So, I've successfully been able to draw $x_2(t)$ (which gives me an inverted triangle whose summits coordinates are : $\left(-\frac{1}{2}, 0\right), \left(\frac{1}{2}, 0 \right)$ and $\left(0, -1 \right)$).

My problem starts when I have to draw $x_1(t)$. After doing some research, I have found that we represent the dirac delta function using "arrows", but I still don't know how to draw $x_1(t)$ .

Also, if we assume I was able to draw it (i.e I would have a graph with arrows), how am I supposed to draw the convolution (using the flip and drag method)?

Convolution is a linear operator, so let me explain just how to convolve $x_2(t)$ with $\delta(t-t_0)$. The only thing you would have to do afterwards is to convolve the same signal $x_2(t)$ with the other Dirac delta, multiply the function you get by some amplitude, and then add both convolution results.

The Dirac delta has what is known as sifting property:

$$\int_{-\infty}^\infty f(t)\delta(t-t_0)\,\mathrm{d}t=f(t_0)$$

This means that if we multiply a function by an impulse located at $t=t_0$, the area under that differential along the $t$-axis corresponds to the value of the function at $t_0$.

To solve your convolution graphically, imagine the following. Leave the impulses where they are. Then flip your convoluting signal, $x_2(t)$, and take it to start at $-\infty$. Start sliding this function to the right. At first, the convolution (the integral) returns $0$ for all $t$ because the product of $x_2(t)$ and $\delta(t+1)$ is $0$, as they are not intersecting anywhere. But, at some point, they both come across each other. This will continue to happen as long as the sliding version of $x_2(t)$ takes some non-zero value in $t=t_0$. When the triangle goes through entirely, then the convolution will be zero for every other $t$ again, as they will not encounter again.

I'll put some (extremely rudimentary) drawings to make this easier to see.

The first picture corresponds to the cases where the triangle and the impulse don't take any non-zero values at common points. The second one, where the triangle and the impulse are both non-zero at some value (in this case, always at $t_0$). The third picture corresponds to the moment where the triangle has passed through entirely, and thus the convolution is $0$ again.

Try to think about it, and do it on your own. As a rule of thumb, convoluting with deltas is by far the easiest convolution that there can be. In general, if we have something like this (with $A$ being some constant that represents the amplitude of the impulse): $$g(t)=f(t)*A\delta(t-t_0)$$ Then $g(t)$ will just be $f(t)$ but centered at $t_0$ and multiplied by $A$. That means, you just shift $f(t)$ and invert it or scale it depending on the value of $A$.


Methods for aligning FFT bins with an FSK signal Say I want to demodulate some FSK modulation. My (novice) understanding is a pretty good way to do this is by taking the FFT, with the length chosen such that the bin spacing is equal to the tone spacing. In the common case of the tones being spaced at the reciprocal of the symbol duration, this just means taking the FFT over the symbol period. Symbol detection could then be as simple as looking for the bin with the highest magnitude.

Let's assume timing issues are already solved: there's still the problem of frequency alignment. If the receiver's frequency is a little bit off then some of the tone power will bleed into adjacent bins, degrading accuracy.

What are some basic algorithms for addressing this?

Before I address your question, I'll suggest that using a frequency domain technique (FFTs) to demodulate inherently time domain data (a time series of sequential data symbols) won't be very practical.  Symbol timing recovery is inherently a time domain process. ISI encountered in common modulations like GFSK/GMSK means you really have more than 2 frequencies for 2FSK symbols at the receiver.

Now on to your question of how to perform frequency alignment.

You requested algorithms for what is normally called Carrier Frequency Offset (CFO) estimation and correction.

In communications receivers, there are sometimes two stages of CFO estimation and correction: coarse and fine.  For 2-FSK burst communications, a single stage of coarse CFO correction is usually all one would need for non-coherent FSK demodulation.

CFO correction is generally performed as a (cyclic) digital spectrum shift by multiplication with a time-varying complex exponential in the time domain:

$$x_c[n] = x[n]\space e^{j\pi \frac{-f_{offset}}{F_s/2}n}$$

where $f_{offset}$ is the estimated CFO, and $F_s$ is the sample rate.

Although for FSK, if one has already computed the baseband waveform from the FSK-modulated symbols (i.e. the output of an FM discriminator), a simple level shift of the basedband waveform, to have it's peaks and troughs centered about 0, is equivalent to CFO correction.

CFO estimation algorithms can be divided up into two broad classes: Open-loop feed-forward algorithms and closed-loop feedback algorithms.  Closed-loop algorithms work a sample at a time, are usually used for fine frequency offset estimation, and usually involve some arrangement of PLLs.  Open-loop algorithms work on a block of samples at a time and are used for coarse frequency offset estimation.

You probably just want to use an open loop instantaneous frequency estimate for computing your CFO estimate.  The instantaneous frequency estimate of an FM/FSK waveform at any given sample is given by

$$y[n] = \dfrac{F_s/2}{\pi}\mathrm{arg}\left(x[n]x^*[n-1]\right)$$

Normally to generate an CFO estimate, one computes the instantaneous frequency over several symbols and looks at the average peak and trough levels to estimate how much the waveform level needs to be shifted for the waveform to be centered about 0.

In your proposed case of having all the samples for one symbol, assuming you don't have a terribly large CFO, computing the instantaneous frequency for that one symbol will tell you what 2FSK symbol you have without performing the FFT.


What channel model to choose while simulating a system in passband with frequencies over GHz range? I was wondering as to what kind of channel model do we use to test for BER vs SNR in a passband RF system. The system that I designed is basically a wideband frequency hopped OFDM system where I use Digital up-converter (DUC) to up convert the signal to a frequency in GHz range and use the Digital down converter (DDC) to down convert the signal back to baseband in the receiver side. Now my system works fine. However, when I add an AWGN channel, I get pretty low BER of 10^-5 in passband as compared to around 10^-4 BER in baseband under same AWGN channel of SNR 10 dB. Why is it happening? I thought that there would be more attenuation in passband as attenuation increases with respect to frequency. Or is it that my way of modelling the channel isn't good. Btw, I have tried adding Rayleigh SISO fading channel along with AWGN but in that case I achieve around 0.4 BER which is quite high. I just don't know anything about BER vs SNR plots in passband. All the plots that I found are mostly in baseband.

So my question is

How do I model the channel in passband?
How does a typical BER vs SNR look like in passband?

Thank you!

However, when I add an AWGN channel, I get pretty low BER of 10^-5 in passband as compared to around 10^-4 BER in baseband under same AWGN channel of SNR 10 dB. Why is it happening?

You've got a bug in either passband channel or baseband channel model, or in how you transform your noise or your signal – end of story!

  I thought that there would be more attenuation in passband as attenuation increases with respect to frequency

If your passband model does that, then your baseband model has to do the same – they are mathematically equivalent.

If your passband model is frequency-selective, then your baseband model has to do the same.


Bartlett's and Welch's Method for PSD Lets say. I have an image of 100 samples, and I want to find the presence of smaller image of 24 sample using cross-correlation in the fourier domain. I use Bartlett's and Welch's method for PSD, where I split up the original 100 samples into 4 segments of 25 samples.

Now I want to find the maximum or peak response,

Do I look at all 4 segments to find the maximum response?

This is just a 1D example that shows the idea (I can't seem to get scilab to cooperate with real images; at least not yet).

The top plot in blue shows the "image".  The "patch" is indicated by the red dots.  I've selected the location of the "patch" randomly throughout the "image".

The bottom plot shows the cross-correlation between the patch and the image (in green), and the red dot indicates the actual location (start position) of the patch.

The code below does no prewhitening because the "image" is generated using white noise.  Apart from the mean-correction, that is.

// 25803
clf;
NIM = 100;
IM = rand(1,NIM,"uniform");
IM = (IM - min(IM))/(max(IM) - min(IM))*255;
subplot(211)
plot(IM)

NPT = 24;
idx = ceil(rand(1,1,"uniform")*75);
PT = IM(idx + [0:NPT-1]);
plot(idx+[0:NPT-1],PT,'r.')
title('Original and patch')
NFFT = NIM + NPT - 1;

IMFFT = fft([IM-mean(IM),zeros(1,NPT-1)]);
PTFFT = fft([PT-mean(PT),zeros(1,NIM-1)]);

subplot(212)
CORR = ifft(IMFFT.*conj(PTFFT));
plot(CORR,'g')
plot(idx,CORR(idx),'r.')
title('Cross correlation and actual patch position')


Am I handling offline FFT correctly? I need some help clarifying FFTs and what they represent. I have a buffer containing compressed audio. Due to limitations, I can't handle the full uncompressed audio but can decompress small segments at a time.

Lets say I take 10 seconds of uncompressed samples, I would need to store this in  PCMaudiobuffer of size 10 * 44100 * sizeof(float) and I'd have 441000 samples. I then loop through these 441000 samples with a length of 1024 (with a 512 overlap).

for(int i =0; i < 441000; i += 512)
{
   //code to copy (i + 1024) from PCMaudiobuffer to tempbuffer. (not shown)

   //perform FFT
   InPlaceFFT(tempbuffer);

   //copy the first 512 values of the tempbuffer across to
   //a new buffer so tempbuffer can be reused in the loop
   CopyArray(512, finalbuffer);
}

Assuming all of this is correct, this is where I'm confused. Am I right in copying only the first 512 values (N/2) to the finalbuffer? (I'm using the Accelerate framework's FFT method). Is this an efficient/correct way of performing FFTs offline? Lastly, I want to be able to create a frequency spectrum from finalbuffer. To do this do I simply loop through every 44100 (1 second) and calculate the magnitude, or is 1 second generally too inaccurate for a frequency spectrum?

I may be misunderstanding all of this, so feel free to tell me to throw all this out the window :) Thanks for any help, this is quite confusing!

Edit
The Apple docs on this aren't easy to follow, but I'm using the method outlined here in vDSP_fft_zrip. (The formatting of the Apple site sometimes messes up the anchor position, so you may have to scroll down slightly). I found this post to be a helpful guide on how to use it.
Thanks.

You are missing one important piece of information in your question and that is the prototype definition for the InPlaceFFT() function. The most important thing to be aware of is that the FFT returns complex valued numbers. It can be that the function expects the data to be formatted such that real and imaginary values are interleaved or it could be that the function expects the real parts to be placed first in the buffer and the imaginary part to be placed in the end of the buffer. It could also be that the function returns data bit reversed (see section "Data reordering, bit reversal, and in-place algorithms" here http://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm).

Another potential issue is the length of the buffer you pass as argument, your temp buffer. You might have to pass a buffer with a length that is at least 1024*sizeof(float) or 2048*sizeof(float), if you are sure that the InPlaceFFT() function expects a float* as input. If the InPlaceFFT() function is constrained to work on real input, a buffer size of 1024 might be sufficient. But you should take a better look at the documentation for how the InPlaceFFT() works in order to find out how large the size of you temp buffer should be and how the output is laid out in the temp buffer after it returns. Another thing that looks strange is that you don't need to pass the FFT size to InPlaceFFT().

If you are sure that the sample rate is always 44.1kHz it is fine to hardcode 44.1kHz into the code, otherwise you might want to make your code more flexible by loading the sample rate into a variable and use that in your code.

A 1024 point FFT returns 1024 values (complex). You only need the first 513 output values (these values are typically called bins), the rest are in principle redundant (although the can be handy to have). If you then want to compute the magnitude spectrum you need to compute the absolute value of the complex numbers (so your abs() function should work on complex numbers) and you are in principle done. The only thing left to do is to correct the values from index 1 to index 511 by multiplying them with 2. The values at index 0 and index 512 are real and need no correction. If you are going to plot the magnitude spectrum in a GUI you probably want to plot the magnitudes in dB. So you need to compute the log (base 10) of the magnitude values and multiply by 20.

It would be useful if you could update your post with the documentation for InPlaceFFT().


Difference between Fourier transform of $\left(\frac{1}{3}\right)^{-k}$ with $k\le -1$ and $\left(\frac{1}{3}\right)^{k}$ with $k\ge 1$ When I want to calculate the discrete Fourier transform of $\frac{1}{3}^{-k}$ with $k\le-1$ I get that:
$$
\sum_{k=-{\infty}}^{-1}\left(\frac{1}{3}\right)^{-k}e^{-j\omega k} = \sum_{k'=1}^{\infty} \left(\frac{1}{3}e^{j\omega}\right)^{k'} = \frac{\frac{1}{3}e^{j\omega}}{1-\frac{1}{3}e^{j\omega}}
$$
But if I change the variable before finding the Fourier transform I have the function $\left(\frac{1}{3}\right)^k$ with $k\ge1$ (both functions are the same). The discrete Fourier transform of the new function is as follows:
$$
\sum_{k=1}^{\infty}\left(\frac{1}{3}\right)^{k}e^{-j\omega k} = \sum_{k=1}^{\infty}\left(\frac{1}{3}e^{-j\omega}\right)^{k} = \frac{\frac{1}{3}e^{-j\omega}}{1-\frac{1}{3}e^{-j\omega}}
$$
Are these two the same?

As Atul has pointed out, your both sequences are not the same, but they are time-reversed versions of each other:

k = np.arange(-15, 15)

s1 = (1/3.)**(-k)*(k <= -1);
s2 = (1/3.)**(+k)*(k >= +1)

plt.stem(k, s1, 'r')
plt.stem(k, s2, 'g')

Accordingly, what you see in your calculations is the time-reversal property of the Discrete-Time Fourier Transform:

$$x(t)=y(-t) \Leftrightarrow X(f)=Y(-f)$$

I.e. the frequency changes its sign in your calculations, as is expected from the time-reversal property.


Feature extraction/reduction using DWT For a given time series which is n timestamps in length, we can take Discrete Wavelet Transform (using 'Haar' wavelets), then we get (for an example, in Python) -

>>> import pywt
>>> ts = [2, 56, 3, 22, 3, 4, 56, 7, 8, 9, 44, 23, 1, 4, 6, 2]
>>> (ca, cd) = pywt.dwt(ts,'haar')
>>> ca
array([ 41.01219331,  17.67766953,   4.94974747,  44.54772721,
        12.02081528,  47.37615434,   3.53553391,   5.65685425])
>>> cd
array([-38.18376618, -13.43502884,  -0.70710678,  34.64823228,
        -0.70710678,  14.8492424 ,  -2.12132034,   2.82842712])

where ca and cd are approximation and detailed coefficients. Now if I use all of them I can construct my original time series back using inverse DWT. But instead I want to use a fewer coefficients (like in Fourier Transform if we use only first few coefficients, we can approximately reconstruct the original time series). If I just use ca or just use cd I don't get the desired results. If I use only we coefficients from each of them (like first 4), I get only half of the time series.

How should I select the coefficients (from ca and cd) so that I can approximately create the original signal from them (i.e. most of its energy)?

I think it is kind'a similar to soft and hard thresholding using in wavelet de-noising. Have you come across this topic? pywt has already an in-built function for this purpose. Please take a closer look at this code and try to play with it:

import pywt
import matplotlib.pyplot as plt
import numpy as np

ts = [2, 56, 3, 22, 3, 4, 56, 7, 8, 9, 44, 23, 1, 4, 6, 2]

(ca, cd) = pywt.dwt(ts,'haar')

cat = pywt.thresholding.soft(ca, np.std(ca)/2)
cdt = pywt.thresholding.soft(cd, np.std(cd)/2)

ts_rec = pywt.idwt(cat, cdt, 'haar')

plt.close('all')

plt.subplot(211)
# Original coefficients
plt.plot(ca, '--*b')
plt.plot(cd, '--*r')
# Thresholded coefficients
plt.plot(cat, '--*c')
plt.plot(cdt, '--*m')
plt.legend(['ca','cd','ca_thresh', 'cd_thresh'], loc=0)
plt.grid('on')

plt.subplot(212)
plt.plot(ts)
plt.hold('on')
plt.plot(ts_rec, 'r')
plt.legend(['original signal', 'reconstructed signal'])
plt.grid('on')
plt.show()

This will produce following - hope this will guide you.

Note on running the code with the latest versions:

if you get an error like this:

module 'pywt' has no attribute 'thresholding'

use on line 9 and 10:

 cat = pywt.threshold(ca, np.std(ca)/2, mode='soft')
 cdt = pywt.threshold(cd, np.std(cd)/2, mode='soft')

if you get an error on:

plt.hold('on') then comment this line out:

  # plt.hold('on')


Why do we really split pulse shaping between the transmitter and the receiver? EDIT: There was a small error in my code when plotting the Squared RC. Below is the correct frequency response, and clearly RC² doesn't fulfill Nyquist's criterion.

---------------------------------- Original question (to discard)

So first let me start by asking if using a Raised Cosine filter at both the transmitter and receiver (matched filter) gives a response that is ISI free? (Assuming an ideal channel)

In other words is the SQUARED (not root-squared) frequency response of RC $H_{RC}(f)^2$ still fulfills Nyquist's ISI criterion?

I avoided a direct derivation since the mathematical expressions are quite intimidating so I turned to Matlab and quickly plotted the frequency response of both the RC & RC squared and tried to calculate: $\sum H(f-k/Ts)$ and check if it's constant.

The graphs show that both $H_{RC}(f)$ and $H_{RC}(f)^2$ add up to a constant after performing the previous operation. Moreover I performed an ifft on $H_{RC}(f)^2$ and plotted it and it seems to have zeroes at $kTs$ (for $k\neq 0$) just like the RC time domain function.

So bets are, if my Matlab code is correct, that RC² fulfills Nyquist's criterion. Am I missing something here?

Now if the RC² is also ISI free, why wouldn't we use an RC at both the transmitter and receiver? Why do we turn to a root-raised cosine instead?

Why did we go from an RC to a Root-RC? Most textbooks say it's to split the filtering between TX & RX and to get an overall ISI free response but the same could be done with a RC filter.

Thanks in advance.

There must be something wrong in your code because it's quite straightforward to show that the squared frequency response of a raised cosine pulse does not satisfy the Nyquist criterion.

A raised cosine pulse satisfies the Nyquist criterion simply because a raised (co)sine function and an inverted raised (co)sine function add up to a constant:

$$(1+\sin(x))+(1-\sin(x))=2$$

where one of the two expressions corresponds to the rising edge of the frequency response whereas the other corresponds to the falling edge of a shifted version of the frequency response.

Now if you square the frequency response you get the following expression in the region where two shifted spectra overlap:

$$(1+\sin(x))^2+(1-\sin(x))^2=2(1+\sin^2(x))\neq const$$

which is not constant, and which proves that the squared frequency response of a raised-cosine pulse does not correspond to a Nyquist pulse.


Equalizing $1/f$ noise I am recently acquiring data through ADC on signals ranging from $0.1\textrm{ Hz}$ to $50\textrm{ kHz}$ ($F_s = 128\textrm{ kHz}$). I noticed that apart from the white noise, I have something similar to $1/f$ noise with a cutoff frequency near $1\textrm{ Hz}$. There isn't much that I found on equalizing this $1/f$ noise on google. I am wondering what type of filters or algorithms that I can use to equalize the noise and have a flat noise bed in my frequency range as post-processing, so that I can see my input signals without the decay of magnitude caused by $1/f$ noise.

Also, I am wondering if I am able to equalize this noise, will it cause distortion in my signal in the time domain, or phase response?

Any advice, links to references or help would be greatly appreciated.

I guess that your problem may not really be 1/f noise, but rather the effects from using the FFT directly on your time signal (without any windowing). If that is true, then try applying the Hann window to the time signal before the FFT and check the result.

There is plenty of information available about the window functions and how to use them.


Why is pitch detection necessary in pitch correction? I am trying to understand the process of changing the pitch of a signal. After performing the FFT on a windowed signal, we have the frequency information for the signal. Why can we not simply change these frequencies to change the pitch of the signal? Instead, we must first use pitch detection methods to determine pitches that were present in the time-domain function. How does the FFT not give you this information? What is the frequency information given by the FFT if not the frequencies of the time-domain signal?

I don't know what algorithm you are talking about exactly, as you didn't state that in your question, but my idea of pitch correction consists of detecting the actual pitch and shifting it to the right one.

You ask why pitch detection is necessary. Well, to be able to write a program that automatically corrects the pitch of a recording (i.e. in opposition to a scenario where you are manually checking what frequency is the right one for each instant of time) you should first have an array containing all possible frequencies corresponding to musical notes. That is, a list containing all the values that would be correct if detected, so that they wouldn't have to be changed. If the singer has accurately hit a $\mathrm{C_3}$, then why would you auto-tune it?

Then, you should perform pitch detection using STFT. For each instant of time, there will be a dominant frequency. One (primitive) way to do this would be to find the frequency that shows the maximum amplitude. When you find that maximum (pitch detection), you should compare the frequency where it's located to the array of musical notes you created before. Then, assuming the singer is not that bad and the tone he's hitting is approximate to the desired one, one should take the detected frequency to the closest musical note (pitch shifting). For example, if the maximum is detected at $\mathrm{124 \ Hz}$, one would go for a $\mathrm{B_2}$ ($\mathrm{123.47 \ Hz}$) instead of a $\mathrm{C_3}$ ($\mathrm{130.81 \ Hz}$). I believe there are some things about logarithmic perception that should be taken into account here, but that's not relevant in this question. Take into account that there are pitch detection algorithms, I just used this easy method as an example to make it illustrative.

That's why detection has to be made. You have no way to automatically "round" the pitch to the right one if you don't know what actual pitch the singer is in.


Normalized cross correlation via FHT - how can I get correlation score? I'm using the 2D Fast Hartley Transform to do fast correlation of two images in the frequency domain, which is the equivalent of NCC (normalized cross correlation) in the spatial domain.

However, with NCC, I can get a confidence metric that gives me an idea of how strong the correlation is at a certain offset. In the frequency domain version, I end up with a peak-finding problem in the inverse FHT after doing the correlation, so my question is:

Can I use the value of the peak that I find in the correlation image to derive the same (or similar) confidence metric that I can get from NCC? If so, how do I calculate it?

Ok, I finally found the reason for my confusion - thanks to this answer.

There were a number of reasons why I couldn't get my spatial NCC and frequency domain peak values to agree numerically. Some obvious, some subtle, all stupid. For the sake of completeness, here they are:

I had missed an offset and scale factor that was introduced when
normalising my image data after importing the pixel values in
preparation for FHT.
I was using a predicted maximum value for the
    peak, from my sub-pixel shift calculation instead of the actual
    peak.
I had a low-pass filter in the pipeline, which - for
    performance - I was applying concurrently with the correlation
    operation. So my peak value was scaled by the filter.

Not being mathematically inclined, I was wrongly attempting to "explain" my problem by assuming that using FHT was somehow inherently different to using FFT. After being pointed in the right direction I found the above issues and finally got the results I was expecting.


wavelet denoising routine for environmental data series I have a time series of water temperature, e.g.

y = 1+(30-10).*rand(1,365);

I have previously used the wavelet denoising routine in the wavelet toolbox by matlab to remove unwanted noise from a signal, e.g.

[s_denoised, ~, ~] = wden(y, 'minimaxi', 's', 'sln', 2, 'db4');

However, I am unsure whether I am using this method correctly! If my time series is of water temperature which is characteristic of an environmental data series i.e. high power at low frequency (e.g. seasonal cycle) and low power at high frequency (e.g. diel cycle) and I wish to remove noise caused by internal wave motion in lakes (which cause sudden increases or decreases in temperature), should I be using the wavelet shown here or is there a better wavelet to use?

First, a comment - before you denoise, you are basically going to be converting your data from the (time)-domain into the wavelet domain. This is nothing but a series of projections of your data unto user-picked basis functions. (The wavelets).

When you denoise, you will be zeroing out, or shrinking, coefficients which (ostensibly) belong to the noise.

Now, imagine that you picked a wavelet type, and when you transformed your noiseless signal, the very first coefficient was some number, and the rest were zero. This means that you know, that your pure signal, will always give you a non-zero as the first coefficient in the wavelet domain, with other coefficients always being zero.

Now imagine that you transformed your noisy signal. You would probably still see a high value for that first coefficient, (belonging to your true data), but all the other coefficients that you were expecting to be zero have some small value. If you zero them out and inverse transform, you have denoised your signal.

The question you want to answer then becomes, "Which wavelet type, when projected against my data, gives me the sparsest representation in the wavelet domain?" In other words, which wavelet, when used against your (pure, noiseless) data, gives you coefficients that are mostly zero?

Mathematically, what you want is a wavelet that has as many vanishing moments as your data template.

Example: For example, let us say that know that your data template (in your case, the temperature without internal wave motion corruption) is a polynomial of 0th order, (a constant), then you want to use the Haar wavelet, since it has 1 vanishing moment. This means that if you transform your pure-signal into the Haar domain, most of the coefficients will be zero. However when you transform your data-with-internal-wave-corruption, the coefficients that were normally zero now have some value, that you can then clobber and inverse transform.

The same thing can be said for other templates. Let us say that your pure-data was a polynomial with maximum degree of 7. That means then that you want a wavelet that makes my polynomials of maximum degree 7 disappear, or in other words, has 8 vanishing moments. This would make your signal look quite sparse in the transform domain, and thus any other noise of higher order would make it less sparse, and you would be able to denoise by removing those suspicious coefficients that were supposed to be zero.


Least Squares with Blocks / Updates I have a continuous-time system that I want to fit via least squares. I just send $N$ digital samples $x[n]$ through the system and receive (via analog signal chain, ADC etc) $N$ digital samples $y[n]$. My system can be described as

$$
\mathbf{y} = \mathbf{X} \mathbf{c}
$$

where $\mathbf{X}$ is a matrix composed of $x[n]$. Hence I solve for the coefficients via:

$$
\hat{\bf c} = (\mathbf{X}^{H} \mathbf{X})^{-1} \mathbf{X}^{H} \mathbf{y}
$$

which works wonderful and it's straight forward and easy. Note - $(\cdot)^H$ denotes the Hermitian transpose

Now $\mathbf{y}$ are heavily contaminated by noise so I need many samples. But on each receive call I can only receive a certain maximum number of samples. Even if I would receive an infinite amount - calculating the pseudo inverse becomes very expensive.

What is the easiest way and direct generalization of the above least squares set up so that I can iteratively receive blocks and refine the estimate after each block? The method should match a hypothetical, full least squares solution as closely as possible.

There are a couple of concepts which may be related but I am not sure how they relate to each other:

I could first receive, say 1000 batches of repetitive $\mathbf{y}$ and average them and solve one least squares. Is this equivalent to solving a larger least squares system with all samples? Disadvantage: I need to send repetitive data.
I am not looking for LMS or similar sampled-based methods. Ideally I would just use Least Squares ... if the number of measurements would not be as excessive
What about RLS oder IRLS? To my understanding, they are also sample based (not block based) methods.
I have heard about gradient descent + Newton type:

$$
\mathbf{c}_{n+1} = \mathbf{c}_{n} = \mu \left(\frac{\partial^2 J}{\partial \mathbf{c}^{*} \partial \mathbf{c}}\right)^{-1} \frac{\partial J}{\partial \mathbf{c}^{*}}
$$

which, in my linear model would expand to:

$$
\mathbf{c}_{n+1} = \mathbf{c}_{n} = \mu (\mathbf{X}^{H}\mathbf{X})^{-1} \mathbf{X} \mathbf{e}
$$

where $\mathbf{e}$ is the error vector between actual measurements and model with the current estimate $\mathbf{c}_n$. It requires that the unknowns are smaller than the number of samples in each iteration. This looks very close to what I am looking at ... but I am not sure if it is the right thing and how it relates to the other methods.

RLS is derived using the Woodbury Identity for rank-1 update.

I would go back and re-derive the RLS approach in light of Woodbury and generalize it to the block update case. I vaguely recall doing this before, but I can't recall when/what-for at the moment.

Here's the identity:

$$
(A+UCV)^{-1} = A^{-1} - A^{-1}U(C^{-1}+VA^{-1}U)^{-1}VA^{-1}
$$

[Aside:  I will try to use upper case for matrices, lower case with underbar for column vectors, and simply lower case for scalars - though I don't think we will have any scalars. I also use superscript $H$ to denote Hermitian transpose]

Now, consider your data vector and matrix at time 'k' and the resulting LS solution:

$$
\underline{c}_k = arg\min_{c_k}||\underline{y}_k - X_k\underline{c}_k||_2^2
$$

Then, at time 'k+L' - after another block of 'L' samples, we have

$$
\underline{c}_{k+L} = arg\min_{\underline{c}_{k+L}}||\underline{y}_{k+L} - X_{k+L}\underline{c}_{k+L}||_2^2
$$

Next, define some terms for the block update

$$
\underline{y}_{k+L} = \left[\begin{array}{l}\underline{y}_k\\ \underline{\tilde{y}}_{k+L}\end{array}\right]
$$
$$
{X}_{k+L} = \left[\begin{array}{l}X_k\\ {\tilde{X}}_{k+L}\end{array}\right]
$$

Now, solve for $\underline{c}_{k+L}$ using the Woodbury identity (in line 3), state, and the solution from time $k$:

\begin{eqnarray}
\underline{c}_{k+L} &=& (X^H_{k+L}X_{k+L})^{-1} X^H_{k+L}\underline{y}_{k+L}\\
&=&(X^H_{k}X_{k}+\tilde{X}^H_{k+L}\tilde{X}_{k+L})^{-1} (X^H_{k}\underline{y}_{k} + \tilde{X}^H_{k+L}\underline{\tilde{y}}_{k+L})\\
&=&[(X^H_{k}X_{k})^{-1} + (X^H_{k}X_{k})^{-1}\tilde{X}^H_{k+L}(I+\tilde{X}_{k+L}\tilde{X}^H_{k+L})^{-1}\tilde{X}_{k+L}(X^H_{k}X_{k})^{-1}][X^H_k\underline{y}_k + \tilde{X}^H_{k+L}\underline{\tilde{y}}_{k+L}]\\
&=&\underline{c}_k+W^H_{k+L}\underline{\tilde{y}}_{k+L} + [W^H_{k+L}(I+\tilde{X}_{k+L}\tilde{X}^H_{k+L})^{-1}W_{k+L}](\underline{v}_k+\tilde{X}^H_{k+L}\underline{\tilde{y}}_{k+L})
\end{eqnarray}

where

$$
W_{k+L} = \tilde{X}_{k+L}(X^H_{k}X_{k})^{-1}
$$
$$
\underline{v}_k = X^H_k\underline{y}_k
$$
$$
\underline{v}_{k+L} = \underline{v}_k + \tilde{X}^H_{k+L}\underline{\tilde{y}}_{k+L}
$$

I would keep state for the latter intermediates and for $(X^H_{k}X_{k})^{-1}$. The expensive part here is usually the calculation of $(I+\tilde{X}_{k+L}\tilde{X}^H_{k+L})^{-1}$, since it is a rank $L>1$ update.

Apologies if I have typos above. I was deriving as I typed, and eventually I couldn't see the wysiwyg output anymore. ;) Hopefully, you can see the jist of the approach, nonetheless. There certainly could be further intermediate and state optimization. I didn't go back to RLS to remind me of how it was wrung to minimize operations, but it probably goes somewhere along these lines.

Finally, given the above block update approach, the question for you is - which is preferred:  $L$ rank 1 updates (RLS) or 1 rank $L$ ("block RLS" as above) update? This depends on your situation and environment. A block approach may be better if the rank-1 updates would be done in an interpreted loop (e.g. Matlab), and one may be better than the other for accumulated precision errors.


2D Frequency Domain Convolution Using FFT (Convolution Theorem) In the time domain I have an image matrix ($256x256$) and a gaussian blur kernel ($5x5$). I've used FFT within Matlab to convert both the image and kernel to the frequency domain as zero padded $260x260$ matrices      ($N + M -1 = 256 + 5 -1 = 260$)
I then multiply the image matrix by the kernel and use IFFT to convert the result back to the time domain. When I try to display the result, it is just junk and doesn't resemble the original image with a gaussian blur like it should.
Here is the Matlab code I am using, where image = $256x256$ and kernel = $5x5$:
imagefreqdomain = fft2(image,260,260)

kernfreqdomain = fft2(kernel,260,260)

filtimagefreqdomain = imagefreqdomain * kernfreqdomain

filtimage = ifft2(filtimagefreqdomain)

What am I doing wrong? Thanks

Similar to your question Applying 2D Image Convolution in Frequency Domain with Replicate Border Conditions in MATLAB the issue is what happens when you multiply in 2D in frequency domain.
So few remarks about that:

Multiplying in frequency domain for discrete signals with finite support is equivalent to applying convolution in spatial domain under the assumption of cyclic / periodic boundary conditions.
In image processing we usually define per kernel the anchor pixel of the kernel. Usually it is marked as (0, 0) of the kernel. We also mostly set it as the center pixel (In Image Processing most kernels have odd length). When we pad the kernel to the size of the image we usually add zeros on its bottom and right. Which means its (0, 0) isn't aligned with the image.

The misalignment with the circular boundary extension yields the following for the naïve code:
clear();
close('all');

gaussianKernelStd       = 0.5;
gaussianKernelRadius    = ceil(5 * gaussianKernelStd);

mI = im2double(imread('cameraman.tif'));
mI = mI(:, :, 1);

numRows = size(mI, 1);
numCols = size(mI, 2);

vX = [-gaussianKernelRadius:gaussianKernelRadius].';
vK = exp(-(vX .* vX) ./ (2 * gaussianKernelStd * gaussianKernelStd));
mK = vK * vK.';
mK = mK ./ sum(mK(:)); %<! The Gaussian Kernel

mIFiltered = ifft2(fft2(mI) .* fft2(mK, numRows, numCols), 'symmetric');
figure();
imshow([mI, mIFiltered]);

As you can see at the top and left filtered image (The right) has artifacts which are the result of the circular extension and the misalignment. How to fix it?
Well, padding the image correctly and padding with circular extension the kernel.
I showed it in Applying 2D Image Convolution in Frequency Domain with Replicate Border Conditions in MATLAB.


Is There a Sparse Representation for Noise? Is there sparse representation for stationary noise and nonstationary noise?
How can I learn dictionary for each noise class?
(my mean of noise is noises with which speech signals are often contaminated such as white gaussian noise, car noise, babble noise, impulsive noise....)

Let's think about it in a different way - Generate Noise from a Dictionary.

Let's create a Dictionary $ A \in \mathbb{R}^{m \times n} $ where each of its rows is normalized (Has Euclidean Norm of $ 1 $) and generated by a Gaussian Random Vector.

Now, let's create $ N $ random vector $ {\left\{ {r}_{i} \right\}}_{i = 1}^{N} $ by:

$$ {r}_{i} = A {g}_{i} $$

Where $ {g}_{i} $ are random vectors $ {g}_{i} \in \mathbb{R}^{n} $ with only $ k \ll n $ non zero elements (Let's say they are pre defined and deterministic otherwise we'll have to deal with multiplication of random variables)..

Now, the output vectors are valid Gaussian Random Variables and given enough of them and some properties of $ A $ you'd be able to build a Dictionary which they are sparse with relation to it.

Yet it means the dictionary is random by itself and next time when you get new set of random data it won't be able to represent it correctly (It doesn't generalize).

While the above is not rigorous (It means to develop some intuition, not analytic analysis) it tells you why when we deal with finite amount of data it can be done but for noise as in general it won't work.


Problem in verifying zero blocking property of a zero of a multiple input multiple output system I have a multiple input multiple output control system, and in particular a two inputs two outputs system. I want to find the zero state direction and the zero input direction of the system, which has a zero at $s=0$. So, I have foud them with the following code:

s = tf('s');

P = 1/(s+5);
C = 6/s;

S = 1/(1+P*C);
T = P*C/(1+P*C);

G_2 = [T S; S -S]; %transfer matrix

S = ss(G_2)
Smr = minreal(S)
A = Smr.A
B = Smr.B
C = Smr.C
D = Smr.D
trz = tzero(A,B,C,D)
RSM_1 = [0*eye(2)-A B; -C D]
rRSM_1 = rank(RSM_1)

null = null(RSM_1)

%verifying zero blocking property
x0 = [-0.4388; -0.1482];
u0 = [-0.0000; 0.8863];
t = linspace(0,5);
u = exp(t);
[y,x] = lsim(A,B*u0,C,D*u0,u,t,x0);
plot(t,y)
axis([0 5 -1 1])
xlabel('time, seconds')
title('response to x_{0} and u(t)=u_{0}e^{t}')

and I want to verify the zero blocking property, which means that if I apply an input at $u(t)=u_0e^{t}$ starting from $x(0)=x_0$, where $u_0$ is the input direction and $x_0$ is the state direction, the output is zero.

I have found it analytically in this way:

I know that the since the system zeros of a MIMO transfer function are those values $s = z$ which make the system matrix $S(z)$ loose rank, there exists at least a direction $u_0$ (zero input direction) and a n-dimensional vector $x_0$ (zero state direction), not simultaneously zero, such that

$\begin{pmatrix}
zI-A & -B\\
 C&D
\end{pmatrix}$ $\begin{bmatrix}
x_0\\
u_0
\end{bmatrix}$ $=$ $\begin{bmatrix}
0\\
0
\end{bmatrix}$

And by solving this system we obtain the zero state direction $x_0$ and the zero input direction $u_0$.

In order to make the question even clearer I can make an example using a different transfer function which has a simpler minimim realization. The example is taken from this : Here, and it can be found at the bottom of the paper.

The system in this case is:

$P(s)=\begin{bmatrix}
\frac{2}{s^2+3s+2} & \frac{2s}{s^2+3s+2}\\
 \frac{-2s}{s^2+3s+2}& \frac{-2}{s^2+3s+2}
\end{bmatrix}$

Its minimum realization is:

$A=\begin{bmatrix}
-1 &0  &0 \\
 0&-2 &0 \\
 0& 0 & -2
\end{bmatrix}$

$B=\begin{bmatrix}
2 &-2 \\
 -2& 4\\
 -4& 2
\end{bmatrix}$

$C=\begin{bmatrix}
1 &1  &0 \\
1 &0  &1
\end{bmatrix}$

$D=\begin{bmatrix}
0 &0 \\
 0& 0
\end{bmatrix}$

and the after using the same code, which can also be found in the paper, the simulation for verifying the zero blocking property is:

s = tf('s');
P = [2/(s^2+3*s+2) 2*s/(s^2+3*s+2); -2*s/(s^2+3*s+2) -2/(s^2+3*s+2)]
S = ss(P)
Smr = minreal(S)
A = Smr.A
B = Smr.B
C = Smr.C
D = Smr.D
trz = tzero(A,B,C,D)
RSM_1 = [eye(3)-A B; -C D]
rRSM_1 = rank(RSM_1)

null(RSM_1)

x0 = [0.5345; -0.5345; -0.5345];
u0 = [0.2673; -0.2673];
t = linspace(0,5);
u = exp(t);
[y,x] = lsim(A,B*u0,C,D*u0,u,t,x0);
plot(t,y),grid
axis([0 5 -1 1])
xlabel('time, seconds')
title('response to x_{0} and u(t)=u_{0}e^{t}')

and the output of this is:

which is correct, since the output coes to zero.

But, what I obtain with my code is:

I am still working on this, and by using the Matlab command isstable(T) it tells me that the closed loop is not stable. But if I run the step response, step(T) ,I find:

what am I doing wrong?

The transmission zeros of a MIMO system are obtained by solving

\begin{align}
\dot{x}^* = A\,x^* + B\,u^*, \\
y = C\,x^* + D\,u^*,
\end{align}

for $x^*$ and $u^*$ such that $y = 0$. By using the Laplace variable $s$ for the time derivative that system of equation can also be written as

$$
\begin{bmatrix}
0 \\ 0
\end{bmatrix} =
\begin{bmatrix}
A - s\,I & B \\ C & D
\end{bmatrix}
\begin{bmatrix}
x^* \\ u^*
\end{bmatrix},
$$

which can be seen as the generalized eigenvalue problem $\mathcal{A}\,v = \lambda\,\mathcal{B}\,v$, with $\lambda$ equivalent to $s$,

$$
\mathcal{A} =
\begin{bmatrix}
A & B \\ C & D
\end{bmatrix}, \quad
\mathcal{B} =
\begin{bmatrix}
I & 0 \\ 0 & 0
\end{bmatrix}, \quad
v =
\begin{bmatrix}
x^* \\ u^*
\end{bmatrix}.
$$

I believe that tzero in matlab also solves such generalized eigenvalue problem to obtain the transmission zeros (however by using eig yourself the associated eigenvectors are more accurate than the nullspace calculations). The associated solution in time for both $x(t)$ and $u(t)$ can be obtained using

$$
\begin{bmatrix}
x(t) \\ u(t)
\end{bmatrix} = e^{\lambda\,t}\,
\begin{bmatrix}
x^* \\ u^*
\end{bmatrix},
$$

so you can't just pick any value for $\lambda$ when generating such solution. You could now simulate the system using $x(0)=x^*$ as initial condition and the expression for $u(t) = e^{\lambda\,t}u^*$ with lsim. However, I believe that slim uses zero-order-hold for the input so the actual applied input will deviate slightly causing the output $y$ to deviate away from zero using that approach. Instead, if you look at the analytical solution, both the output and the error in state dynamics should remain practically to zero (for non-positive $\lambda$):

  n = length(A);
  m = length(D);
  [V,L] = eig([A B; C D], diag([ones(1,n) zeros(1,m)]));
  ind = 1;    % Index of associated eigenvalue (avoid infinite eigenvalues)
  lambda = L(ind,ind);
  x0 = V(1:n, ind);
  u0 = V(1+n:n+m, ind);

  N = 1e3;
  t = linspace(0, 5, N);
  x = @(t) x0 * exp(lambda* t);
  u = @(t) u0 * exp(lambda* t);
  y = @(t) C * x(t) + D * u(t);
  e = @(t) A * x(t) + B * u(t) - lambda* x(t);

  Y = zeros(n, N);
  E = zeros(m, N);
  for k = 1 : N
      Y(:,k) = y(t(k));
      E(:,k) = e(t(k));
  end
  figure, plot(t, E)
  figure, plot(t, Y)