How to implement IDFT using the positive frequencies of DFT, modulation and resampling on MATLAB? This question is a continuation of the discussion at How can a linear operator on DFT vector produce the same vector using only half of the DFT vector?, where the higher level motivation is better explained, in case there is interest.

Question:

How to take the positive frequencies from an hermitian symmetric DFT vector, take the IDFT of this (which is complex), multiply by a carrier to shift the spectrum, resample, take the real part of it, and obtain the same IDFT result as the result that would be obtained by the IDFT of the original hermitian symmetric symbol?

Details:

I would really appreciate to have my understanding and the concepts I mention checked.

Consider that the DFT vector $\mathbf{X}$ of a certain signal is hermitian symmetric and represents a bandwidth of $B$ Hz. If I consider a signal represented by the positive frequencies (from DC to Nyquist), it is the same as shifting the positive side of the spectrum, and centralizing it at frequency 0. Now, the new spectrum is not hermitian symmetric anymore, yields a complex IDFT and represents a Bandwidth of $\frac{B}{2}$.

To use a numerical example, for $N=8$, considering from DC to Nyquist, that would give 5 positive frequencies and 3 negative. Then, I would shift the 5 DFT samples (correspondent to the positive frequencies) to the center of the spectrum. In this case, the new DFT size would be odd, and the Nyquist frequency (at k= 5/2 = 2.5) wouldn't fall on a DFT sample

At this point, I have a signal bandlimited to $\frac{B}{2}$. Then, if I multiply this signal with a carrier of frequency $\frac{B}{2}$, I shift the spectrum again to the positive and negative side, such that the original spectrum is reconstructed. However, since this operation would increase the bandwidth back to $B$ Hz, the product $x[n]*c[n]$ would need to have its sampling frequency increased by a factor of 2, so that the frequency spectrum can actually represent the spectrum of the signal without aliasing.

The carrier should be an exponential:

\begin{align}
c[n] = e^{j2\pi n k_0/N}
\end{align}
where $k_0$ is the number of DFT samples to shift. In this case, $k_0 = \frac{N}{2}$, such that:
\begin{align}
c[n] = e^{j \pi n}
\end{align}

Finally, after doing resampling, I should take the real part of the signal and that should be equal to the original IDFT (the one obtained by simpling taking the IDFT of the original hermitian symmetric symbol).

However, I am confused mainly by the following issues:

From DC to Nyquist there are 5 DFT samples. That would yield a complex IFFT with size 5, Then, I would multiply point-wise by 5 samples of the carrier and upsample by 2, resulting in 10 samples. On the other hand, the original IDFT has length 8.
I'm not sure about the rationale behind taking the real part.
I'm not sure if the exponential is the correct carrier.

MATLAB implementation:

Here is the code I'm trying to make work.

N=4;
A=rand(N-1,1) + j*rand(N-1,1); % Random positive carrier
X=[0; A; 0; flipud(conj(A))]; % Hermitian symmetric DFT vector

x = ifft(X)

% Positive side of the spectrum shifted to the center
Xp = [X((N/2 + 1):N+1); X(1:N/2);];

% Complex IFFT:
xp = ifft(Xp);

% Carrier
n=(0:N);
k0 = (N+1)/2;
carrier = exp(j*2*pi*k0*n/(N+1)).';

% Modulation
x_up = xp.*carrier;

% Increase the sampling rate by a factor of two
x_re = resample(x_up,2,1);

real(x_re)

I basically want the two vectors that are printed to the workspace to be equal: x and real(x_re)

Here is the IDFT with normal and positive only frequencies if it helps:

N=4;
A=rand(N-1,1) + j*rand(N-1,1);   % Random positive carrier
X=[0; A; 0; flipud(conj(A))];    % Hermitian symmetric DFT vector
Xp=[0; 2*A; 0; zeros(size(A))];  % Positive freq only DFT vector

x = ifft(X);         % Normal IDFT
xa = ifft(Xp);       % Analytic signal
xp = real(ifft(Xp)); % Same as normal IDFT


Trignometric Fourier series representation of a continous time signal While learning Fourier series I read the definitions of representation for a continuous time signal $x(t)$ as:

$$x(t)=A_0 + 2 \sum_{k=1}^{\infty} A_k \cos(k \omega_0 t) - B_k \sin(k \omega_0 t) \tag{1}$$

where $A_k$ and $B_k$ are real.

Another definition was:
$$x(t)=a_0 + \sum_{k=1}^{\infty} a_k \cos(k \omega_0 t) + b_k \sin(k \omega_0 t) \tag{2}$$

Why two different definitions?  I mean in $(1)$ there is a coefficient of $2$ for sinusoids and sine coefficients are negated.

alright, let's review a little bit of Euler before we get to the Fourier.

$$ e^{j \theta} \ = \ \cos(\theta) \ + \ j \sin(\theta) $$

from that you can get

$$ \cos(\theta) = \frac{e^{j \theta} + e^{-j \theta}}{2} \quad\quad\quad \sin(\theta) = \frac{e^{j \theta} - e^{-j \theta}}{2j} $$

so now let's look at Eq (1)

$$ \begin{align}
x(t) \ &= \ A_0 \ + \ 2 \sum_{k=1}^{\infty} A_k \cos(k \omega_0 t) - B_k \sin(k \omega_0 t) \\
 &= \ A_0 \ + \ 2 \sum_{k=1}^{\infty} A_k \frac{e^{j k \omega_0 t} + e^{-j k \omega_0 t}}{2} \ - \ B_k \frac{e^{j k \omega_0 t} - e^{-j k \omega_0 t}}{2j} \\
 &= \ A_0 \ + \ \sum_{k=1}^{\infty} (A_k + jB_k) e^{j k \omega_0 t} \ + \ (A_k - jB_k) e^{-j k \omega_0 t} \\
 &= \ \sum_{k=-\infty}^{\infty} c_k \ e^{j k \omega_0 t} \\
\end{align} $$

where

$$ c_k =
\begin{cases}
A_{-k} - jB_{-k}, & \quad \text{for } k < 0 \\
A_0, & \quad\quad k = 0 \\
A_k + jB_k, & \quad\quad k > 0
\end{cases} $$

and, going in the other direction,

$$ \begin{array}{lcl}
A_0 & = & c_0 \\
A_k & = & \Re\{c_k \} \quad\quad \text{for } k>0\\
B_k & = & \Im\{c_k \} \quad\quad\quad k>0
\end{array} $$

note that for real $A_k$ and real $B_k$, then

$$ c_{-k} = c_k^* = \text{"complex conjugate of } c_k \text{ "}$$

so, for Eq (1), the simplicity of having the real and imaginary parts of the $c_k$ coefficients be simply $A_k$ and $B_k$ (for $k \ge 0$) is the motivation for the convention of the leading "2" before the summation and for the minus sign.  it's just a convention.  this convention is useful because it's much easier to derive the coefficients $c_k$ than it is to derive the coefficients $A_k$ and $B_k$.

$$ c_k \ = \ \frac{\omega_0}{2 \pi}\int_{t_0 - \pi/\omega_0}^{t_0 + \pi/\omega_0} x(t) \ e^{-j k \omega_0 t} \ dt $$

where $ -\infty < t_0 < +\infty $ can be any convenient real value.

the convention for Eq. (2) comes about from looking at Fourier series first as a real analysis problem without the use of complex variables.  then this is a simple first statement:

$$x(t) \ = \ a_0 \ + \ \sum_{k=1}^{\infty} a_k \cos(k \omega_0 t) + b_k \sin(k \omega_0 t)$$

as you can see $a_k$ and $b_k$ are related to $A_k$ and $B_k$ in a very simple and straight forward manner.  you can derive that simple relationship.  but the formulae for getting $a_k$ and $b_k$ from $x(t)$ (and $\omega_0$) is less straight forward to derive and to express.  if fact, it's two equations, not one simpler equation.


When convolving two functions that are constants in a region and 0 everywhere else, where does the integration start? Heads up, this is for homework.  I never took a signals and systems course, so I'm behind on this stuff.

I want to compute the convolution of two rectangular regions.  I know the standard equation for convolution, where the integrals are from $-\infty$ to $\infty$, but it is useful to restrict the integration to the areas that matter.

I have one function that has a value of $1/3$ when $-0.5 \le x \le 0.5$ and $0 \le y \le 0.4$.

I have a second function that has a value of $1$ when $-1/50 \le x,y \le 1/50$.  Where would the bounds for integration be for convolving functions like this?

For simplicity, let's talk about one-dimensional convolutions only.

You need to understand that the result of a convolution is a function and the alleged convolution integral is actually an infinitude of integrals, one for each time instant.  By this I mean that to compute the output $y(t)$ of a system with impulse response
$h(t)$ to an input signal at time $t = 5$, we write
$$\begin{align*}
y(5) &= \int_{-\infty}^\infty h(\tau)x(5-\tau)\mathrm d\tau,\\
\end{align*}$$
and compute $y(5)$.  For $t = 6.5$, we write
$$\begin{align*}
y(6.5) &= \int_{-\infty}^\infty h(\tau)x(6.5-\tau)\mathrm d\tau,\\
\end{align*}$$
and compute $y(6)$.  For $t = 10.2$, we write
$$\begin{align*}
y(10.2) &= \int_{-\infty}^\infty h(\tau)x(10.2-\tau)\mathrm d\tau,\\
\end{align*}$$
and compute $y(10.2)$.and so on.  After doing this a while, some smart-ass
said, "Hey Ma, I am beginning to see a pattern here!  Lookit, I can
write
$$\begin{align*}
y(t) &= \int_{-\infty}^\infty h(\tau)x(t-\tau)\mathrm d\tau, ~\text{for all} ~t\\
\end{align*}$$
and the world was never the same again.

Suppose that $h(t)$ is $0$ for $t < 0$ or $t > 12$.  Then all the integrals
exhibited above simplify in the sense that the limits can be changed to
$0$ and $12$ since the integrand, being a multiple of $h(\tau)$, is $0$
when $\tau < 0$ or $\tau > 12$.  Suppose in addition that $x(t)$ is nonzero
only when $|t| < 3$.  Now look at the first integral (for $y(5)$ with the
limits changed to $0$ and $12$).  At $\tau = 0$, the integrand is
$h(0)x(5-0) = 0$ since $x(5) = 0$.  As $\tau$ increases towards $12$, the
argument of $x$ in the integrand decreases from $5$ downwards till the
argument hits $-7$ at $\tau = 12$.  Thus, the integrand is nonzero
only for $\tau \in (2,8)$, having value $h(2)x(5-2) = h(2)x(3)$
at one end of the interval and value $h(8)x(5-8) = h(8)x(-3)$ at
the other.  We conclude that
$$y(5) = \int_{-\infty}^\infty h(\tau)x(5-\tau)\mathrm d\tau
= \int_{0}^{12} h(\tau)x(5-\tau)\mathrm d\tau
= \int_{2}^{8} h(\tau)x(5-\tau)\mathrm d\tau.$$
Similarly,
$$y(6.5) = \int_{-\infty}^\infty h(\tau)x(6.5-\tau)\mathrm d\tau
= \int_{0}^{12} h(\tau)x(6.5-\tau)\mathrm d\tau
= \int_{3.5}^{9.5} h(\tau)x(6.5-\tau)\mathrm d\tau.$$
However,
$$y(10.2) = \int_{-\infty}^\infty h(\tau)x(10.2-\tau)\mathrm d\tau
= \int_{0}^{12} h(\tau)x(10.2-\tau)\mathrm d\tau
= \int_{7.2}^{12} h(\tau)x(10.2-\tau)\mathrm d\tau$$
(can you see why?)

So, there is not a simple answer to your question about the
limits of integration; you have to work them out for each
value of $t$ for which you want to evaluate the convolution
integral. But you don't have to do an infinite amount of
work!  There are Look Ma moments here and there.  Can you
combine two of the results above to write that for each $t$,
$3 \leq t \leq 9$, we have the useful pattern
$$y(t) = \int_{t-3}^{t+3} h(\tau)x(t-\tau)\mathrm d\tau$$
shown above?

I hope this will help you gain the understanding that
you need in order to complete your homework.


Upscaling detection approach I'm trying to implement a solution that would, to some extent, evaluate a picture telling what are the chances that it was upscaled from the smaller image.

My approach is as follows:

Take a picture of unknown origin (scaled or not)
Perform sharpening using some technique which sharpens only evident edges (1)
Extract edges using Difference of Gaussians
Measure entropy of resulting image

Then downscale the image and upscale it back to the same size and compare entropies.
OR
Then upscale the image and downscale it back to the same size and compare entropies.

Does it sound reasonable? What technique should I use as (1)? Also, maybe I should do it in small steps so I don't introduce too many artificial points? I mean if the image was upscaled (preferably instantly, to target size) it won't matter a lot - but it the image weren't upscaled this incremental improvement will be better for the picture than taking it and upscaling instantly?

Here is the poor man's approach:

var histogram = []
foreach pixel
    histogram.append(abs(pixel.value - pixel.south.value))
    histogram.append(abs(pixel.value - pixel.east.value))
sort(histogram)
sharpness = histogram[histogram.count * 0.95]
upscaled = sharpness < X // You choose X experimentally


historical context of impulse response-convolution Hi All: I have zero background on the history of LTI systems so I was wondering when the concept of LTI and the fact that an LTI signal can be be completely described by its impulse response came about. I ask because, in econometrics, there are models called auto-regressive distributed lags (ARDLs). A specific ARDL called the koyck distributed lag came out in the early 1950's and was a pretty big deal in the econometrics community.

https://www.reed.edu/economics/parker/312/tschapters/S13_Ch_3.pdf
( koyck is on page 12 ).

But, now that I understand LTI, impulse response and convolution, my guess is that the explanation and theory behind the Koyck ADl could have been much more easily explained and understood using the LTI framework ( by using the impulse response and convolution). Yet, as far as I am aware, an explanation of the Koyck distributed lag using this approach does not exist. That's very odd to me since the LTI framework makes the understanding of the Koyck ADL so much clearer and easier to understand. Thanks for any historical context of the LTI concept.  If  it's after the 1950's, then it's understandable that Koyck didn't know about it yet. If it was well known before the 1950's, then my guess is that he was unaware of it. Otherwise, he would have explained his Koyck model using that approach.

Stanley: Thanks for your input. Your answer led me to the "brilliant" idea of googling myself. Based on the link below, I think one can define the birth of convolution as sometime between 1760-1900 depending on how one defines birth. https://pulse.embs.org/january-2015/history-convolution-operation/.

Either way, it seems that Koyck was unaware of it since  he came  up with his model in the 1950's. He would have had a much easier time using the DSP framework in his derivation but he came up with the model without impulse response or convolution which is quite impressive in itself. I'm gonna check this answer but your mention of Mr. Heavside is appreciated.


Conceptual Questions on Colored Noise Process I am having a tough time finding answers to some specific questions and finding references where there is information regarding Brownian noise or Red Noise. I'm referring to white and colored noises where white noise can be Gaussian or Uniform Noise where the power spectrum is flat. Along the same line we have colored noise whose frequency is $1/f^{\alpha}$ where $\alpha =1$ for pink noise and $\alpha =2$ for brown or red noise. Matlab Code to generate colored noise presents the code to generate colored noise by filtering white noise simulated using randn.
This question is inspired from a previous question asked here: Stationary vs non-stationary signals?
Shall be obliged for the answers to the  following specific questions. Please correct me wherever I am wrong.

Is brown noise stationary or non-stationary? If stationary then, the signal does not change with time. Then how come it appears to be random?

If a signal is corrupted with brown noise and pink noise respectively, how would the noisy signal's characteristics differ? I think since brown and pink noise samples are correlated, the resulting noisy signal's samples if it was uncorrelated when clean, would now become correlated? Not sure about this. It's important for me to know so that I can apply specific noise filtering techniques for correlated and uncorrelated signals.

An example of non-stationary process?

Do the terms "noise" and "process" mean the same?

Once you talk about the Spectrum of noise / process you implicitly says it is stationary in the wide sense.
What does it mean to have a signal with uncorrelated samples? Do you understand it means you can't using linear predictor, to have any information from all past samples ion the current one? Usually this is not how signals behave. This is exactly why correlated noise is hard. It behaves like a signal hence it is harder to separate them apart. One way to deal with the model of Signal + Colored Noise is to apply whitening transform on the data, apply White Noise methods and then apply the inverse transform.
Any process where the probability for a specific value depends on the index value.
Here are some specific examples (From Investopedia - An Introduction to Stationary and Non-Stationary Processes):

Noise is used when dealing with signals. We call noise to the non deterministic added (Usually, we have multiplicative noise as well) term of the signal. Since we have some index (Time / Spatial, etc...) for the signal we have it for the noise as well which means it is a process.


FFT Frequency bin relationship I guess this is a stupid question or it is not clear.  I will try again

  I have a situation where I take the fft of a combined source, both at the same frequency but at different locations. I am trying to figure out what the relationship between the source and the combined.

  Maybe I am using the wrong terminology here.  A combined source in my application is at least two instances of the same frequency from different locations in the focal plane, combined in the same waveform.

  Is the phase of the combined complex value the sum of the phase of the source? so one could write something like:

  tanCombined = tan(combinedIm / combinedRe) = tan(s1Im/s1Re +
  s2Im/s2Re + ...);

Second Attempt.

Say I have a waveform waveOne that is a sin wave at a specific frequency that has an origin at a point in the focal plane relative to a single measuring microphone.

If I take the fft of this waveform I can calculate the phase angle and amplitude from the fft data, phaseOne and amplitudeOne.

Say I have another waveform waveTwo that is a sin wave of the same frequency that has its origin at a different point in the focal plane.

If I take the fft of this wave form I can calculate the phase angle and amplitude from the fft, phaseTwo and amplitudeTwo.

If I combine these waveforms into a single wave, waveThree and take its fft the same bin is excited as in the above two samples and both the amplitude and phase angle are different, as one would expect.

The question is, for the combined experiment, is there a relationship between the amplitude and phase of the combined fft value with those of the individual waveform.

Yes, simply use a trigonometric addition formula.  The results can vary from complete constructive interference to most destructive interference.  Depending on the relative phases at your microphone.

Here is the math:

$$ W_1 = A_1 \cos( F t + P_1 ) $$
$$ W_2 = A_2 \cos( F t + P_2 ) $$

$$ W_1 = A_1 \cos( F t ) \cos( P_1 ) - A_1 \sin( F t ) \sin( P_1 ) $$
$$ W_2 = A_2 \cos( F t ) \cos( P_2 ) - A_2 \sin( F t ) \sin( P_2 ) $$

$$ W_3 = W_1 + W_2 = [ A_1 \cos( P_1 ) + A_2 \cos( P_2 ) ] \cos( F t ) - [ A_1 \sin( P_1 ) + A_2 \sin( P_2 ) ] \sin( F t ) $$

Now

$$ W_3 = A_3 \cos( F t ) \cos( P_3 ) - A_3 \sin( F t ) \sin( P_3 ) $$

so

$$ A_3 \cos( P_3 ) = A_1 \cos( P_1 ) + A_2 \cos( P_2 ) = C $$
$$ A_3 \sin( P_3 ) = A_1 \sin( P_1 ) + A_2 \sin( P_2 ) = S $$

Calculate $C$ and $S$.

$$ A_3 = \sqrt{ C^2 + S^2 } $$
$$ P_3 = \operatorname{atan2}( S, C ) $$

If they are in phase $(P_1 = P_2)$ then $ A_3 = A_1 + A_2 $.

If they are completely out of phase $( P_3 = P_1 = P_2 + \pi )$ then $ A_3 = A_1 - A_2 $.

Also, the FFT is a linear operator, so the FFT of your sum is the sum of your FFTs.

Hope this helps.

Ced


Why zero padding the 2-d DFT interpolates images in spatial domain? I was applying different image interpolation techniques and I came know to about interpolation in frequency domain. In this technique we first take 2d DFT of an image, pad it with zeros and take the inverse DFT. The magnitude of inverse DFT gives interpolated image. It is working quite well!

Now my question is that how can it interpolate image in spatial domain when we are just padding zeros in frequency domain? How it constructs data points? And why only use zero padding?

Now I know about Shannon interpolation theorem (sinc interpolation) in which we can ideally reconstruct continuous signal from its discrete counterpart (if the sampling frequency is greater than Nyquist frequency). We can construct unknown values by sinc interpolation, but I can not understand how inverse DFT is constructing values for interpolated signal?

For the sake of simplicity, I'll explain the 1-D case; the 2-D case is completely analogous. Let $x[n]$ be a finite length sequence with $n\in[0,N-1]$. Its discrete Fourier transform (DFT) is

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N} kn}\tag{1}$$

The sequence $x[n]$ can be obtained from $X[k]$ via the inverse DFT (IDFT):

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N}kn}\tag{2}$$

Note that the sequences defined by $(1)$ and $(2)$ are periodic with period $N$. Let's use $\tilde{X}[k]$ to denote the periodic sequence of DFT coefficients. For odd $N$, $\tilde{X}[k]$ just equals (the periodic continuation of) $X[k]$. For even $N$, we split the coefficient $X[N/2]$ in two equal terms $\frac12 X[N/2]$, such that for even $N$, the periodic sequence $\tilde{X}[k]$ is defined as

$$\tilde{X}[k]=\\\big[\ldots ,\underbrace{X[0],X[1],\ldots,\frac12 X[N/2],\frac12 X[N/2],\ldots,X[N-1]}_{N+1\textrm{ elements}},\ldots\big]\tag{3}$$

With that definition, the period of $\tilde{X}[k]$ is always odd. Let's denote that odd period by $2K+1$. For odd $N$ we have $K=(N-1)/2$, and for even $N$ we get $K=N/2$.

Using $\tilde{X}[k]$ and $K$ as defined above, we can rewrite $(2)$ as

$$x[n]=\frac{1}{N}\sum_{k=-K}^{K}\tilde{X}[k]e^{j\frac{2\pi}{N}kn}\tag{4}$$

With these preparations we can now interpret the coefficients $\tilde{X}[k]$ as Fourier coefficients of a continuous periodic function $x_c(t)$ with period $NT_s$, where $T_s$ is an arbitrary sampling interval:

$$x_c(t)=\frac{1}{N}\sum_{k=-K}^{K}\tilde{X}[k]e^{j\frac{2\pi}{NT_s}kt}\tag{5}$$

The sequence $x[n]$ can be obtained by sampling $x_c(t)$:

$$x[n]=x_c(nT_s)\tag{6}$$

But we can also obtain other, more densely sampled sequences by sampling $x_c(t)$. We choose the new sampling interval such that the new sequence is also periodic with period $M>N$:

$$\hat{x}[m]=x_c\left(m\frac{NT_s}{M}\right)=\frac{1}{N}\sum_{k=-K}^{K}\tilde{X}[k]e^{j\frac{2\pi}{M}km}\tag{7}$$

The new length $M$ sequence $\hat{x}[m]$ can be interpreted as an interpolated version of the original length $N$ sequence $x[n]$.

Note that $(7)$ can be written as a length $M$ IDFT with a zero-padded sequence $\tilde{X}_0[k]$:

$$\begin{align}\hat{x}[m]&=\frac{M}{N}\left(\frac{1}{M}\sum_{k=0}^{M-1}\tilde{X}_0[k]e^{j\frac{2\pi}{M}km}\right)\\&=\frac{M}{N}\cdot \textrm{IDFT}_M\big\{\tilde{X}_0[k]\big\}[m]\tag{8}\end{align}$$

with

$$\tilde{X}_0[k]=\big[\tilde{X}[0],\ldots,\tilde{X}[K],\underbrace{0,\dots,0}_{M-2K-1\textrm{ zeros}},\tilde{X}[K+1],\ldots,\tilde{X}[2K]\big]\tag{9}$$

In terms of the original DFT coefficients $X[k]$, the zero-padded sequence $\tilde{X}_0[k]$ is given by

$$\tilde{X}_0[k]=\big[X[0],\ldots,X[(N-1)/2],\underbrace{0,\dots,0}_{M-N\textrm{ zeros}},\\\qquad X[(N+1)/2],\ldots,X[N-1]\big]\tag{10}$$

for odd $N$, and by

$$\tilde{X}_0[k]=\big[X[0],\ldots,\frac12 X[N/2],\underbrace{0,\dots,0}_{M-N-1\textrm{ zeros}},\\\qquad \frac12 X[N/2],\ldots,X[N-1]\big]\tag{11}$$

for even $N$.


Optimum matched filter without ISI Given a filter used to shape the digital signal, $p(x)$, and given that we do not want the filter combination to cause any ISI, what "matched" filter, $q(x)$ will maximise the SNR?

Matched filters are used in digital communications to maximize the signal to noise ratio.  Often a root-raised-cosine filter is used to shape the signal, since it is bounded in frequency space and the same filter can be applied to the received signal to improve the signal-to-noise ratio (SNR) without causing inter-symbol-interference (ISI).

However if a less optimum filter is used to shape the signal, then using the same filter at the receiver can introduce ISI.  It is not immediately obvious what the best choice of filter at the receiving end is.

My understanding is that the SNR is maximized by maximizing $\int{p(x)q(x)dx}$, so I want to maximize this while satisfying the constraint that filters cause no ISI ($p(x)*q(x) = 0$ for $x=kT$, $k$ is an integer, $T$ is symbol width).

Presumably one could do this by solving a Euler-Lagrange equation with some lagrange multipliers for the constraints.  Is there an easier way, or am I making a mistake, or going in the wrong direction?

For the case of linear modulation on the AWGN channel with equiprobable symbols (a very common case) the optimum approach is to truly use a filter that is matched to the symbol waveform, i.e.:

$$
q(x) = p(x)
$$

Using a matched filter provides the optimum signal-to-noise ratio at the filter output at each decision instant. This is easy to see when you remember that a matched filter acts like a sliding cross-correlator between its input signal and the expected symbol waveform, correlating the two at all possible lags. At optimum decision instants, the filter's impulse response (typically scaled to have unit energy) lines up exactly with a transmitted symbol, analogous to a zero-lag condition on the cross-correlation operation. At this time value, the filter's output is equal to the amount of energy in the received symbol, scaled by a data-dependent factor (e.g. for BPSK, the matched filter would output $ E_s $ or $ -E_s $), plus a noise term.

The noise energy at the filter output during the sampling instant is not dependent upon the time-domain shape of the filter's impulse response, only the impulse response's total energy (as noted previously, typically unity). Therefore, the signal to noise ratio is maximized by maximizing the amount of signal energy in the filter output at the sampling instant. By choosing the receiver filter to be matched to the symbol shape, we have done so, as the symbol waveform has maximum correlation with a filter impulse response that has an identical shape. Thus, the matched filter provides maximum SNR, for the AWGN channel case.

With that bout of hand-waving out of the way (you can definitely get at it with more mathematical rigor, but I'm an engineer and this is a free service; if you want to dig into the details, check any digital communication theory text), you might be thinking that I forgot that you asked about the non-ideal, ISI case. Fear not, for I assert that if you know the transmitted pulse shape, the matched filter is still the optimum choice for the AWGN channel.

The key: if you know the responses of the pulse-shaping and receiver detection filters $p(x)$ and $q(x)$, and the last "few" transmitted symbols, you can calculate what the ISI induced by those previous symbols would be and account for it accordingly; it is a deterministic quantity. The amount of symbol history that you require is related to the amount of ISI you have, i.e. how many symbol periods the cascaded filter response smears across.

Of course, you typically don't know with certainty what the previous few symbols were; if you did, then you might be at a high enough SNR that your ISI can be neglected. In the more interesting case, you can't make that assumption. Instead, a maximum-likelihood sequence detection approach is employed using the Viterbi algorithm. This process is referred to as Viterbi equalization, because in this model you treat the ISI induced by the pulse shape like a soft-valued convolutional code that is applied to your transmit waveform. The time duration of the ISI in the Viterbi equalizer defines the required number of algorithm states, similar to the constraint length in a convolutional code.

This approach is often used in systems that have the non-optimum pulse shape that you noted; one notable example is GSM (which uses a Gaussian pulse shape that extends across multiple symbol intervals). One great reference on this topic was published by Sklar in 2003:

B. Sklar, “How I learned to love the trellis“, IEEE Signal Processing Magazine, pp. 87-102, May, 2003


Calculating an in-loop signal as part of a hierarchical control loop I've got a control system with two feedback paths, with each path going to a different actuator that corrects the error in the system. One feedback path provides feedback at low frequencies, and the other at high frequencies. Here is the loop:

I want to calculate the signal on one of the actuators (the node at the output of the "slow" block) while the loop is closed, given a disturbance that I inject before the summing junction. The signal I want to calculate is "A". The input is normally zero (I just want to suppress disturbances to zero, not some set-point).

I initially started by calculating the open loop gain provided by the inner plant-slow-sum loop, i.e. A = slow * plant * disturbance, and then divided the disturbance by this open loop gain to give me something that approximates the disturbance suppressed by the closed loop, but this is not correct because in reality some of the disturbance is also suppressed by the outer plant-fast-sum-sum loop. In effect, the disturbance that the "slow" loop will see at high frequencies will be much lower, since the "fast" loop suppresses it there.

How do I calculate the closed-loop signal at A due to the disturbance, including the effect of the outer loop?

UPDATE: My original answer incorrectly assumed the outer loop would not impact the result, but as pointed by Ben, this transfer function is dependent on the fast feedback.

The resulting equations are found by solving for $A/X$ from the block diagram below, with X representing the disturbance input:

Solving this:

$$Y = P(SY+FY+X)$$
$$A/S = P(A + FA/S + X)$$
$$A = PAS + PFA + PXS$$
$$A(1- PS - PF) = PXS$$

$$\frac{A}{X} = \frac{PS}{1-PS-PF}$$

Note that this follows the general equation for the closed loop gain using:

$$G_{CL} = \frac{G_F}{1+ G_{OL}}$$

Where:

$G_{CL}$: Closed Loop Gain

$G_{F}$: Forward Gain (In this case $G_{F}  = PS$)

$G_{OL}$: Open Loop Gain (In this case $G_{OL} = -P(S+F)$

Note how negative feedback is assumed in the generic expression, but since the block diagram does not show negation at the summations, the open loop gain is then negative as given above. Using the general equation and recognizing that the parallel result of S and F as a gain given as S+F allows us to solve this from inspection.


Is there a way in MATLAB to approximate an FIR filter using an IIR Filter? I am a student trying to understand how this is possible. I've done some research on the web with not many sources showing up with solutions to this problem that are usable. What techniques are used to give an approximation of an FIR filter using IIR filters with say the same order?

Sure there is.

The simplest version would be to calculate the transfer function of your FIR filter, sample it on some meaningful frequency grid and then shove it into invfreqz() to design the IIR filter.

Depening on you FIR filter this may work well or be a total disaster. In my experience invfreqz() does often very poorly.

If that's not good enough you can certainly write your own search an optimization algorithm. However, this is NOT a beginner project and requires significant math and coding skills.

See mathworks.com/help/signal/ref/invfreqz.html#mw_3c79f471-99ad-46f8-a680-14c29ff47d07


Understanding unexpected i/q calculation on a downsampled signal I am reading samples with an RTL-SDR, centered on an FM broadcast station:

with closing(RtlSdr()) as sdr:
    sdr.sample_rate = 2.4e6
    sdr.center_freq = 88.9e6
    sdr.gain = 'auto'
    x1 = sdr.read_samples(16384)

The read_samples() function normalizes the data:

>>> sdr.read_samples.__doc__
Read specified number of complex samples from tuner. Real and imaginary
parts are normalized to be in the range [-1, 1]. Data is safe after
this call (will not get overwritten by another one).

This returns an array that I can work with, with numpy, scipy, etc.

I create arrays for the real and imaginary parts:

x3_real = x3.real; x3_imag = x3.imag

For each pair in (for i,q in zip(x_real, x_imag)), I calculate h=math.sqrt(i**2 + q**2).  I don't expect h to ever be larger than sqrt(2), because the real and imaginary parts are normalized to [-1, 1].

Then I downsample to several different frequencies, less than the sampling frequency, using scipy.signal's decimate() function:

x3 = signal.decimate(x3, ds_factor, zero_phase=True)

and repeat the hypotenuse calculations.  In x2 I do not see numbers > sqrt(2), however I see numbers > sqrt(2) in the decimated signal (x3).  Why is this?

Because the decimate function does not try to keep power constant.

Remember that decimation is not inherently an LTI operation; so, probably, the authors of that function simply did not care too much for keeping the amplitude of a passband signal constant.

Note that a decimation step typically includes filtering to $\frac1{\text{ds_factor}}$ of the original bandwidth to avoid aliases. That filter might or might not have a passband gain $\le 1$.

A general remark:

  For each pair in (for i,q in zip(x_real, x_imag)), I calculate h=math.sqrt(i**2 + q**2). I don't expect h to ever be larger than sqrt(2), because the real and imaginary parts are normalized to [-1, 1].

In fact, if coming right out of a receiver with good IQ balance, don't expect $h\ge1$, even!

Remember that I and Q aren't independent; they are the inphase and quadrature component of a signal whose max power, no matter which phase, is limited! Thus, as long as you don't risk clipping in your receiver, $i^2 \le {1-q^2}\iff q^2\le 1-i^2$!

If your I and Q would be independently limited to an amplitude of 1, your $\sqrt 2$ would be right. They are not – imagine, for example, the passband signal that gives you $i\equiv1$: it's the cosine with the carrier frequency $1\cdot\cos(2\pi f_c t + 0)$ (note the $0$ phase). Now, any cosine with larger power would clip, right? Obviously, cosines of the same amplitude would also be at the extreme power limit, no matter whether they have zero phase. The moment the phase is not zero, you'll see some of the cosine's energy appear on Q, and missing on I. However, it's the same power.


What is the difference of each frequency response of partial filters in a polyphase method I try to understand how works polyphase method in term of frequency response of partial filters.

Of course, the common is that the filters operate at low frequency, but what differency between them. How it works?

The short answer is the polyphase filter converts a low pass filter into a series of all pass filters each with a different time delay. So it is a series of delays at even fractions of the time between samples of the lower sampling rate of the polyphase filter. By getting an output of the same signal at different fractional delays, we can combine these to have a higher sampled version of the same signal. Polyphase means many phases, and the ideal implementation would be all-pass filters modifying the delay only- in good polyphase filter design we approach this ideal.

The longer answer is the polyphase filter is simply an array of filters formed by decimating the coefficients of a FIR filter, resulting in each filter output being identical to each other in their decimated (reduced) frequency response and each with a fractional delay of the other given the decimation relationship:

$$z^{-m}D\downarrow  = z^{-m/D}$$

Meaning a delay of $m$ samples at the higher rate will be a delay of $m/D$ samples at the lower rate after decimating by D.

This turns a low pass filter that would have a passband covering the frequency range after decimation into all pass filters, each with a different delay, hence the name poly-phase meaning many phases (it is actually many delays).

Consider the following 96 tap FIR filter with coefficients generated from MATLAB/Octave:

coeff = firls(96, [0 .2 .3 1], [1 1 0 0])

Which has this frequency response:

The time delay in samples is given by the slope of the phase vs frequency consistent with using $z = e^{j\omega}$ for $z^{-1}$: one cycle delay will have a phase that extends negatively to $\pi$ as the normalized angular frequency extends to $\pi$. Linear phase FIR filters (such as this one) have a delay that is half the number of taps, so in this case the filter delay is 48 samples, so the linear phase slope would be given by $z^{-48} = e^{-j48\omega}$ which is consistent with the phase slope we see in this filter's passband.

If we decimate this 96 tap filter into four separate 24 tap filters, the sampling rate for each filter will be at 1/4 of the rate above. Decimation is resampling and causes all the higher frequency spectrum to fold into the new sampling rate, while maintaining the original spectrum that was at the lower rate (if that is not clear, see this post The condition in order to not loose information after upsampling and downsampling). So if we decimated the impulse response of the filter (the coefficients), the frequency response below $0.25\pi$ rad/sample will remain the same other than an amplitude scaling by 4 or -12 dB, while the response above that will fold in, distorting the primary response if we didn't reject the higher bands sufficiently. For example, here is the frequency response if we take the first sample and every 4th sample after that in our original 24 tap filter, to create a 6 tap filter:

freqz(coeff(1:4:end));

So now the filter is nearly an "all-pass" filter modifying only the phase with a delay. It is passing all signals within it's new sampling frequency that is 1/4 of the original rate, and it's delay is 12 samples instead of 96: The phase at $\omega = \pi$ would be given by the phaes of $e^{-12j\omega}$; $12\pi  180/\pi = 2160°$, consistent with the maximum phase in the plot below.

If we select the next sample, and every 4 samples after that, we will end up with the same frequency response, since we are decimating the same filter, except the phase will be advanced by one sample at the higher sampling rate, or as we first explained $1/4$ sample at the lower output rate. Thus we would expect the phase slope to decrease to $11.75$ samples (with the associated smaller phase slope).

Superimposing all four filters given by coeff(1:4:end), coeff(2:4:end), coeff(3:4:end), and coeff(4:4:end) reveals that at least out to $0.8\pi$ radians/sample, these are now all pass structures with incrementally varying delay. The magnitude is completely flat out to $\omega = 0.8\pi$ which would be the usable range of our filter, we assume there is no signal of interest above $\omega = 0.8\pi$:

And the group delay (which is the negative derivative of phase with respect to frequency, and the time delay for linear phase filers) result is 12, 11.75, 11.5 and 11.25 samples:

The four delays are very conveniently at 1/4 samples of the higher 4x sampling rate! So each filter passes our signal within the passband with no distortion but provides 4 different outputs each at a time offset of the underlying waveform of 1/4 of a sample.

Given a signal exemplified by the sinusoid in the plot below, each color group would represent the output of one of the four polyphase filters.


Is this system invertible or not? Prove that the following system is invertible.
$$y(t) = \mathcal{T}\{x(t)\} = \int_{-\infty}^{3t} x(\tau) \,\mathrm d \tau$$

Answer: yes, the system is invertible. I need some hint here, not the full solution.

In this problem of finding the inverse system (if it exists) its intuitive to try differentiating the integral as the system input/output is given by:

$$y(t) = \mathcal{T}\{x(t)\} = \int_{-\infty}^{3t} x(\tau) d\tau$$

Before differentiating te integral however, I would like to make this little change which is quite clear I assume:
$$y(t/3) = \mathcal{T}\{x(t)\} = \int_{-\infty}^{t} x(\tau) d\tau$$

Then I differentiate both sides of the equality:
$$ \frac{d}{dt} \left( y(t/3) \right) = \frac{d}{dt} \int_{-\infty}^{t} x(\tau) d\tau $$

Which proceeds as:
$$ \frac{1}{3} y'(t/3) = x(t) $$ which yields the inverse system.

Note that differentiation of an integral with variable limits is known as the Leibnitz Rule and can be summarized as below:

given
$$F(x) = \int_{\alpha(x)}^{\beta(x)} g(x,t) dt$$
then
$$ \frac{d}{dx} F(x) = \frac{d}{dx} \int_{\alpha(x)}^{\beta(x)} g(x,t) dt $$
$$ F'(x) = g(x,\beta(x)) \beta'(x) - g(x,\alpha(x)) \alpha'(x) +  \int_{\alpha(x)}^{\beta(x)} \frac{d}{dx} g(x,t) dt $$


Why do SDR receivers emit I/Q samples? To my understanding hardware receivers for software-defined radio applications basically take the input signal, mix it with the tuning frequency to remove the carrier frequency and then sample the resulting voltage with a sampling rate that is just high enough for the bandwidth of the payload signal. They emit those samples to the demodulation software in the form of I/Q value pairs. I'm assuming they obtain the Q value by taking another sample $1/4$ cycle (with respect to the tuning frequency) later, effectively doubling the sample rate.

Why do they use the I/Q representation?

I can see how I/Q is a nice representation (in hardware) when synthesizing signals because you can for example do frequency or phase modulation just by varying amplitudes, but this reason doesn't seem to apply to the case of SDR receivers.

So, is there anything gained by using I/Q for output instead of I at twice the sampling rate? Or is it just a matter of convention?

The SDR (or any general digital signal processing system) takes the received RF signal, and downconverts it from the carrier frequency to the baseband.

Now, the real bandpass signal from the antenna does not necessarily have a symmetric spectrum around the carrier frequency, but it can be arbitrary. If the downconverter now shifts the spectrum to the center frequency, the corresponding time-domain signal becomes complex. So, the I and Q samples you get from the SDR are the real and imaginary part of the complex baseband signal, that corresponds to your real passband signal around the carrier frequency.

More details can e.g. be found at the wikipedia site for digital down conversion.

To answer your question:

I/Q representation does not correspond to different sampling points of the signal. Instead, it corresponds to the real and imaginary part of the digital complex-valued baseband signal. These parts are obtained by separately multiplying the RF signal with a sine and a cosine and sampling both streams after low-pass filtering.

Sampling with double frequency can yield the same information as I/Q. It would be neccessary to homodyne the signal to $f_s/4$ to get all information that would have been in the baseband IQ signal to be in the passband signal at $f_s/4$ (where $f_s$ is the sampling rate).


Inverse Z-transform mystic simplification I have the following expression:

$$X(z) = \frac{16}{15}\frac{1}{1-\frac14z^{-1}} - \frac{16}{15}\frac{1}{1-4z^{-1}}$$

According to my understanding this should become:

$$x(n) = \frac{16}{15}\left(\frac14\right)^n  u(n) - \frac{16}{15} 4^n  u(n)$$

But according to my source it is:

$$x(n) = \frac{16}{15}\left(\frac14\right)^n  u(n) + \frac{16}{15} 4^n  u(-n-1)$$

Are these expressions equal?

If they are, how can $-u(n) = u(-n-1)$ ?

If not, why?

I used $a^nu(n)\Longleftrightarrow  \frac{1}{1-az^{-1}}$.

This is probably homework, so I'll just give you a few hints to get you on the right track. First of all, given an expression $X(z)$, the corresponding sequence is usually not unique. Your solution is correct if the region of convergence (ROC) of $X(z)$ is $|z|>4$, which gives you an exponentially increasing sequence. If this were the impulse response of a discrete-time system, then the system would be causal but unstable.

The other solution assumes that the ROC is $\frac14<|z|<4$, which includes the unit circle. This corresponds to a two-sided decaying sequence. As an impulse response it would correspond to a non-causal but stable system.

So, to answer your question, both solutions are correct but they are not identical, because they are inverse transforms of $X(z)$ assuming different regions of convergence.


convolution vs correlation? What is difference between convolution and correlation?In simple words or in a nut shell?

As far as i am able to study is that Convolution is a mathematical operation used to express the relation between input and output of an LTI system. It relates input, output and impulse response of an LTI system while Correlation is a measure of similarity between two signals

Convolution:

$$ y(t) = h(t) \circledast x(t) = \int\limits_{-\infty}^{\infty} h(u) \, x(t-u) \ \mathrm{d}u $$

Cross Correlation:

$$ R_{xy}(t) = \int\limits_{-\infty}^{\infty} y(u) \, x(t+u) \ \mathrm{d}u $$

The difference between the two is effectively the sign on $u$ in $x(t-u)$ in the integral.  That correlation is like convolution but with one of the signals flipped left to right is essentially the basis of the Matched Filter.


How to Enhance this finger print? I am trying to extract finger print from this image by using green LED backlight but when it convert to binarize image most of finger print is gone.

gray= green

I use this fomular to convert from this color image

to this gray image

I use otsu for binarization. So Which preprocessing I should use before convert to binarization to emphasis the finger print.

Two steps,

do a simple histogram equalization to make the brightness a little
bit more even.
Then use canny edge detector (as suggested by Marcus Muller).

matlab code:

I=imread('Your Image');
G1=rgb2gray(I);
G2=histeq(G1);
E2=edge(G2,'canny');
figure imshow(E2,[]);

Result (I hope it is what your are looking for):


Why are there two phases for one dominant frequency? I was playing with a set of time history, trying to extract its phases corresponding to the each distinguishable dominant frequency.

Just look at the figures: I got one of the dominant frequency in Fig.1, but in Fig. 2 there are two different phases close to the dominant frequency (f=0.05). Here I am confused; I do not understand why there are two phases near the dominant frequency.

Edit: it seems when I increase the threshold of masking noises, the negative phase disappeared.

Fig. 1

Fig. 2

Appendix: the code I used for phase extraction

    %     figure;plot(f,P1);
%     xlabel 'Frequency (Hz)';xlim([0 2.5])
%     ylabel 'A' %ylabel 'degrees' %
%     grid

    Yb=fft(y);
    ly = length(Yb);
    fb = (0:ly-1)/ly*Fs;

%     X3=Yb;
%     threshold = max(abs(Yb))/2; %tolerance threshold
%     %X3(abs(Yb)<threshold) = 0; %maskout values that are below the threshold
%     phs = unwrap(angle(X3));
%     %phs = (angle(X3));
%     figure;plot(fb,phs*180/pi)
%     %figure;plot(fb,(mod(phs*180/pi,360)))
%     xlabel 'Frequency (Hz)';xlim([0 2.5])
%     ylabel 'Phase (degrees)' %ylabel 'degrees' %
%     grid

    X1=Yb;%store the FFT results in another array
    %detect noise (very small numbers (eps)) and ignore them
    [maxYb, maxin]=max((abs(Yb))/ly);
    threshold = max((abs(Yb)))/2; %tolerance threshold
    X1(abs(Yb)<threshold) = 0; %maskout values that are below the threshold
    phase=atan2(imag(X1),real(X1))*180/pi; %phase information

    %phase(sum(fb<=0.05))% need to adjust manually
    %fb(sum(fb<=0.05))   %need to adjust manually

%     figure;plot(fb,phase); %phase vs frequencies
%     %disp(['at ',num2str(),])
%     xlabel 'Frequency (Hz)';xlim([0 2.5])
%     ylabel 'Phase (degrees)' %ylabel 'degrees' %

    disp(['Amax=',num2str(maxYb)])
    fmax=fb(maxin);disp(['f for Amax=',num2str(fmax)])
    phsmax=phase(maxin);disp(['phs for Amax=',num2str(phsmax)])
    phs=phsmax;

    %inverse fft test
%     X1i=ifft(X1);
%     figure;     plot(t,X1i)
%     hold on;    plot(t, y)

If the dominant frequency is not exactly integer periodic in the FFT width, then most of that frequency peak's energy will be split between the two nearest FFT result bins (two different frequencies).  The phase of your frequency peak might be slightly ahead or behind the phase of the slightly lower FFT frequency bin. If so, it will have the opposite phase relationship (slightly behind or ahead) of the the FFT frequency bin slightly on the higher or opposite side in frequency difference.

So it's not two different phases.  It's two opposing phase relationships when your dominant frequency spectra is correlated against the two nearest frequency FFT basis vectors, that are both slightly different (and one integer period of the FFT length apart), but on opposite sides in frequency, and thus opposite in phase difference relationship to your dominant frequency peak.

Unless, of course, your dominant peak represents a signal that is exactly integer periodic in the FFT aperture, in which case only one FFT result bin (per harmonic) will be involved, and you won't see two phases (for any non-zero FFT result bins that are above the noise floor)


Autoconvolution vs Autocorrelation I have a serious doubt about Autocorrelation and AutoConvolution.
My understanding is,

Autocorrelation

Autocorrelation is nothing but correlation of the signal and it's delayed version. Example : Autocorrelation of 2 functions f(t) & f(t+1).

Autoconvolution

Autoconvolution is combination of  the two signals where the output signal will be a modified version of a signal in term of another signal with every point/unit.

My Question

Are two signals are the same if their auto-correlation functions are the same?
Are two signals are the same if their auto-convolution functions are the same?

Are two signals are the same if their auto-convolution functions are the same?

Almost. Look at the autoconvolution in the frequency domain where the autoconvolution of $x$ (with itself) gives us $(X(f))^2$ in the frequency domain while the autoconvolution of $-x$ (with itself) gives us $(-X(f))^2 = (X(f))^2$. So, given an autoconvolution function, there are two (very related) signals $x$ and $-x$ that have the same autoconvolution function.

  Are two signals are the same if their auto-correlation functions are the same?

Not quite.  Now we are given $|X(f)|^2$ in the frequency domain and there are many different factorizations possible. For example, if $y(t)$ is a signal such that values taken in by its Fourier transform always lie on the unit circle in the complex plane (for every $f$, $|Y(f)|=1$) then $|X(f)Y(f)|^2 =  |X(f)|^2$ and so $x\star y$ has the same auto-correlation function as $x$. Note that $x(t)$ and $x(t-\tau)$ (which is a delayed version of $x(t)$) have the same autocorrelation function ($Y(f)$ happens to be $\exp(-j2\pi f \tau)$ here). Another factorization replaces $X(f)$ by $X^*(f)$ which tells us that $x(t)$ and $x(-t)$ (which is just $x(t)$ running backwards in time) have the same autocorrelation function.


Why power and not amplitude? I'm reading 'The Scientist and Engineer's Guide to Digital Signal Processing'. It says for a noisy signal, the important parameter is not the deviation from the mean, but the power represented by the deviation from the mean.

So, if I were to consider Amplitude it would just be voltage deviations from the mean, but for power it becomes Voltage^2 deviations. Then an average is taken and a square root to give the standard deviation [It would just be an average deviation in the case of the Amplitude].

I'm trying to understand why considering just the Amplitude is not good enough?
Why is it that when I'm dealing with noisy signals, I need to consider power calculations to try and reduce the noise?

It's generally not true that the average deviation from the mean is not important or much less important than the power of the deviation. For sampled data, the mean deviation would be defined as

$$\epsilon=\frac{1}{N}\sum_{n=0}^{N-1}|x_n-\mu_x|\tag{1}$$

where $\mu_x$ is the mean of the data. This quantity is definitely relevant. However, the absolute value in $(1)$ makes it much less convenient to use than the equally relevant standard deviation

$$\sigma=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}|x_n-\mu_x|^2}\tag{2}$$

which is differentiable, and, consequently, easier to handle analytically (such as in least squares fitting, etc.). (Note that in Eq. $(2)$ you might as well use a factor $1/(N-1)$ instead of $1/N$, but that is not relevant here).

Apart from the above argument based on mathematical convenience, it is true that in the case of uncorrelated noise sources, it is the noise variances (i.e., the squared standard deviations) that add up, and not the mean deviations.


Signals - question due to sampling I don't have an idea how to start it. Help pls :)

$$m = z(t) + \frac{\operatorname{Sa}(t-T)}{4}$$

We have some signal $z(t)=\left(\operatorname{Sa}(\frac{t}{T})\right)^2$. How often the m signal should be sampled to be sure of correct signal reconstruction

This is not an answer.  I am just trying to restate the question so that the equation is dimensionally homogeneous:

I don't have an idea how to start it. Help pls :)

$$m(t) = z(t) + \tfrac{1}{4}\operatorname{sinc}\left(\tfrac{t-T}{T}\right)$$

We have some signal $z(t)=\left(\operatorname{sinc}(\tfrac{t}{T})\right)^2$. How often should the $m(t)$ signal be sampled to insure correct signal reconstruction?

this question can be meaningfully answered.

let's define:

$$ \operatorname{sinc}(u) \triangleq \begin{cases}
\frac{\sin(\pi u)}{\pi u} \qquad & u \ne 0 \\
1 \qquad & u = 0 \\
\end{cases} $$

$$ \Pi(u) \triangleq \begin{cases}
1 \qquad & |u| < \tfrac12 \\
\tfrac12 \qquad & |u| = \tfrac12 \\
0 \qquad & |u| > \tfrac12 \\
\end{cases} $$

$$ \Lambda(u) \triangleq \begin{cases}
1 - |u| \qquad & |u| \le 1 \\
0 \qquad & |u| > 1 \\
\end{cases} $$


Reduce Jitter in Live Kalman Filter I have a stream of data over time that I have fed into a Kalman filter. This data represents the vertical displacement in time of a vehicle's rearwheel as it travels over a bump in the road. The wheels are being detected through video processing.

How can I fine tune my algorithm to reduce the sawtooth jitter that the Kalman filter gives me?

code:

import numpy as np
import pandas as pd
from pykalman import KalmanFilter
from numpy import ma
import matplotlib.pyplot as plt

# The numpy.ma module provides a convenient way to address the issue of dropout, with masked arrays.
# When an element of the mask is False, the corresponding element of the associated array is valid and is said to be unmasked.
# When an element of the mask is True, the corresponding element of the associated array is said to be masked (invalid).
rearwheel_y = ma.masked_values(rearwheel_y, 0)

# time step
dt = 1
# initial_state_mean
X0 = rearwheel_y[0]
# initial_state_covariance
P0 = 0
n_timesteps = len(rearwheel_y)
n_dim_state = 1
filtered_state_means = np.zeros((n_timesteps, n_dim_state))
filtered_state_covariances = np.zeros((n_timesteps, n_dim_state, n_dim_state))
# Kalman-Filter initialization
kf = KalmanFilter()
# iterative estimation for each new measurement
for t in range(n_timesteps):
    if t == 0:
        filtered_state_means[t] = X0
        filtered_state_covariances[t] = P0
    elif t != 0:
        filtered_state_means[t], filtered_state_covariances[t] = (
        kf.filter_update(
            filtered_state_means[t-1],
            filtered_state_covariances[t-1],
            observation = rearwheel_y[t])
        )

# plot of the resulting trajectory
plt.figure()
plt.plot(rearwheel_y, 'k.', label = 'Observations')
plt.plot(filtered_state_means[:, 0], "g-", label="Filtered Positions", markersize=1)
plt.ylim(150, 205)
plt.grid()
handles, labels = plt.gca().get_legend_handles_labels()
by_label = dict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys())

Let's start the right way (as Tim's comment says): let's start with the signal model of what you're doing.
Let's assume that the measurements at time $k$ are just positions $x_k$ and $y_k$ and that these are also part of the state, which also includes $\dot{x}_k$, $\dot{y}_k$, $\ddot{x}_k$, and $\ddot{y}_k$:
$$
\mathbf{x}_k = \left[   x_k \ \dot{x}_k \  \ddot{x}_k \ y_k \  \dot{y}_k \
 \ddot{y}_k   \right]^T
$$
Then the output equation is
$$
\mathbf{z}_k = \mathbf{H} \mathbf{x}_k + \epsilon_k = \left [ x_k\ y_k \right ]^T + \epsilon_k
$$
with
$$
\mathbf{H}  = \left [  \begin{array}{cccccc} 1 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 0 & 0
\end{array} \right ]
$$
and $\epsilon_k$ is the measurement noise and it is distributed jointly normally with zero mean and covariance $\mathbf{R}$ (called the observation_covariance in pykalman). The observation_matrices is what $\mathbf{H}$ is called in pykalman.
The state update equation will then just be
$$
\mathbf{x}_{k+1} = \mathbf{A} \mathbf{x}_k + \mathbf{B} \eta_k
$$
where $\mathbf{A}$ is the state transition matrix and $\eta_k$ is the process noise.
For a standard displacement, velocity, acceleration model driven by noisy jerk inputs, this means:
$$
\begin{align}
x_{k+1} &=& x_k + \Delta T \dot{x}_k + \frac{1}{2} \Delta T^2 \\
\dot{x}_{k+1} &=& \dot{x}_k + \Delta T \ddot{x}_k \\
\ddot{x}_{k+1} &=& \ddot{x}_k + \eta_k
\end{align}
$$
with similar expressions for $y$.
That means
$$
\mathbf{A} = \left[
\begin{array}{cccccc}
1 & \Delta T & \frac{1}{2} \Delta T^2 & 0 & 0 & 0\\
0 & 1 & \Delta T& 0 & 0 & 0\\
0 & 0 & 1& 0 & 0 & 0\\
0 & 0 & 0 & 1 & \Delta T & \frac{1}{2} \Delta T^2\\
0 & 0 & 0 & 0 & 1 & \Delta T\\
0 & 0 & 0 & 0 & 0 & 1
\end{array}
\right]
$$
and
$$
\mathbf{B} = \left [
\begin{array}{cc}
0 & 0 \\
0 & 0 \\
1 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 1 \\
\end{array}
\right ]
$$
Then the process noise covariance is $\mathbf{B}\mathbf{B}^T \Sigma_\eta$ where $\Sigma_\eta$ is the covariance of $\eta_k$. This is called the transition_covariance in pykalman.
If I do that and take some guesses for the process and measurement noise scales, then I get the result below.

As you can see, the early part of the plot gets perturbed by missing data.  I'm not quite sure the best way to deal with that.
New Code below
import numpy as np
import pandas as pd
from pykalman import KalmanFilter
from numpy import ma
import matplotlib.pyplot as plt

# Reading data
observations = pd.read_csv('Q66687.csv')
rearwheel_y = observations['rearwheel_y'].to_list()
rearwheel_x = observations['rearwheel_x'].to_list()

# Locate and clip relevant snippet of motion, removing unecessary data at the beginning and end
rearwheel_y = rearwheel_y[385: 600]
rearwheel_x = rearwheel_x[385: 600]

#rearwheel_y = np.sin(20*time) + 10

# The numpy.ma module provides a convenient way to address the issue of dropout, with masked arrays.
# When an element of the mask is False, the corresponding element of the associated array
# is valid and is said to be unmasked.
# When an element of the mask is True, the corresponding element of the associated array
# is said to be masked (invalid).
rearwheel_x = ma.masked_values(rearwheel_x, 0)
rearwheel_y = ma.masked_values(rearwheel_y, 0)
measurements = np.array([rearwheel_x, rearwheel_y])

# initial_state_mean
Xbar = np.array([rearwheel_x[0],0,0,rearwheel_y[0],0,0]).transpose();
# initial_state_covariance
P0 = 100*np.identity(6)
n_timesteps = len(measurements[0])
n_dim_state = 6
filtered_state_means = np.zeros((n_timesteps, n_dim_state))
filtered_state_covariances = np.zeros((n_timesteps, n_dim_state, n_dim_state))

# Kalman-Filter initialization
DeltaT = 0.01668335
A = np.array([[1, DeltaT, 0.5*DeltaT**2, 0, 0, 0],
              [0, 1, DeltaT, 0, 0, 0],
              [ 0, 0, 1, 0, 0, 0],
              [0, 0, 0, 1, DeltaT, 0.5*DeltaT**2],
              [0, 0, 0, 0, 1, DeltaT],
              [0, 0, 0, 0, 0, 1]])
B = np.array([[0, 0],
             [0, 0],
             [1, 0],
             [0, 0],
             [0, 0],
             [0, 1],])
SigmaEta = np.array([[1,0], [0,1]]) # Assume eta covariance is unit...
#Q = B @ B.transpose() * 0.0001
Q = np.identity(6) * 0.00001
R = np.array([[1,0],[0,1]])*0.00005
H = np.array([[1,0,0,0,0,0],[0,0,0,1,0,0]])
kf = KalmanFilter(transition_matrices = A,
                  observation_matrices = H,
                  transition_covariance = Q,
                  observation_covariance = R,
                  initial_state_mean = Xbar,
                  initial_state_covariance = P0)
# iterative estimation for each new measurement
for t in range(n_timesteps):
    if t == 0:
        filtered_state_means[t] = Xbar
        filtered_state_covariances[t] = P0
    elif t != 0:
        filtered_state_means[t], filtered_state_covariances[t] = (
        kf.filter_update(
            filtered_state_means[t-1],
            filtered_state_covariances[t-1],
            observation = measurements[:,t])
        )

# plot of the resulting trajectory
plt.figure()
plt.plot(rearwheel_y, 'k.', label = 'Observations')
plt.plot(filtered_state_means[:, 3], "y-", label="Filtered Positions", markersize=1)
plt.ylim(180, 200)
plt.grid()
handles, labels = plt.gca().get_legend_handles_labels()
by_label = dict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys())

Older version below
OK, so if you change your Kalman filter to something like this:
kf = KalmanFilter(transition_matrices = np.array([[0.1,0.9], [0, 1]]),
                  observation_matrices = np.array([[1., 0]]),
                  transition_covariance = np.array([[1, 0],[0,1]]),
                  observation_covariance = np.array([[100]]),
                  initial_state_mean = np.array([[0,0]]),
                  initial_state_covariance = np.array([[0.13, 0], [0, 0.13]]))

you get:

You can change the observation_covariance up or down to get smoother or more jagged tracks. The smoother tracks will have more delay.

Original post
To get anything more sensible out of it, you'll have to choose a better model.  The model you chose before your recent edit had that the measurement changed as a simple random walk.  If that's all you know about how the measurement changes, then what you have is as good as it gets.
If you know more, for example if the measurement has some inertia, then you will need to choose different parameters, especially for your transition matrix.
Without understanding a bit more about what your measurements are, I won't be able to help much. Have a look at this answer and see if that helps.  That answer looks at including position and velocity into the state, so that the transition_matrices parameter is a $2\times 2$ matrix.
I couldn't load your data sensibly, so I just generated an offset sinusoid and changed the KF parameters.  This can yield something like the plot below.  However, I'd need to know more about your signal model to get it to work sensibly for your situation.

Updated Code Below
import numpy as np
import pandas as pd
from pykalman import KalmanFilter
from numpy import ma
import matplotlib.pyplot as plt

# Add time
time = np.linspace(-np.pi, np.pi, 201)

# Locate and clip relevant snippet of motion, removing unecessary data at the beginning and end
#rearwheel_y = rearwheel_y[385: 600]

rearwheel_y = np.sin(20*time) + 10

# The numpy.ma module provides a convenient way to address the issue of dropout, with masked arrays.
# When an element of the mask is False, the corresponding element of the associated array is valid and is said to be unmasked.
# When an element of the mask is True, the corresponding element of the associated array is said to be masked (invalid).
# rearwheel_y = ma.masked_values(rearwheel_y, 0)

# time step
dt = 1
# initial_state_mean
X0 = 0 #rearwheel_y[0]
# initial_state_covariance
P0 = 0
n_timesteps = len(rearwheel_y)
n_dim_state = 2
filtered_state_means = np.zeros((n_timesteps, n_dim_state))
filtered_state_covariances = np.zeros((n_timesteps, n_dim_state, n_dim_state))
# Kalman-Filter initialization
kf = KalmanFilter(transition_matrices = np.array([[0.01,0.99], [0, 1]]),
                  observation_matrices = np.array([[1., 0]]),
                  transition_covariance = np.array([[1, 0],[0,1]]),
                  observation_covariance = np.array([[1000]]),
                  initial_state_mean = np.array([[0,0]]),
                  initial_state_covariance = np.array([[0.3, 0], [0, 0.3]]))
# iterative estimation for each new measurement
for t in range(n_timesteps):
    if t == 0:
        filtered_state_means[t] = X0
        filtered_state_covariances[t] = P0
    elif t != 0:
        filtered_state_means[t], filtered_state_covariances[t] = (
        kf.filter_update(
            filtered_state_means[t-1],
            filtered_state_covariances[t-1],
            observation = rearwheel_y[t])
        )

# plot of the resulting trajectory
plt.figure()
plt.plot(rearwheel_y, 'k.', label = 'Observations')
plt.plot(filtered_state_means[:, 0], "y-", label="Filtered Positions", markersize=1)
plt.ylim(8, 12)
plt.grid()
handles, labels = plt.gca().get_legend_handles_labels()
by_label = dict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys())


Audio frequency modulation algorithm I have a digital audio signal (float samples). I need to apply the Vibrato effect to a signal. "Vibrato" is no common term, but it is the same as the frequency modulation.
So, I need to modulate my audio signal by a low-frequency sinusoid. I need to set two parameters for the modulation: "amplitude" and period. In this case, "amplitude" of freq modulation is a maximum shift of signal frequency, for example 10 Hz. Period is just a period of the sinusoid.
So, I need a C/C++ function for such modulation.
The related article, but I cannot quite understand it: https://christianfloisand.wordpress.com/2012/05/01/shaking-it-up-with-vibrato-2/

You need to build a time varying delay, where you can modulate the delay amount over time.

The peak delay modulation is a function of your maximum desired frequency shift and the modulation frequency.

This is not trivial since it will require fractional sample delays with some kind of interpolation algorithm. You can't round to the nearest integer delay without "zipper" noise. The best interpolation method depends on the specific requirements of your application (tolerance for noise, input spectrum, MIPS & memory available, etc.)


Is the sampling period arbitrary for a bilinear transform, and why? I have been given the specifications for a digital highpass filter (stopband, passband, stopband attenuation and maximum passband ripple).
I am expected to design a prototype lowpass filter in the analogue domain, use frequency transformation to convert it to a highpass, then use bilinear transformation to convert it into a digital filter. I am using a Butterworth filter for prototype.
However I have not been given the sampling frequency(or period) for pre-warping. I have seen textbooks take period = 1 as it is arbitrary, however I don't understand why. Any confirmation that this is okay, and explanations will be appreciated!

The sampling frequency is not arbitrary and important to specify for pre-warping. This is clear when reviewing how each analog frequency in the bilinear transform is mapped one-one to a digital frequency according to:
$$f = \frac{F_d}{F_s} =  \frac{1}{\pi}\tan^{-1} \bigg(\pi\frac{F_a}{F_s}\bigg) \tag{1} \label{1}$$
Where:
$f = F_d/F_s$ is the normalized digital frequency in units of cycles/sample,   $F_a$ is an analog frequency in Hz
$F_s$ is the sampling rate in Hz
Observe that $tan^{-1}(\infty) = \pi$ and thus the Bilinear Transform "warps" the entire analog frequency domain from $F=0$ to $F= \infty$ to the digital frequency domain of $f=0$ to $f=1$ (where $f=1$ is the sampling rate in units of normalized frequency), but the specific mapping of each analog frequency over that range is dependent on the value of $F_s$ that is used, as given by \ref{1}.
Normalizing the discrete-time frequency domain is commonly done such that the resulting spectrum associated with DC to the sampling rate extends from $f =0 \text{ to } 1$ cycles/sample, which may be the source of the OP's confusion.


Linear vs Non linear distortion I'm learning about distortion and I got stuck at

"What is the difference between linear and non-linear distortion?"
(Also in frequency spectrum i.e. Amplitude and phase spectrum)

The Wikipedia article on "distortion" is actually pretty good. https://en.wikipedia.org/wiki/Distortion

Partially it's about semantics and how you exactly define the term. By most definition "linear distortion" is used equivalent to "non ideal frequency response" where magnitude and/or phase are a non-trivial function of frequency.

Non-linear distortion is caused by any type of mechanism that violates linearity. Typical examples are limiting or clipping or compression.

An LTI does not cause any type of non-linear distortion, however they typically do cause linear distortion if defined as "non-flat frequency response"


Determining invertibility of a system Given a system which is either non-linear and/or time variant, Is there any systematic way to determine if this system is invertible?
For example how to determine invertibility of this system?
$$\int_{-\infty}^\infty x(5\tau)d\tau$$

In general there is no systematic way and you simply have to analyze the given system. In the case of the system in your question, it's easy to see that it can't be invertible, because the output is just a constant, namely the integral over the input function (assuming this integral exists). There are infinitely many functions which will result in the same value when integrated from $-\infty$ to $\infty$.


What does an image of Fourier Transformation of an image tell us? First time studying image processing...
I just don't understand what does fourier transformed image of an image describe?
For example consider given following pictures, The first one is the image, and the second one is the fourier transformation of the image.

Now my question is:
By given the fourier transformation image, what am i suppose to comprehend from that?

I would really appreciate any help here, i just cannot proceed with my studies without understanding this.

Thanks in advance.

The Fourier transform image is a convenient representation to perform other functions on the images such as filtering. By itself, there's probably not much you can intuit by just looking at it.

In the example you posted you can see how most of the energy in the Fourier image is near the center. This means that most of the energy in the image is in low frequencies. In other words, there are many areas of the picture with smooth texture like her shoulder, the hat and the blurred background.

The black areas in the fourier image are away from the center. These represent the high frequencies, i.e. rough texture on the image or portions of the image that change rapidly. The high frequencies are black in the fourier image, which means that there is not a lot of high frequency energy in the Lena image. If you took a picture of static on an old tv screen then it would have a lot of high frequency energy in the image. There would be no dark regions in the fourier image.

Each pixel in the fourier image corresponds to a 2d-wave that is used to make up the original image. There is no spatial correspondence between the regions of the fourier image and the original image. Each pixel in the fourier image corresponds to a 2d-wave across the whole of the original image.


Normalizing vectors to be the same size to retain PSD data Context:
PSD=periodogram(y) in MATLAB gives the PSD (power spectral density) estimate of a sound-source signal. I want to find the local maximum within 3 different x ranges of PSD.

How do I make a vector that maps 1:1 to the PSD -y value vector?
The x values that MATLAB generates are Hz - frequency values, so I don't want to change the values per se, but create a vector that goes from x=0:3500, but have the same number of elements as the PSD vector, PSD= 32769*1 double.

How do I do this while retaining the same data?
    PSD plot. acquired with plot(PSD) <-- x values generated by MATLAB

  [1]: https://i.stack.imgur.com/PiYGB.png

    %% formantdaddy will try to do all the steps of the process
% 1. grab sound from folder
% 2. get PSD = Power Spectral Density
% 3. look for local maxima within prefined ranges
    % 200 - 800 Hz
    % 800 - 1800 Hz
    % 1800 - 3500 Hz
%store data into a struct

clear all;
close all;
clc;
%% audioread .wav file
[y Fs]= audioread('100-daddy1.wav');
%% y=source signal
sourceFig=figure(1);
plot(y);
xlabel('milliseconds'); ylabel('amplitude'); title('spectrogram = voice source signal');

%% PSD - Power Spectral Density <-- peaks here should give formants
filterfcn=figure(2);
PSD=periodogram(y); %calculate the power spectral density of the source signal
plot(PSD);
ylabel('magnitude || intensity of signal'); xlabel('frequency in Hz'); title('PSD of Source');
xlim([0, .35*10^4]);
%% F1 range
subset1=zeros(1, 601); %create a vector to store PSD values within the range of interest
for x=200:800;
    subset1(x)=PSD(x);% input in the values of interest
end

for x=200:800;
    Y_max=max(subset1);
    F1=x;
end
disp(F1);
disp(Y_max);
% find the max value within the range
%PROBLEM: wanted to find the x value which corresponded to the max y <-- but this actually is the max y value
%PROBLEM: x vector not the same size as PSD values within the range, so many points not considered.

%% F2 range
subset2=zeros(1, 1001);
for x=800:1800;
    subset2(x)=PSD(x);
end

for x=800:1800;
    Y2_max=max(subset2);
    F2=x;
end
disp(F2);
disp(Y2_max);
%% F3 range
subset3=zeros(1, 1701);
for x=1800:3500;
    subset3(x)=PSD(x);
end
for x=1800:3500;
    Y3_max=max(subset2);
    F3=x;
end
disp(F3);
disp(Y3_max);
%% store values into a struct
Ftable=struct("FileName", "100-daddy1.wav", "F1_Hz", F1, "F2_Hz", F2, "F3_Hz", F3);

For function call [pxx, f] = periodogram(y_soundsource_data,window,nfft, Fs);
  you can get the pxx = PSD points= vector of y values, and f= frequency points tied 1:1 to the pxx values.

... this requires that you specify nfft & Fs

You can get y_soundsource_data and Fs<-- sampling frequency, from the audioread function. file[y Fs]= audioread('100-daddy1.wav');
window and nfft are not understood by the periodogram function UNLESS you define them.

y_soundsource_data and Fs are understood, because audioread returns the data metrics.
nfft - refers to the number of points used in the DFT. A larger nfft value will lead to a better estimation. nfft MUST be an integer. You can't have half a point.
window = vector of the size of the y_soundsource_data. 1:length_y_soundsource_data
both window and nfft need to be defined before you call the periodogram function.

keep getting errors with the [pxx,f] = periodogram(_,fs) version of the call. mathworks.com/help/signal/ref/periodogram.html


Measured signal from tachometer has ripples, why and how to minimize? I don't know if it would be straightforward enough to chalk it up to this being a real system. Thus there's friction and physical mechanics of gears and teeth pulling.

Well, regarding the why: don't know the underlying mechanism of your tachometer; there's many, many ways to build one.

Minimization: compared to your signal of interest, these seem to be rather high in frequency. So, a reasonably designed low-pass filter will help. Don't forget to account for group delay, should you need to feed back the resulting signal in e.g. a control system.


Algorithm to pan audio I'm looking for a way to pan the PCM audio I have.
At the moment the audio is in a stereo format and I would like to pan it to the desired side (or more/less to a side).

I've recently figured out how to do mixing here so I guess it would be something similar (adding stuff to each other and such). What I was thinking is to add the audio data from one side to the other side.
For example if I were to pan hard left, I'd get the data from the right channel and add it to the left channel. I would then set the right channel to a constant value so that it would be silent.

Am I thinking in the right way or am I missing something obvious?

Thanks in advance!

The basic technique to place a mono source in stereo is called constant power panning.

If you want to place a mono source at angle $\theta$ you can just use $A_\mathrm{amp}$ and $B_\mathrm{amp}$ as amplitudes for your channels:

$A_\mathrm{amp} = \frac{\sqrt{2}}{2} (\cos{\theta} + \sin{\theta})$

$B_\mathrm{amp} = \frac{\sqrt{2}}{2} (\cos{\theta} - \sin{\theta})$

Typically $\theta$ ranges from $-45^{\circ}$ to $45^{\circ}$ with $0^{\circ}$ being the center.

Hope, that helps.

Reference: Roads, Curtis (1996). The Computer Music Tutorial. Cambridge: The MIT Press, pp. 457–461.

In addition:

You can put it in code directly. E.g. if you use Python and you have mono audio in a numpy array, you can convert it into stereo by a given angle like this:

def panner(x, angle):
    """
    pan a mono audio source into stereo
    x is a numpy array, angle is the angle in radiants
    """
    left = np.sqrt(2)/2.0 * (np.cos(angle) - np.sin(angle)) * x
    right = np.sqrt(2)/2.0 * (np.cos(angle) + np.sin(angle)) * x
    return np.dstack((left,right))[0]

You can use it like panner(np.array([1,2,3]), np.radians(20)). (Of course [1,2,3] is a quite pointless audiobuffer.)


How to calculate the zero state response? In the following question and answer I kind of get it how it is calculated at surface but I don't get it fully. Can someone please explain the answer in depth?

The answer to your question is essentially contained in Dilip's comment, but I guess as a beginner you need some more help. First, it's important to realize that the response of an LTI system to a scaled and shifted impulse
$$x[n]=a\delta(n-k)$$ is simply a scaled and shifted version of its impulse response:

$$y[n]=ah[n-k]$$

Furthermore, the response to a sum of signals equals the sum of the responses to the individual components of the input signal. Consequently, the response to the given input signal

$$x[n]=2\delta[n]-3\delta[n-2]$$

is

$$y[n]=2h[n]-3h[n-2]\tag{1}$$

Now you just have to express the impulse response $h[n]$ in terms of $\delta$-impulses. By noting that
$$u[n]-u[n-3]=\delta[n]+\delta[n-1]+\delta[n-2]$$ you get

$$h[n]=\delta[n]+\delta[n-1]+4\delta[n-2] \tag{2}$$

Plugging (2) into (1) gives the desired result.


Why does the Y component in YCbCr represent luma/intensity? As far as I know, the RGB to YCbCr conversion is just a linear transformation, yet the Y channel successfully captures the intensity value of each pixel. How?

In other words, where do the following magic numbers come from?

$$
\begin{align}
E'_Y &= 0.299E'_R + 0.587E'_G + 0.114E'_B \\
E'_R - E'_Y &= 0.701E'_R - 0.587E'_G - 0.114E'_B \\
E'_B - E'_Y &= -0.299E'_R - 0.587E'_G + 0.886E'_B
\end{align}
$$

Focus on the first equation for EY. Back in the day when color television was being developed, the color signal had to be compatible with black and white TVs and vice versa. So the compatible brightness signal (luma Y) has to be calculated from the three primary color signals (R, G B) for transmission. Human visual system does not perceive brightnesses of different wavelenghts equally. So the actual coefficients depend on which specific RGB phosphors are used on the CRT, the spectrum that each phosphor emits, their brightness differences and how human visual system perceives them. See Color FAQ by Charles Poynton for better explanation.


How do I calculate the noise variance in this problem? In a paper about radar, the transmitted signal expressed in low-pass equivalent form is given by $\sqrt{E}s(t)$, where $E$ is the transmitted energy, $s(t)$ satisfies
$$
\int_T | s(t) |^2 dt = 1,
$$
and $T$ is the pulse duration.

Here
$$
s(t) = \sqrt{\frac{1}{T}} \exp(j2\pi f_c t)
$$
and the received signal is
$$
r(t) = \alpha s(t) + n(t)
$$
and the SNR in this paper is
$$
\frac{E |\alpha|^2}{\sigma^2}
$$
where $\alpha$ is the reflectivity, and $\sigma^2$ is the variance of the zero-mean Gaussian noise process $n(t)$.

So if I want to generate $r(t)$ in Matlab with SNR = 10, I sample $s(t)$ with a sample period $T_S$ so I get
$$
s(nT_s) = \sqrt{\frac{1}{T}} \exp(j2\pi f_c nT_S)
$$

My question is:  how do I calculate the noise variance, $\sigma^2$ from the above so that I can calculate:
$$
r(nT_S) = \alpha s(nT_s) + \sigma {\tt randn(size(}s{\tt ))}
$$

??

Since the signal $r(t)$ is expressed as the low-pass equivalent of a bandpass signal, in general it will be complex. That means that the noise is also complex.

White Gaussian noise has a power spectral density commonly denoted by $N_0/2$. This means that if you filter the noise with an ideal filter of bandwidth $B$, the noise at the filter's output is $\sigma^2=BN_0$.

Note that this is the power "per dimension". That is, the in-phase and the quadrature components of the low-pass signal $r(t)$ each suffer from noise of power $\sigma^2$.

You define the SNR as $E|\alpha|/\sigma^2$ and you need it to be equal to 10. One way to do this is to configure your simulation so that $E|\alpha|=1$, and then choose $\sigma=\sqrt{0.1}$. If you later need a different SNR, you can just change one value ($\sigma$) in your simulation.

You can generate the complex noise as follows:

N = length(s);
noise = sqrt(0.1)*(randn(1,N)+1i*randn(1,N));
r = alpha*s + noise;

Note that in this setup, the total noise power is $2\sigma^2$. The in-phase noise and the quadrature noise each have power $\sigma^2$.


How can I know if these two signal are orthogonal to each other So the question is to determine if this statement is true or false. The question can be seen below.

I know that two signals are orthogonal if their inner product is equal to 0. I know how this calculation is done with vectors but not with signals. Would be really helpful i someone could explain this to me. Is it also possible to just directly see it in the two graphs if they are orthogonal to each other?

The inner product of two continuous-time signals is defined by

$$\langle\mathbf{s}_1,\mathbf{s}_2\rangle=\int_{-\infty}^{\infty}s_1(t)s^*_2(t)dt\tag{1}$$

where $*$ denotes complex conjugation. In your example the conjugation is irrelevant because both signals are real-valued. So you just need to multiply the signals and determine the area of the result.

With $(1)$ it should be easy for you to show that the two given signals are indeed orthogonal.


what is the output of BoW after an image has been trained with SIFT algorithm and k-means SIFT algorithm provides a 128 dimensional feature vector that is used for image classification.When all the interest points(key points) are taken together and K-means clustering is applied,the image is then clustered into many clusters.My question is how does this clustering/grouping of feature points together help in creating a dictionary(of equivalent image words).Does the output after applying this BoW algorithm yield a single vector for an image.

You essentially got it right: the final purpose of the whole BoW clustering algorithm is to somehow produce a single image descriptor for every image.

In case of BoW clustering (either K-means, or hierarchical K-means, or some other clustering), this image descriptor is a histogram of visual words for that image often normalized by the number of local features for each image.

Let me offer a step-by-step explanation of the process:

first you get all the local features from all the images (or just from your learning dataset), and cluster all those features by a clustering algorithm of choice
after completing the clustering, for every cluster we calculate the "representative member" by some kind of mean or average over all the local feature samples from that cluster

This "representative member" is actually your visual word -- it's like the "core" of all the features of the cluster (if you would be working with actual words, the samples "work", "working" and "workaround" might all have a common representative member "work").

This step might actually be skipped -- we don't care about exact visual word, as long as we can distinguish between the clusters (e.g. knowing how to classify a new local feature in to "cluster 1", "cluster 2" or "cluster 3" is enough).
now, for every image in our database we construct the "image descriptor" by making the histogram by clusters for that image.

E.g., if there is 4 clusters in final clustering, and our image has 4 features in cluster 1, 10 in cluster 2, 4 in cluster 3 and 2 more in cluster 4, the descriptor might be: 4|10|4|2.
finally, we can normalize the histograms

Following with the same example from above, the image had a total of 20 features. We could thus normalize the 4|10|4|2 histogram by dividing everything by 20 and getting: 0.2|0.5|0.2|0.1.

With this, you ended your "database description" -- you have a numeric descriptor for all your images. Now, when a new image arrives, you do the following:

extract local SIFT (or any other kind) of features
determine for all the local features of the new image to which cluster they belong
make a histogram for the new image (+ normalize histogram)
compare the histogram of the new image to the histograms of database images

You have methods other than BoW to create the "image descriptors". One alternative to BoW is VLAD (Vector of locally aggregated descriptors), which is a bit more compact but has the same purpose.


Is it possible to hack camera calibration without having access to the camera? Many times we have a stream of video to process without access to camera. Having access to camera matrix would be beneficial for various processing techniques. Is it possible to hack camera calibration without having access to the camera?
I have a stream of video from a single camera mounted on a moving car recording the road (Hence multiple parallel lines on the ground plane, corners from lane markers but no circles). I want to create a top-down view of this but I do not have access to the camera. Is it possible? If so how?
I understand from Learning OpenCV: Computer Vision with the OpenCV Library book that I need the following matrices:

Intrinsics, and
Distortion

.. to cv2.undistort() the image, compute homography cv2.getPerspectiveTransform() and cv2.warpPerspective() to finally get the topview.
Question

How can I compute/approximate/guess Intrinsics or Distortion Matrices/parameters? Are all parameters important?
Would it be okay to copy parameters from other cameras (like OpenCV source code samples)?

Prior Research:

OpenCV Python Camera Calibration Tutorial - Requires Access to the Camera
Camera AutoCalibration - Gives hope "calibration may be obtained if multiple sets of parallel lines or objects with a known shape (e.g. circular) are identified"
Attempting to understand camera calibration related answers[1][2][3] on SO trying to find answers to my problem.

Any ideas? Appreciate your help in advance!

Update 1: Perspective Transform Experiment
I had attempted using getPerspectiveTransform(pts1,pts2) and dst = cv2.warpPerspective(img,M,(x,y), flags=flags) to match four points of input image to get a sort of top-down view. But I am not sure how to handle the distortions:
Selecting Points: First, I zoomed in on the input and tried to precisely select matching exterior points on the lane markers to create the Homography Matrix and previewed the perspective transformed image with warpPerspective. I saw the lanes were distorted but didn't know how bad. To get an idea I chose points further out by delta (50px) flat on each end. This is what I get:

Can I fix the distortions without having access to the camera? Is there any other way to fix this.
Few input images to play with:

Update 2: Non-Parallel (Scattered) Optical Flow in Top-Down View
Is this due to distortion or something else?:

If I understand correctly, you don't need the intrinsics or extrinsics to achieve that, if a top-down view is all you want. You could basically define 4 points on your parallel lines and then warp the entire image into a canonical view (say $\{\{0,0\}, \{480,960\}\}$).

To do that in OpenCV, all you need to do is compute the homography using findHomography as described here. Then, use warpPerspective (look here) to warp your images.

The only thing you need to care about is the selection of the dimensions of the fronto-parallel view that you are after. It should more or less match the aspect ratio of the original image region.

Regarding the second part of your question: Yes. First of all, let me remind you that many works such as this choose to ignore the distortion estimation. But, I can think of two ways to accomplish that:

Use the road lanes itself to correct the distortion: In the fronto-parallel view, create another perspective warp so that the lanes themselves would be vertical. In this setting, think about the lane lines as the vanishing lines in the image (as in a perspective view, they always do). From that, you can even calibrate your camera. Look here, here and here for example.
The cars should be following more or less a straight/smooth trajectory. If you compute the motion of all the cars, they should give you similar motions, parallel to the lanes. That information is again useful in constraining and generating a transformation, which would warp your image to a reduced distortion one.


The precise meaning of 'coefficient' in signal processing when I read papers on signal processing, I often see the term 'coefficient' being used to refer directly to 'samples', i.e. "We can reduce the number of coefficients, by simply subsampling the even samples." or in presentations on the lifting scheme (wavelets), it seems the output of the low pass filter is also called coefficients.

To me, the output of a low pass filter is a 'smoothed' version of the original input sequence. Further, to me, a coefficient is e.g. a, b in  y[n] = ax[n] + bx[n-1]. In this sense, the coefficients a and b certainly is a representation of the signal-- but i would not call y[n] the coefficients.

I am probably very confused about the terminology. Could someone clarify this with some intuition?

There's really not much I can clarify for you here:

A coefficient really is just a factor in series. Whether that series is e.g. a Fourier series, a sum over powers of 2, filter delay taps or just a series of sampling instants can't be said. Context!


sketching magnitude of frequency response of H(z) I am trying to plot the magnitude of H(z), I got it to factore, and would like to sketch its plot of magnitude. But I am have trouble evaluating the function.
$$
H(z)= \frac{(1-2z^{-1})(1+0.5z^{-1})(1+0.9z^{-1})}{(1-z^{-1})(1-0.7jz^{-1})(1+0.7jz^{-1})}
$$

I know we replace z with $e^{jw}$ and try to evaluate it at various w such as pi, pi/2 etc. But that is really complicated and when I tried it I couldn't get the right answer as the text book. Is there an easier way to evaluating magnitude and phase by hand?

Thanks

Since you have the transfer function factored into poles and zeros, it's easiest to draw the pole-zero diagram to get an idea of the overall shape of $H(z)$. The first thing you notice is the pole at $z=1$, i.e. at DC. So your filter is unstable and $|H(e^{j\omega})|$ goes to infinity at DC. Furthermore there's a complex conjugate pole pair on the imaginary axis, i.e. at half the Nyquist frequency. So you can expect a small bump around that frequency. Finally, I would calculate the value of $H(z)$ at Nyquist, i.e. at $z=-1$:

$$H(-1)=\frac{3\cdot 0.5\cdot 0.1}{2\cdot 1.49}\approx 0.05$$

A more detailed sketch than that is indeed difficult, but that's probably not the idea of the exercise.


Why does altering path length in radio transmission alter the phase of the received signal? I do have a very naive yet bothering question:
Assume I send a sinusoid signal

\begin{equation}
A*cos(2\pi ft+ \phi_0)
\end{equation}

over a radio channel. Assume there is only a LOS-path for each receiver antenna and no interference or other fading effects.

Next, I receive the signal at two different locations with different distances from the receiver. Will the phase of the received signals be different for the two sites?

In order to deepen my understanding:

Let's assume I record my received signal and I set $t=0$ for each receiver individually when the signal arrives at the respective receiver.
Will the phase of the received signal at each receiver be different for different distances? Obviously, for one observation instant (time on reference clock equal for all receivers) the phases vary at the different sites. However, what I mean is:
Does the signal phase vary for the situation I explained above (individual receiver clocks)?
This question boils down to the understanding of the physical reality of an EM wave propagating. Assume I follow the wave and I 'look' at the wavefront at different times; will I see the same phase at all times? Or does the 'front' of the signal wave oscillate in time and space?
I really hope that does not sound too stupid; it's really a question that bothers me...

Yes indeed it will be (although we could argue that it may not be necessarily for all distances since the phase is cyclical!).

In free space the signal propagates at the speed of light, therefore this sets the wavelength in distance based on the frequency transmitted according to:

$$\lambda = c/f$$

Where $c$ is the speed of light in meters/second (or other units of choice), and $f$ is the frequency of the signal in units of  Hz if time is in units of seconds. A fixed tone at a given frequency with have a length in phase of $2\pi$ over the distance of one wavelength.

If the channel has a dielectric constant higher than free space given by $\epsilon_r$, the wavelength will be decreased by $1/\sqrt{\epsilon_r}$.

For example the wavelength of a 1 GHz signal in free space (or air which has a dielectric constant that is very close to free space) is approximately $3e8/1e9 = 0.3$ meters, while over a cable with a polyethylene dielectric ($\epsilon_r = 2.3$) would be $3e8/(1e9\sqrt{2.3}) = 0.198$ meters.

For this case of a 1 GHz signal in free space given a wavelength of $0.3 m$ (corresponding to the phase over one wavelength of $2\pi$ radians), if you receive the direct signal with no multi-path at two different locations with a $1$ cm path difference between the two, the expected phase difference will be $0.21$ radians.

As a follow up to the OP's question; if you set $t=0$ upon arrival of the wavefront, then the phase everywhere at the new $t=0$ will be equal to the phase when the signal was initiated at the transmitter. Although this is somewhat of a meaningless result: it is equivalent to measurement the distance of two objects and resetting the distance to 0 once you determine the objects length.


Finding baseline (reference) of nosiy data and lag/lead between 2 signals I have 2 datasets which vary with time as shown in the figure - Sa – Red  ,  Sb – BLACK . I am doing the analysis in R. The entire dataset and a segment of data from 30 - 31 s are shown below -
 -

As you can see from the figures, the datasets appear similar but one dataset  leads (or lags) the other dataset in time. The downward peaks appear to be slightly shifted in time for one of the signals wrt the other.  I would like to determine the following –

The baseline value (reference) for the two signals –
Is there a statistical method to find out this? Is there some signal processing techniques that need to applied before determining the baseline? For Sa, the baseline is more noisy and shifts up and down – so it seems difficult to assume a constant reference line
The time lag / lead between these two signals –
I tried to make the reference line constant by removing data points above an approximate baseline and replacing it with the approx baseline value. Then I found out values sa1, sa2, sa3, sa4, ..etc.  Similarly, I found sb1 ,sb2, sb3, sb4,… etc.

2a. Is there a way I can map sa1 to sb1, sa2 to sb2..and so on? The idea is to calculate (sa1-sb1), (sa2-sb2),… and find out the average time lag/ lead for the two signals over the entire dataset. I understand that this method is very approximate, since the resolution is lost while assuming a baseline based on just eyeballing.

2b. Is there another method for calculating the time lag/ lead? Or the exact times sa1,sa2.. etc  at which the signal drops below the reference line and comes back to its reference line?    Will a  cross correlation work? How do I use it for noisy signals?

3.The frequency of the downward peaks

You can try the function 'ccf' in R.
The way I would setup your data is:
df=
time | red | black | red diff | black diff |

** make sure your time interval's are constant. i.e every 5 seconds/10 seconds etc

red diff and black diff being the difference between the adjacent values in the red and black data columns, respectively.
i.e. diff(df$red) = red diff

Then do a ccf on the 2 difference columns and adjust 'lag.max' as necessary to view the plot that is generated. There should be a peak in the graph which is the lag (that occurs the most) between red and black.


Frequency spectrum of a sinc function I am doing one example from my book as a preparation for exam.

The assignment is:

  It is given that:
  $$\mathbb{rect}(t)=pf(t) \leftrightarrow PF(f)=2AT_0 \cdot \mathbb{sinc}(2\pi fT_0)$$
  you need to calculate frequency spectrum of:
  $$pf(t)=\mathbb{sinc}(\omega t)$$

Truthfully, I have no idea where to start. I presume that the solutions is absolute value of sinc function, because I read it from solution, but in the solution there was only diagram.

I tried to solve directly using Fourier transformation on sinc function, but I got very messy equation at the end.

So, my question is: How can I solve this assignment?

Thanks you very much!!!

EDIT:
I found this pdf, on the 7(162) side it explains what I want to do, but this is only with pictures. I want to understand it.
http://ultrasound.ee.ntu.edu.tw/classnotes/ckt2/Chapter12.pdf

Look at the forms of the functions, and remember the duality property.

I invite you the read the Wikipedia article: http://en.wikipedia.org/wiki/Fourier_transform


Which Anisotropic Diffusion Method is Best for Denoising Now? The best standard may be varying from PSNR to human's visual feeling. Any one is OK.

I have tested J.Weickert's tensor diffusion and P&M's, they all can't satisfy me.

Is there any anisotropic diffusion method now is best for denoising, especially in tensor diffusion, others like the KAD, MAD and so on is ok as well.

Thanks for reading~

I think "Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE’s" by David Tschumperlé is really good.

I started implementing it using MATLAB once, yet I didn't finish.
If you do implement it, it will be great if you updated us.

Update

Have a look on my implementation at Fast Anisotropic Curvature Preserving Smoothing.


How to encode binary sequence if the number of bits are odd to e.g. 4-PAM? (Beginner) I'm new to signal constellation and trying to understand how M-PAM works. If there is a binary sequence that has on odd number of bits and I want to send a 4-PAM signal, there will be one bit left after encoding a symbol by two bits ($k=2$). For example: the sequence $1001001$ will be encoded as $10\rightarrow3$, $01\rightarrow-1$, $00\rightarrow-3$ and the last bit will be by itself. How do I encode the last bit? What about 16-PAM? Will I need to express the sequence in more bits?

This problem is solved in the upper layers of the communication system, not at the PHY.
For example, in the data link layer (DLL), the transmission unit is a frame. Frames are designed to be compatible with the PHY, so that problems like the one you describe don't occur. Also, many DLL protocols allow for fixed-size frames to be padded, or for variable-size frames, which provide some flexibility.