Usefulness of Practical filter ouptput If a practical filter can't remove all unwanted frequency components like ideal filter, does that  mean unwanted frequency component are still present after filtering?

how can we use such distorted output in practical use ?

A real filter wont perfectly remove all unwanted frequencies. But then, a real $X_\mathrm{anything,whatever}$ won't do its job perfectly. The real world isn't perfect, you only get perfection in mathematics.

What's important though is that a filter can approximately remove the unwanted frequencies. For example if you have some interference at 100 kHz, then a simple 2nd-order lowpass filter at 20 kHz will not be able to remove that completely, but it will reduce its amplitude by more than 24 dB. If the signal you're interested in is strong enough, this smaller inteference level may then not bother you anymore. Or maybe it will still bother you and you therefore need to switch to a higher-order filter. Judging something like this and deciding which filter is most appropriate for the application is the bread and butter of linear signal processing.


kernels to compute second order derivative of digital image For an image $I$, its first order derivatives can be computed using several oprators, such as
$$K_{sobel} = \left[ \begin{array}{ccc} -1 &0 &1 \\ -2 &0 &2 \\ -1 &0 &1 \end{array}\right]$$
$$I_{x} = I * K_{sobel},\ I_{y} = I * K_{sobel}.T$$.
And to compute the second order derivatives $(I_{xx}, I_{yy}, I_{xy})$, I can apply $K_{sobel}$ twice on $I$, however, I wonder are there any kernels which can be applied once to get the second order derivatives.

Here is what I tried, I made up 2 kernels,
$$K_{xx} = \left[ \begin{array}{ccc} 0 &0 &0 \\ 2 &-4 &2 \\ 0 &0 &0 \end{array}\right]$$
$$K_{yy} = \left[ \begin{array}{ccc} 0 &2 &0 \\ 0 &-4 &0 \\ 0 &2 &0 \end{array}\right]$$
and $I_{xx} = I * K_{xx},\ I_{yy} = I * K_{yy}$, they look pretty much the same as twice-applied-sobel results, but I failed to find a good kernel for $I_{xy}$, I need your help, any idea is appreciated.

Convolution $(*)$ is associative, so

$I_{xx} = K_{sobel} * (K_{sobel} * I) = (K_{sobel} * K_{sobel}) * I$

meaning the convolution of a 1st order Sobel filter kernel with another 1st order Sobel filter gives the filter kernel for a second order filter. So e.g. a 2nd order horizontal Sobel filter kernel is:

$\left(
\begin{array}{ccc}
 -1 & 0 & 1 \\
 -2 & 0 & 2 \\
 -1 & 0 & 1 \\
\end{array}
\right) * \left(
\begin{array}{ccc}
 -1 & 0 & 1 \\
 -2 & 0 & 2 \\
 -1 & 0 & 1 \\
\end{array}
\right)=\left(
\begin{array}{ccccc}
 1 & 0 & -2 & 0 & 1 \\
 4 & 0 & -8 & 0 & 4 \\
 6 & 0 & -12 & 0 & 6 \\
 4 & 0 & -8 & 0 & 4 \\
 1 & 0 & -2 & 0 & 1 \\
\end{array}
\right)$

And for $I_{xy}$:

$\left(
\begin{array}{ccc}
 -1 & 0 & 1 \\
 -2 & 0 & 2 \\
 -1 & 0 & 1 \\
\end{array}
\right) * \left(
\begin{array}{ccc}
 -1 & -2 & -1 \\
 0 & 0 & 0 \\
 1 & 2 & 1 \\
\end{array}
\right)=\left(
\begin{array}{ccccc}
 1 & 2 & 0 & -2 & -1 \\
 2 & 4 & 0 & -4 & -2 \\
 0 & 0 & 0 & 0 & 0 \\
 -2 & -4 & 0 & 4 & 2 \\
 -1 & -2 & 0 & 2 & 1 \\
\end{array}
\right)$

If you want a 3x3 kernel, you can use the 1st order derivative filter $\left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)$, which leads to the second order filters:

$\left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)*\left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)=\left(
\begin{array}{ccc}
 1 & -2 & 1 \\
 2 & -4 & 2 \\
 1 & -2 & 1 \\
\end{array}
\right)$

and

$\left(
\begin{array}{cc}
 -1 & -1 \\
 1 & 1 \\
\end{array}
\right)*\left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)=\left(
\begin{array}{ccc}
 1 & 0 & -1 \\
 0 & 0 & 0 \\
 -1 & 0 & 1 \\
\end{array}
\right)$


Is Mask R-CNN a good fit for what I want to do? I'm looking at isolating the foreground from a set of photos, and I think I want to use instance segmentation so that I can try to gauge how well-focused the main subject of each photo is.
I have rudimentary experience of Yolo/darknet (but want to move beyond simple bounding boxes). I've never used python, and am mainly a PHP developer. Happy to learn whatever's needed though, but I don't want to spend a day messing around only to find I've picked the wrong tool for the job, or that it's far too complex and I can't get my head around it!
I have a server with a 1080 GPU in it. Speed is kind of important, but so long as it's less than a couple of seconds per photo, that should be ok. Is that feasible? The photos are 10MP+ in size, but I'm guessing it's standard practice to size down to (say) 1MP first, or whatever.
I have millions of photos that I want to run this against, to help weed out poor-quality shots.
Where should I be looking to get started? I see a lot of mentions of Mask R-CNN, but I don't know if that's a good fit or if there's better/faster/easier options out there right now?
Furthermore, should I use existing yolo/darknet bounding boxes, and somehow import them into Mask R-CNN et al, or is it better to start afresh and let Mask R-CNN deal with all of it?
Thank you!

I'm looking at isolating the foreground from a set of photos

This is called semantic segmentation. Instance segmentation is where the network actually can tell the difference between different instances of the same type. For example, a few apples on a table will be instance segmented and labels the three apples individually as apple 1, 2, and 3. Semantic segmentation is an easier problem than instance segmentation and they are the same if there is only one instance of the object.

Where should I be looking to get started?

I do not have experience with Yolo/darknet but have used Mask RCNN. Mask RCNN is written in Python using Tensorflow/Keras and the base version of the code can be found here: https://github.com/matterport/Mask_RCNN. It is fairly straight forward to get running out-of-the-box with little changes needed. Mask RCNN does instance segmentation but I recall it was straightforward to add a background class.
I'd suggest choosing a tool that is meant for the problem instead of hacking a tool to match your problem. You can look into the semantic segmentation network called DeepLab, https://github.com/sthalles/deeplab_v3. This also provides a base version of the code written in Python using Tensorflow/Keras. If you pick either of these tools, I assure you that you will not waste time and believe the setup to be minimal as they both work as is.


Autocorrelation strange results? I'm trying to do autocorrelation with numpy in python but i'm getting strange results :

For a simple 440Hz sine wave sample :

This is the result of autocorrelation :

Is this normal ? Maybe I'm not understanding everything

UPDATE :

Here is my code :

def autocorrelation(signal):
  res = np.correlate(signal, signal, mode='same')
  res = res[len(res)//2:]
  plt.plot(signal)
  plt.show()
  plt.plot(res)

  plt.show()

I think you have some overflow/underflow issues? At about 32768, there are signal glitches in the autocorrelation results. Could it be that your calculations are done on 16 bits?


Beamforming techniques - pros and cons I would like some advice on which beamforming method to use from someone with experience in the field.

I have radar data from the ocean which has been transmitted by an omnidirectional transmitter and received by a linear array of 8 antennas with half wavelength spacing.

The signals come from all directions and I need to separate them using beamforming.

I started by using a simple phase shift beamformer which I wrote in Python but it wasn't filtering the signals from other directions well.

I've found the Capon beamforming technique to be an improvement. There's a toolbox on MATLAB which provides good functions for doing this.

However, my questions are:

What beamforming method would be best for this job?
What are the pros and cons of the different methods?

If any more details are needed, please ask.

For spatial filtering there are numerous beamforming techniques and it will be difficult to compare all here, there are numerous resource available on internet.

As related to your specified beam forming techniques,

Phased shift beamformer also known as conventional beamformer:
This beamforming techniques tries to add constructively the signal from the provide direction. This algorithm thus shifts the signal such that they align perfectly at the provided direction and thus sums up constructively. This algorithm has the worst angular resolution amongst different beam forming techniques as all signals from other directions are also included into the result. The resolution is in the range of the beam width of linear array.
Minimum variance distortionless response (MVDR) beamformer also known as CAPON beamformer: This beamformer is based on matching the array response matrix with the inverse of the sample correlation  matrix, thus enhancing the signal from considered direction and also suppressing the effects from the interfere signals from different direction other than the direction under consideration. The resolution of MVDR is far better than phased shift beamformer.

Among these two algorithm stated above the MVDR has benefit related to separation between two signals but on the expense of more computations.
It would be of worth to also look into the Linear constraint minimum variance beamformer (LCMV).

The implementation of MVDR and LCMV weights for beamformer is present in MATLAB MVDR weights and LCMV beamformer.


Eigenfunction property for LTI sinusoidal and the sinusoidal steady-state response All LTI systems possess the eigenfunction property for complex exponential inputs. That is (restricting our attention to periodic complex exponentials), if $e^{j\omega_k t}$ is an input to the LTI system, then the output is $H(j\omega_k) e^{j\omega_k t}$, where $H(j\omega_k)$ is the frequency response of the systsem evaluated at $\omega_k$. When a linear constant coefficient differential equation with zero initial conditions can represent a causal LTI system. Why is it that when we talk about systems, specifically systems with dynamics described by differential equations, we speak of a “steady-state” sinusoidal output? This seems to imply that transients are always present. But if the system is LTI to begin with, would their be transients in the output if a sinusoidal signal, or complex exponential, is applied? It can be shown mathematically that the output of an LTI system when a complex exponential is applied is just a scaled version of the input, with no transient terms. What’s going on here?

It can be shown mathematically that the output of an LTI system when a complex exponential is applied is just a scaled version of the input, with no transient terms.

that's for an input that has been going on since forever ago (and will continue forever to the future):

if
$$ x(t) = e^{j \omega_0 t} \qquad \forall t \in \mathbb{R} $$

then

$$ y(t) = H(j \omega_0) e^{j \omega_0 t} = H(j \omega_0) x(t) $$

  What’s going on here?

problem is, sometimes we have to apply an input at some time we might call "$t=0$", and before the LTI settles down to the $y(t)$ shown above, there are transients.  but if the LTI is stable (and not "marginally stable"), the transients will go away eventually.


Proof $GM$ (in dB) of $2nd$ order system is $\infty$ How can we prove that the Gain Margin $[GM]$ of $2nd$ order system is $\infty\quad ?$

My Approach:
Let us consider a $2nd$ order open loop system : $$G(s)= \frac{k}{(s+1)(s+2)}$$
Now, we know: $GM(in\quad dB)=-20log_{10} |G(jw_p)|$
$$ For \quad w_p \quad : \quad \angle GM(jw)=-180$$  $$\implies -tan^{-1}(w)\quad-tan^{-1}(w/2)=-180$$ $$\implies tan^{-1}(w)\quad+tan^{-1}(w/2)=180$$ $$\implies w=0 \quad or \quad w_p=0$$ $$\therefore |G(jw_p)|=k/(1 \times 2)=k/2$$ $$\implies GM(in\quad dB)=-20log_{10} (k/2) \neq \infty$$

so where is my mistake and how can i proof this results....??please help

For a second order transfer function

$$G(s)=\frac{k}{s^2+2\zeta \omega s +\omega^2}$$

What are the real part of the poles?

$$\Re\{s_{p1},s_{p2}\}=-\zeta\omega$$

The closed loop transfer function from a constant feedback $h$ is

$$G(s)=\frac{k}{s^2+2\zeta \omega s +\omega^2+kh}$$

$$s_{p1},s_{p2}=\frac{-2\zeta\omega\pm\sqrt{-4\omega^2(1-\zeta^2)-4kh}}{2}$$

$$s_{p1},s_{p2}=-\zeta\omega\pm\sqrt{-\omega^2(1-\zeta^2)-kh}=-\zeta\omega\pm\sqrt{\Delta}$$

What are the real part of the poles?

For a system with poles which have imaginary parts
$$\Re\{s_{p1},s_{p2}\}=-\zeta\omega$$

It means that it doesn't matter how big $k$ is.

But if the poles are pure real numbers, as $h>0$, by increase of the k, $\Delta$ decreases and only the upper bound of the poles reduce. So, the poles do not go to the right hand side of s-plane.


Single pole IIR filter, fixed point design We want to do a fixed point implementation of the single pole IIR filter:

$y[n] = a\cdot x[n] + (1 - a)\cdot y[n-1] \quad ;\qquad 0<a<1$

What are the main design considerations to think about ? In particular, is the fixed point design more challenging for smaller a or larger a ?

One thing to consider when implementing an IIR filter, whatever the order, is quantization and limit cycles.

Let me show you with a quick example with your original filter

$y[n] = a*x[n]+(1-a)*y[n-1]$

Let a = 0.005 and say that we use 16-bit signed coefficients.

$a_{fixedpoint} = a * 32768 = 164$

Let's assume that the input and output are 16-bit integers. They will go from -32768 to 32767. The previous output $y[n-1]$ is also a 16-bit integer in this scenario. The implementation would look like this, with temp being a 32-bit integer. All other variables are 16-bit integers.

$ temp = a_{fixedpoint}*x[n] + (32768 - a_{fixedpoint})*y[n-1]$

$y[n] = round(temp/32768)$

Let's try this implementation with a step response

It looks good, right? Perhaps we should zoom in.

We don't reach 32767 even though the gain of the original filter is 1. This is called a limit cycle. For order-1 IIR filters, limit cycles take the form of a DC offset that you can't get rid of. For IIR filters with complex poles, limit cycles take the form of a small oscillation that you can't get rid of. The problem is caused by the fact that $y[n-1]$ is stored on 16 bits. If we increase the resolution of $y[n-1]$ to 32 bits, it would solve the problem. Another solution would be to use fraction-saving.

If the input would go from 32767 to 0, we would have the same problem, the output would get close to 0 but would not actually reach 0.

Increased resolution

If instead of storing the output on 16 bits, we store it on 24 bits and name this variable acc (increased resolution of y[n]).

$ temp = a_{fixedpoint}*x[n]*256 + (32768 - a_{fixedpoint})*acc$

$ acc = temp/32768$

$y[n] = round(acc/(256))$

As we can see, with the increased resolution we can reach 32767 even though the coefficients are still on 16 bits. However the price to pay is having intermediate variables with greater resolution, the calculations might be slower.

fraction-saving

Finally there's a clever approach to minimize the added resolution needed for an IIR, especially when the poles are close to the unit circle. It is called fraction-saving. The idea is to memorize the round-off error of the previous sample and apply it to the next calculation to reduce the quantization effects. An additional variable is needed called error, but you don't need to increase the resolution of the previous output. So if x is represented on 16 bits, then y[n] can be stored on 16 bits too. The error variable would be represented with 16 bits too in this example. As Robert Bristow-Johnson pointed out, fraction-saving is a form of noise-shaping

$ temp = a_{fixedpoint}*x[n] + (32768 - a_{fixedpoint})*y[n-1] - error[n-1]$

$y[n] = round(temp/32768)$

$ error[n] = y[n]*32768 - temp$


Quadratic Programming with Linear Equality Constraints I need to solve an equality constrained minimization problem as give below
$$\min_{\textbf{w}} \mathbf{w}^TR\mathbf{w} $$
such that
$$X\mathbf{w} = \mathbf{1}$$
where $R\in \mathbb{R}^{n\times n}$ is covariance matrix (hence positive semi-definite)

$X\in \mathbb{R}^{m\times n}$ matrix with $m\gg n$,

$\mathbf{1}$ is column vector of size $m$ with all $1's$

$\mathbf{w}\in \mathbb{R}^n$ is an unknown vector.

Is there any closed form solution for this? If yes, can anyone provide it?

If no can we solve it by GD or any other algorithm?

Thanks in advance.

Let's solve a more general problem (Least Squares with Linear Equality Constraints):

$$
\begin{alignat*}{3}
\arg \min_{x} & \quad & \frac{1}{2} \left\| A x - b \right\|_{2}^{2} \\
\text{subject to} & \quad & C x = d
\end{alignat*}
$$

The Lagrangian is given by:

$$ L \left( x, \nu \right) = \frac{1}{2} \left\| A x - b \right\|_{2}^{2} + {\nu}^{T} \left( C x - d \right) $$

From KKT Conditions the optimal values of $ \hat{x}, \hat{\nu} $ obeys:

$$ \begin{bmatrix}
{A}^{T} A & {C}^{T} \\
C & 0
\end{bmatrix} \begin{bmatrix}
\hat{x} \\
\hat{\nu}
\end{bmatrix} = \begin{bmatrix}
{A}^{T} b \\
d
\end{bmatrix} $$

The trick here is to have a look at:

$$ \frac{1}{2} \left\| A x - b \right\|_{2}^{2} = \frac{1}{2} {x}^{T} {A}^{T} A x - {x}^{T} A b + \frac{1}{2} {b}^{T} b $$

So if we set $ w = x $, $ b = \boldsymbol{0} $, $ X = C $, $ d = \boldsymbol{1} $ and $ \frac{1}{2} {A}^{T} A = R $ then your solution is given by:

$$ \begin{bmatrix}
2 R & {X}^{T} \\
X & 0
\end{bmatrix} \begin{bmatrix}
\hat{w} \\
\hat{\nu}
\end{bmatrix} = \begin{bmatrix}
\boldsymbol{0} \\
\boldsymbol{1}
\end{bmatrix} $$

This is easy to solve in MATLAB or Python.

Handling Large System

As a request from @user5045 I add some info how to handle this in case $ R $ and $ X $ are large matrices.

We basically need to solve large scale matrix equation:

$$ \begin{bmatrix}
2 R & {X}^{T} \\
X & 0
\end{bmatrix} \begin{bmatrix}
\hat{w} \\
\hat{\nu}
\end{bmatrix} = \begin{bmatrix}
\boldsymbol{0} \\
\boldsymbol{1}
\end{bmatrix} = e = F g $$

The way to solve it is using an iterative solver.
I case $ R $ is a PSD matrix then the solver should be Preconditioned Conjugate Gradient. In MATLAB it is implemented using pcg().
In case $ R $ is only symmetric one should use Minimum Residual Solver. In MATLAB it is implemented as minres(). For a a general square matrix I'd go with cgs().


What Is the Meaning of the Term "Frequency Response'' When We Design a Pitch Detector This is from Klapuri 2003. What is the "frequency response" and why do we model is as multiplied with the vibrating system we want to get the frequency from?

Taking a periodic function with many harmonics, and then filtering it, is a common way of generating (or modeling) periodic signals.

It is how you can model human speech, where S is the sound produced by vocal cords, and H is the frequency response of the cavities formed by larynx, mouth, lungs, etc. And this is how many vocoders analyze and synthesize signals.

Imagine the signal S to be a periodic signal with all harmonics; ideally, S may be an impulse train, so that you have harmonics present at all integer multiples of the fundamental frequency.

Then, H is a filter which will affect each harmonic independently, in order to produce (or model) the desired periodic signal.

Since the variable $k$ corresponds to frequency, we can assume that they are actually DFT coefficients. Applying a filter H to a signal S is a multiplication in frequency (or a convolution in time).


How to see the amplitude of each frequency? do you know software that exports spectrograms to images of 500x22000 pixels?
all the spectrograms that I know show diffused frequencies and it does not look very good
I just want to see the amplitude of each frequency

each vertical line is one second of the sound
each horizontal line is a frequency

the entire audible spectrum enters the height of the image

the amplitude of each frequency is the brightness of the pixel

As I have not had enough reputation to comment, I will provide some of my interpretation of your question (Don't need to accept my answer).

It looks like that you want a spectrogram represented by a 500$\times$22000 matrix. This relates to the number of time and frequency slots. The number of time slots (columns) is determined by the sample number of the given signal. The number of frequency slots (rows) is determined by number of FFT output, which you can adjust by methods such as zero padding.

If you want to get better-looking spectrogram, you might need to adjust the window that is used in performing STFT (or maybe more overlapping in STFT). However, sometimes you just cannot get satisfying resolutions as the time-frequency resolutions are bound to the uncertainty principle. This is decided by your time-domain samples.


What is the meaning of $Ta_k$ of fourier series or transform? What is the meaning of $Ta_k$ of fourier series or transform?
I am taking a course on signal and systems.

In 286 page of my textbook, it says that as T becomes arbitrarily large the original periodic square wave approaches a rectangular pulse.
Also it says that all that remains in the time domain is an aperiodic signal corresponding to one period of the square wave.
(textbook: siganls and systems sencond edition, author: oppenheim)

I have a difficulty understanding this.. I can't connect this idea with fourier transform.

I suggest a link http://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-003-signals-and-systems-fall-2011/lecture-videos-and-slides/MIT6_003F11_lec16.pdf

The idea is that a Fourier series is only defined for periodic signals. In the discussion in the linked slides, the author is considering a rectangular pulse train with period $T$. That is, a pulse of width $2S$ repeats periodically with a spacing of $T$ between them. The pulses are therefore centered at:

$$
[\ldots, -2T, T, 0, T, 2T, \ldots ]
$$

Now, consider what happens as $T \to \infty$: in the limit, the only pulse that remains is the one centered at zero; the others are infinitely far away. When the author makes the claim that:

$$
\lim_{T\to \infty} T a_k = E(\omega)
$$

He or she is trying to show that the Fourier transform, which is defined for suitably well-formed aperiodic signals, can be thought of as the Fourier series of that signal (which typically wouldn't be defined since the signal is not periodic) in the limiting case of an infinite period. Stated a little differently, you can in some way think of an aperiodic signal as a periodic signal with infinite period.

The multiplication by $T$ in the limit is to account for the differences in definition between the Fourier series and Fourier transform: the series representation typically has a factor of $\frac{1}{T}$, while the transform does not. I don't know that there is a lot of insight to be gained via this analysis, but it shows that the series and transform representations are intimately related.


Frequency spectrum of complex equivalent Baseband Signal - Frequency Content around carrier frequency Basically my objective is to generate a complex equivalent baseband signal of an AM modulated wave in MATLAB. Then to view the difference in frequency spectrum of AM modulated RF signal vs complex equivalent baseband of (AM modulated RF signal).

I generated an equivalent complex baseband signal for an AM modulated signal and when viewing the spectra of complex baseband signal, I see frequency content centered at 0 but too many frequency content around carrier frequency as well.

Could anyone let me know why I am seeing frequency content around carrier frequency as well?

Here are the steps I followed:

Generated a carrier signal and message signal.
Did AM modulation.
Generated Inphase and Quadrature components.
Then represented complex baseband signal by $x_{bb}(t)=x_i(t)+jx_q(t)$
View the frequency spectra of AM modulated wave.
View the frequency spectra of Complex baseband signal.

When I compare the spectra, I see that the spectra of complex baseband signal is centered around 0 but has frequency components at the carrier frequency range (both positive and negative). Here are the screenshots:

Spectra - AM modulated signal

Spectra - Equivalent complex baseband signal

Here is my MATLAB code:

fs=15000;
t=0:1/fs:1;
fm=5;
m=sin(2*pi*t*fm);  %message signal
%plot(t,m,'b');
fx=1000;
x=cos(2*pi*t*fx); %carrier signal
xm=m.*x;  %AM
ax=subplot(2,1,1);
plot(ax,t,m,'b');
axm=subplot(2,1,2);
plot(axm,t,xm,'g');
hx=imag(hilbert(xm)); %hilbert transform

xit=(xm.*cos(2*pi*fx*t))+(hx.*cos(2*pi*fx*t));
xqt=(hx.*cos(2*pi*fx*t))-(xm.*sin(2*pi*fx*t));

xbb=xit+(1i*xqt); %equivalent complex baseband signal
figure(2)
plot(t,xbb);

%spectrum
nfft=1024;
F=fftshift(fft(xm,nfft));      %spectra of AM modulated wave
fvals=fs*(-nfft/2:nfft/2-1)/nfft;
figure(4);
stem(fvals,abs(F),'b');

%spectrum
nfft=1024;
F=fftshift(fft(xbb,nfft));  %spectra of equivalent complex baseband
fvals=fs*(-nfft/2:nfft/2-1)/nfft;
figure(5);
stem(fvals,abs(F),'b');

You have a bug in the formulation of your complex heterodyne oscillator signal inside your complex multiplication.  This seems to be creating both sum and difference frequency results.  Thus the extra high-frequency spectrum.


What is the difference between a simple gaussian filter and gaussian filter multiplied by its sum of it elements? 1 2 1                         1 2 1
2 4 2 (A)            (1/16) * 2 4 2  (B)
1 2 1                         1 2 1

Both matrices are the same except that the matrix B is multiplied by the sum of its elements.
When convolved with an image, they give nearly the same blurred image.

What is the exact difference between them and when they should be used ?

The rightmost one (where you divide by the sum) ensures that the output of the filter wil have the same dynamic range as the input.

Actually, the output blurred images are not the same because their pixel values are different. You have the feeling that they are similar because your software (I assume Matlab) does remap the intensity values of the pixels to a range that is acceptable by typical displays, i.e., 8-bits values between 0 and 255. If you used a more "close to the metal" software or library, then you would have to do this operation by yourself.

Data dynamic range

The dynamic range of an image is given either by the min and max intensity values present in the image or by the min and max intensity that can be represented by the digital format used to store the intensity values. (From the context it is usually clear which one of the definitions is in use.)

Typical ranges are $[2,255]$ (with integer spacing) that corresponds to 8 bits quantization and $[0,1]$ for floating point values (either single or double precision). The same convention is used in Computer Graphics, OpenGL, etc.

If you take a constant image with value 1 everywhere:

the left filter will output a constant image with value 16 everywhere,
the second (rightmost) one produces a constant image of value 1.

From this simple example, you can see that one filter preserved the enery of the image, while the second one increased it.

While it is harmless in this case, it can lead to numerous problems in practice.
For example, iterative algorithms may diverge (because of numerical overflows), your computer may run out of bits to represent your data, or you may lose precision (single and double precision floating-point values do not have a uniform precision over the range of the possible values).


Waveform samples and I/Q components I'm quite new to DSP and I guess my problem it's quite silly but I can't find a clear answer. Basically I'm writing some python code for a test. The test requires me to write a memory with the I/Q components of a specific waveform.

As far as I know it's quite easy in hardware to get the I/Q samples. I just need to multiply the input signal ( my waveform ) with $sin(w_ct)$ and $cos(w_ct)$ ($w_c$ should be the sampling frequency) and then sample it with an ADC.

In my case I'm able to create the samples of the waveform ( it's a sawtooth but I would like to understand the general case as well ) but then I'm quite stuck because I don't know how to extract the I/Q components from each sample.

Any idea ?

Cheers,
Stefano.

Assume that signal that comes to you through your antenna is:

$$ s_r[n] = I[n] \cos(\omega_0 n) + Q[n] \sin( \omega_0 n ) $$

which is a bandpass, modulated, narrow band signal where the components $I[n]$ and $Q[n]$ are the lowpass baseband in-phase and quadrature-phase components respectively.

Then in order to obtain the I and Q components from $s_r[n]$ you can apply the following just like the analog case:

$$I[n] = h_{LP}[n] \star ( s_r[n] \cos(\omega_0 n) ) $$
and
$$Q[n] = h_{LP}[n] \star ( s_r[n] \sin(\omega_0 n) ) $$

Where the lowpass filter $h_{LP}[n]$ is will have a cutoff frequency larger than the bandwidth of the signals I and Q and less than any adjacent band (if there is any). You may also apply a scaling such as $2x$ to adjust the magnitudes.

If you would like to form the analytic lowpass signal from the I and Q then you can easily form it like:

$$x_a[n] = I[n] + j Q[n]$$

where it can be shown that the analytic signal is also given by

$$ x_a[n] = A[n] e^{ j\phi(n) } $$

where $A[n] = \sqrt{ I[n]^2 + Q[n]^2 } $ and $\phi(n) = \tan^{-1}(\frac{Q[n]}{I[n]})$

Depending on the application, information in $A[n]$ (the envelope) or in $\phi(n)$ (the phase) is used.


Difficulties in understanding mutual information concept What does mutual information (MI) convey? If 2 signals are independent then MI is zero; What does this imply and mean in the case of mutual information of error, and in general definition of MI and entropy.
What is the interpretation if MI of error is decreasing or getting minimized?

What does mutual information (MI) convey?

It indicates that there is a relationship between the two signals- i.e. that they are not independent.  It could be that they are correlated, but the relationship does not have to be that simple.  The higher the mutual information, the stronger the relationship (i.e. the more one of the signals can "tell" you about the other).

  What does this imply and mean in the case of mutual information of error, and in general definition of MI and entropy.

If the error has mutual information with the observation this implies that the receiver did not use all of the information in the observation- i.e. the receiver is not ideal.  If the receiver was ideal it would have used all of the information in the observation, so what was left in the error signal would be independent of the observation.

  What is the interpretation if MI of error is decreasing or getting minimized?

Decreasing MI of error implies that the receiver is doing a better job.  Once the MI is zero it is impossible to do any better.


How to Classify a Kernel as Low Pass Filter (LPF) or High Pass Filter (HPF)? How to Transform an LPF Kernel into HPF Kernel? The filter impulse response is [1/9 1/9 1/9;1/9 1/9 1/9;1/9 1/9 1/9].
how do i determine whether it is an Low pass filter or a high pass filter. If this is a low pass filter how do i convert it to a high pass filter or the opposite.

General

We assume 2 modes of filters: LPF or HPF.

Classifying Filter Type

Usually, if it is a well planned LPF and well Planned HPF a simple test will do.
Calculate the sum of all coefficients.
The sum of the coefficients is the first element of the DFT of the signal.
It means it is the DC gain and well behaved LPF has gain of 1 and well behaved HPF have DC gain of 0.

So sum the coefficients and:

If the sum is $ 1 $ then this is an LPF.
If the sum is $ 0 $ then this is an HPF.

Pay attention that those method can't differentiate between HPF and Band Pass or between LPF and Band Stop.

The more general method is to look on its DFT Spectrum to see its behavior.

Converting Filter Type

There is very intuitive formula - ALL_PASS_FILTER = LPF_FILTER + HPF_FILTER.
It's not accurate but this is the idea.

Be more accurate in the Frequency Domain All Pass Filter is a filter which has value of 1 for any frequency.
Assuming we have LPF Filter then HPF_FILTER = ALL_PASS_FILTER - LPF_FILTER.
In the spatial domain it means $ h \left[ n \right] = \delta \left[ n \right] - l \left[ n \right] $.

Where $ \delta \left[ n \right] $ is the filter with 1 in its origin and $ 0 $ anywhere else.

In your case $ \delta \left[ n \right] = [0, 0, 0, 0, 1, 0, 0, 0, 0] $.
So the HPF Filter becomes $ h \left[ n \right] = \delta \left[ n \right] - l \left[ n \right] = [-1/9, -1/9, -1/9, -1/9, 8/9, -1/9, -1/9, -1/9, -1/9] $.


Sampled and aliasing signal () = cos(200) + 2cos(320). (), is produced by sampling () at sampling frequency $f = 300 Hz$.
Is the signal sampled or aliased and why?

Another perspective: The sampling process can be viewed as a type of amplitude modulation. PCM stands for Pulse Code Modulation. Let's set aside Code as simply the encoding of an amplitude as a digital value, leaving Pulse Modulation. This is the process of multiplying, or amplitude modulation of, the analog signal with a pulse train.

The results of amplitude modulation is well known—it produces the sum and difference frequencies of the component sinusoids of each. Amplitude modulation of 100 Hz by 9 Hz produces frequencies of 109 Hz and 91 Hz. For signals with more sinusoidal components, the result is the sum of all combinations of AM of the individual components.

So, if we know the harmonic content of a pulse train, we know the results of modulating another signal with it. The harmonic series of a pulse train in time is a pulse train in frequency. That is, a pulse train at frequency Fs (our sampling frequency) has harmonics that are cosines of equal amplitude at integer multiples of Fs.

And fortunately, the problem is made simple since we are primarily concerned with the harmonics of the pulse train that frame the band we are sampling (below half the sample rate). The first non-negative harmonics are 0 Fs and 1 Fs. For your sample rate of 300 Hz, amplitude modulation of your signal by 0 Hz and 300 Hz.

Modulation by 0 Hz (again, a cosine, of non-zero amplitude, so it's simply a constant offset and often referred to as DC) simply gives the spectrum of your signal being sampled.

Your signal has harmonics of 100 Hz and 160 Hz, so that will also be the spectral components we get from AM by DC.

Amplitude modulation by 300 Hz yields a spectrum that is the sum and difference. We can ignore the sum, since we want to know if there is aliasing in our band of interest. For modulation by 300 Hz, the difference frequencies resulting from AM with your signal will be (300 - 100) Hz, or 200 Hz, and (300 - 160) Hz, or 140 Hz. Here is the spectrum with that contrubution, in red:

Accounting for these, we have a spectrum with 100 Hz, 140 Hz, 160 Hz, and 200 Hz. When we satisfy reconstruction by lowpass filtering to pass everything below half the sample rate, or 150 Hz, we end up reproducing 100 Hz and 140 Hz, instead of the original 100 Hz and 160 Hz. Here is the spectrum of the result, reconstructed for the analog domain:

Despite the long explanation, once you know this, the visualization is easy—simply mirror the spectrum of the original signal off the sampling frequency.


3D FFT's and Convolutions using FFTW I am trying to perform a of convolution 2 3-dimensional vectors using the FFTW3 library in C.

My question relates to the order in which the frequency domain appears and if it will make a difference to the final result.

In the FFTW documentation it dictates that in 1D complex FFT's the negative frequencies are in the second half of the array. Does this follow through with a 3D complex FFT as the documentation doesn't say? If this is the case how would it be stored? i.e in 1D it would be something like:

|+|+|+|+|-|-|-|-|

so what would it be in 3D (sorry, terrible as visualising these things).

Finally when convolving 2 3D vectors using the convolution theorem do I need to do anything because of the FFTW layout, or can I just do the pointwise product and then IFFT back?

In 1D, you are likely to obtain, more precisely, $f_0, f_1, f_2,\ldots, f_{N/2}, f_{-N/2+1},\ldots, f_{-1}$. In other words, the $N/2+1$ ones correspond to positive frequencies. In higher dimensions, the same applies. In 3D, you ought to get a "top-left-front" sub-cube of positive frequencies.

To reduce periodic artifacts, you could 3D window the data before, and take care of the values of the imaginary parts after IFFTing back. If the input data is real, they should be small, and you can keep the real part only. If not, something gone wrong.


How to find transfer function by state space representation matrices A state space representation is given by:
$$\dot{x}= \begin{bmatrix}0 & 0 & 0 & 0\\ 1 & 0 & 0 & 0 \\ 0&0&-2&-4\\0&0&1&0\end{bmatrix}x+\begin{bmatrix} 1\\0\\1\\0\end{bmatrix} u$$
$$y=\begin{bmatrix}0 & 1 & 0 & -1\end{bmatrix}x$$
I need to find a transfer function of this system. I don't know how, but I've tried it this way:

I wrote all equations:
$$\dot x_1 = u$$
$$\dot x_2 = x_1$$
$$\dot x_3 = -x_3-4x_4+u$$
$$\dot x_4=x_3$$
$$y=x_2-x_4$$

And I've draw a scheme of such system:

And then tried to find its transfer function.
$$-4X_4-2sX_4+U = X_4s^2$$
$$U = X_4(s^2+4+2s)$$
$$X_4= \frac{U}{s^2+2s+4} $$
$$X_2 = \frac{U}{s^2}$$

$$Y = X_2 - X_4 = \frac{U}{s^2} -\frac{U}{s^2+2s+4}$$
$$G(s) = \frac{Y(s)}{U(s)} = \frac{1}{s^2} -\frac{1}{s^2+2s+4}=\frac{2(s+2)}{s^2(s^2+2s+4)}$$
I don't know, whether this is right answer nor solution.

What solution would you recomend for this task? Is the solution above correct?

Use the fact that $\dot x_1 = s x_1$ and take it from there. It appears that your result is correct.

It should be noted that there is a closed form expression to map a state-space representation to a transfer function. However, it includes inverting a matrix with s variables and may hence be more challenging than just solving it as stated above


Derivation of Inverse Fourier transform from forward Fourier transform Consider the Fourier pairs:
$$\psi(x,t) \stackrel{\mathrm{FT}}{\longleftrightarrow} \Psi(k,t)$$
$$\text{If } \quad \quad\Psi(k,t)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} \psi(x,t) e^{-ikx} \, dx \quad \quad \dots(i)$$
$$\text{then, can we derive: } \quad \quad \psi(x,t)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} \Psi(k,t) e^{ikx} \, dk \quad \quad \dots(ii) \quad ?$$
I have used comparison method to proof $eq.(ii) $
i.e., $$\text{consider:} \quad \quad x(t) \stackrel{\mathrm{FT}}{\longleftrightarrow} X(\omega)$$
$$\text{then, } \quad \quad X(\omega)=\int_{-\infty}^{\infty} x(t) e^{-i\omega t} \, dt \quad \quad \dots(iii)$$
$$\text{and, } \quad \quad x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty} X(\omega) e^{i\omega t} \, d\omega \quad \quad \dots(iv)$$
Now, comparing $eq(iii)$ with $eq(i)$ , we find:
$w \to k$
$t \to x$
$x(t) \to \frac{\psi(x,t)}{\sqrt{2\pi}}$
$X(w) \to \Psi(k,t)$ , putting these values in $eq(iv)$, we get:
$$\frac{\psi(x,t)}{\sqrt{2\pi}}=\frac{1}{2\pi}\int_{-\infty}^{\infty} \Psi(k,t) e^{ikx} \, dk $$
$$\implies \psi(x,t)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} \Psi(k,t) e^{ikx} \, dk$$
But can we derive $eq(ii)$ from $eq(i)$? {without using any comparison}

I'll use the common DSP definition and notation of the Fourier transform and its inverse like given in Eqs $(iii)$ and $(iv)$ in the OP. Note that the argument $t$ in Eqs $(i)$ and $(ii)$ is irrelevant.
First of all, note that the Fourier transform of a delay $x(t)=\delta(t-\tau)$ is given by
$$\int_{-\infty}^{\infty}\delta(t-\tau)e^{-j\omega t}dt=e^{-j\omega \tau}\tag{1}$$
Consequently, we have to accept the following expression (in a distributional sense), following from the inverse Fourier transform of $(1)$:
$$\delta(t-\tau)=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{j\omega (t-\tau)}d\omega\tag{2}$$
Now we can prove $(iv)$ from $(iii)$:
$$\begin{align}\int_{-\infty}^{\infty}X(j\omega)e^{j\omega t}d\omega&=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}x(\tau)e^{-j\omega\tau}d\tau e^{j\omega t}d\omega\\&=\int_{-\infty}^{\infty}x(\tau)\underbrace{\int_{-\infty}^{\infty}e^{j\omega (t-\tau)}d\omega}_{2\pi\delta(t-\tau)} d\tau\\&=2\pi x(t)\end{align}\tag{3}$$


Errors when converting from digital filter object to transfer function I am trying to use the designfilt tools in MATLAB to create a bandpass filter, but to be compatible with other code I need the resultant filter in transfer function coefficient form with B and A. To start with I took an example from the filter design gallery, and just changed the sampling rate to match that of my real data.

Fstop1 = 150;
Fpass1 = 200;
Fpass2 = 300;
Fstop2 = 350;
Astop1 = 65;
Apass  = 0.5;
Astop2 = 65;
Fs=16384;
% Fs=1000;

d = designfilt('bandpassiir', ...
  'StopbandFrequency1',Fstop1,'PassbandFrequency1', Fpass1, ...
  'PassbandFrequency2',Fpass2,'StopbandFrequency2', Fstop2, ...
  'StopbandAttenuation1',Astop1,'PassbandRipple', Apass, ...
  'StopbandAttenuation2',Astop2, ...
  'DesignMethod','butter','SampleRate', Fs);

This filter seems to be behaving itself, isstable returns 1 and the $z$-plane seems acceptable:

However, when I try to get the coefficients using [B,A]=d.tf; the resultant filter is no longer stable and the $z$-plane looks totally different:

The values for B are incredibly small around e-26, so I am imaging that it might be to do with gains of the coefficients? If I set Fs back to the original $1\textrm{kHz}$, this two $z$-planes plots are the same.

Can someone help me with sorting this problem? Is the fact that the transfer function form is unstable speak to the stability of the filter when it is in the digitalfilter object form?

edit

The documentation for the butter function even describes these errors here

Generally you should never use the Transfer Function syntax. Look how your poles are closely placed next to the unit circle. With high order Butterworth filters you will see the effect of round-off numerical errors. What is the high order? Well 5 is already high and your filter has 30.

Instead you should use SOS or Zero-Pole-Gain representation. For example try to visualise the same filter but with SOS matrix:

fvtool(d.Coefficients) % or: zplane(sos2zp(d.Coefficients))

You will see that there are no visible numerical errors.


Does Wavelet Packets (WPT) decompose signal on dyadic scales like DWT? As far a I understand WPT is a generalization of DWT in a sense that you get a whole binary tree of possible orthogonal decompositions where DWT is only a one branch. And WPT also downsamples a signal at every level by the factor of 2.

Am I correct?

But at the same time there are statements that WPT gives better precision for high frequencies compared to DWT.

Why is it so? Is it because it contains possible orthogonal decompositions where we decompose details of the previous level?

It seems that two-scale equations for WPT imply dyadic scales ($\frac{t}{2}$ on the lhs)

$$2^{-\frac{1}{2}}W_{2m}(\frac{t}{2} - k) = \sum^{\infty}_{l=-\infty}h_{l-2k}W_m(t-k)$$
$$2^{-\frac{1}{2}}W_{2m+1}(\frac{t}{2} - k) = \sum^{\infty}_{l=-\infty}g_{l-2k}W_m(t-k)$$

So, is it correct to assume that WPT like DWT works on dyadic scales?

I suppose  you are talking about $2$-channel DWT. WPT may have  slightly different  interpretations. With DWT, one only allows to decompose frequencies "downward ", ie after lowpass and subsampling. WPT allows to further decompose, at each level, after (lowpass or highpass) and subsampling. In that sense, WPT are a generic scheme, of which the DWT denotes a specific branch. The decision to decompose or stop at each level is governed by rules "external" to wavelets: frequency target, entropy, etc.

Another vision is the uniform WPT, or full wavelet packet, where at each level, both outputs are further decomposed, generating a uniform frequency splitting.

So basically, you are correct at question 1. For question 2, when a DWT subband is approximately located (in reduced frequencies) between two successive dyadic bounds $[2^{-j-1},2^{-j}]$, a WPT can be located between any two dyadic rational bounds $[k_1 2^{-j_1},k_2 2^{-j_2}]$, with $j, j_1,j_2$ nonnegative integers, $k_1 \le 2^{j_1}$ and $k_2 \le 2^{j_2}$.

For instance, the upper DWT subband only (approximately) localizes stuff in $[1/2,1]$, while you can reach  $[5/8,3/4]$ with a WPT. However, the aliasing comes into play, and this $[5/8,3/4]$ is disturbed by aliased components from other frequencies.  Aliasing is considered more severe for high frequency wavelet packet coefficients.

With $M$-band DWT and WPT, both structures remains $M$-adic at each level, but the WPT tends to split the subband axis onto $M$-adic fractions, with form $k/M^j$.


Calculating smoothed derivative of a signal by using difference with larger step=convolving with rectangular window I have a signal sampled at $\Delta t: fi(ti=i\Delta t)$ where i = 0..n-1. I want to find the first derivative of the signal: f'(t).

My first thought was to estimate this by a central difference:

$f'(t_i) =\frac{f(t_{i+1})−f(t_{i−1})}{2\Delta t}$

However the signal may have a lot of high frequency noise that may cause quick fluctuations in f'.
I guess the proper thing might be to smooth the signal by convolving with a window function e.g., Hann and then find the derivative from the differences.

A colleague suggested a faster way of finding a smoothed estimate of the derivative:
use a central difference over 2n samples, where n >> 1:

$f'(t_i) =\frac{f(t_{i+n})−f(t_{i−n})}{2n\Delta t}$

This would of course be computationally faster than first convolving with a window function but is it a good solution?

If we form the sum:

$S=2\Delta t[f'(t_{i-n+1})+f'(t_{i-n+2})+..+f'(t_{i+n-1})]$

and expand each derivative by the central difference with step $\Delta t$:

$S=f(t_{i-n+2})-f(t_{i-n})+f(t_{i-n+3})-f(t_{i-n+2})+..+f(t_{i+n})-f(t_{i+n-2})$

all terms except two cancels out:

$S=f(t_{i+n})-f(t_{i-n})=2n\Delta tf'(t_i)$

Therefore:

$f'(t_i)=\frac{1}{n}[f'(t_{i-n+1})+f'(t_{i-n+2})+..+f'(t_{i+n-1})]$

So taking the central difference over 2n samples is equivalent with first convolving by a rectangular window of size 2n - 2 and then taking a central difference over +/- 1 sample.

How "bad" is it to smooth with a rectangular window?

If we take the FFT this will cause "ringing", but we do not need to take the FFT.

Thanks in advance for any answers!

This is a hard question to handle generally. Smoothing with a rectangular window is used all the time (often called a "moving average"), so that's not necessarily a problem. I'm not sure what ringing you're referring to, perhaps the sidelobes of the rectangular window's frequency response.

Differentiation is inherently a highpass operation; the ideal continuous-time differentiator has a transfer function of:

$$
H(s) = s
$$

So that its magnitude response is:

$$
|H(j\omega)| = \omega
$$

The gain of a differentiator therefore increases monotonically with frequency. If your signal contains high-frequency noise, then it could be amplified by applying a differentiator. To combat this, two approaches are apparent:

Design a more sophisticated differentiator filter that has the desired linear magnitude response over the portion of the band that covers your signal of interest, then sharply attenuates higher frequencies. You could design such a filter using least-squares methods or the frequency-sampling method, for example.
Use a cascaded approach where you first suppress all of the high-frequency noise that you can using a lowpass filter, then follow with the differentiator. The frequency coverage of the differentiator need not be as tight since the lowpass filter will eliminate the out-of-band noise.

The methods should be roughly equivalent if you're using linear filters; you could think of the first single-filter approach as just a cascade of a differentiator and a lowpass filter. As you noted, the central-difference approach can be modeled in this way. It's hard for anyone to say that it's "bad" without any knowledge of your application. My main thought would be that it is "bad" if the smoothing operation tangibly attenuates your signal of interest, such that the derivative estimate is no longer useful. However, if the parameters of the signal are such that you can smooth the noise without noticeably distorting the signal (i.e. if the signal is well-oversampled), then it could be a win.


Finding the auto-correlation sequence $r_{xx}[k]$ for an AR(2) process Consider the following recursive difference equation of a LTI system, where $v[n]$ is a white noise, zero-mean process with $\sigma_v^2 = 1$.

$
x[n] = v[n] + 0.75x[n-1]-0.25x[n-2]
$

I want to calculate the auto-correlation sequence $r_{xx}[k]$ for $k=0,1,2$. I started as follows:

$
r_{xx}[k] = E\{(v[n] + 0.75x[n-1]-0.25x[n-2])(v[n-k] + 0.75x[n-k-1]-0.25x[n-k-2])\} \\
r_{xx}[0] = E\{ v[n]^2 \} + 0.75^2 E\{ x[n-1]^2 \} - 2 \times 0.25 \times 0.75 E\{ x[n-1] x[n-2] \} + 0.25^2 E\{ x[n-2]^2 \} = \sigma_v^2 + \frac{5}{8} r_{xx}[0] - \frac{3}{8}r_{xx}[1]
$

As you can see, in the equation for $r_{xx}[0]$, a term with $r_{xx}[1]$ remains. If I then want to determine $r_{xx}[1]$, a term with $r_{xx}[2]$ remains, and so on. How is it possible to solve this problem?

We have
\begin{equation}
x[n] = v[n] + 0.75x[n-1]-0.25x[n-2]
\end{equation}
You started correct, you get
\begin{equation}
 \begin{split}
  r_{xx}[k] = E\{(v[n] + 0.75x[n-1]-0.25x[n-2])(v[n-k] + 0.75x[n-k-1]-0.25x[n-k-2])\} \\
 \end{split}
\end{equation}
You get 9 terms
\begin{equation}
 \begin{split}
  r_{xx}[k] &= Ev[n]v[n-k] +0.75E x[n-k-1]v[n] - 0.25E x[n-k-2]v[n]\\
   &+0.75 Ex[n-1]v[n-k] + 0.75^2Ex[n-1]x[n-k-1]+0.75(0.25)Ex[n-1]x[n-k-2]\\
   &-0.25Ex[n-2]v[n-k] -0.75(0.25)Ex[n-2]x[n-k-1] + 0.25^2Ex[n-2]x[n-k-2]
 \end{split}
 \tag{1}
\end{equation}
Since the noise is white then
\begin{equation}
 Ev[n]v[n-k]  = \delta(k)\sigma_v^2 = \delta(k)
\end{equation}
Since $k\geq 0$, write down $x[n-k-1]$ and $x[n-k-2]$ do not contain lags that align with $v[n]$, so
\begin{align}
 E x[n-k-1]v[n] &= 0\\
 E x[n-k-2]v[n] &= 0
\end{align}
On the other hand, $Ex[n-1]v[n-k]$ and $Ex[n-2]v[n-k]$, will agree for $k = 1$ and $k=2$ according to
\begin{align}
x[n-1]v[n-k] &= \overbrace{Ev[n-1]v[n-k]}^{\delta(k-1)} + 0.75\overbrace{ Ex[n-2]v[n-k] }^{\delta(k-2)}- 0.25\overbrace{Ex[n-3]v[n-k]}^{0} \\
x[n-2]v[n-k] &= \underbrace{Ev[n-2]v[n-k]}_{\delta(k-2)} + \underbrace{E(0.75x[n-3]-0.25x[n-4])v[n-k]}_{0}
\end{align}
So
\begin{align}
 x[n-1]v[n-k] &= \delta(k-1) + 0.75\delta(k-2) \\
 x[n-2]v[n-k] &= \delta(k-2)
\end{align}
The other four terms follow from the AR definition, i.e.
\begin{align}
 E x[n-1]x[n-k-1] &= r_{xx}[n-1-n+k+1]= r[k] \\
 E x[n-1]x[n-k-2] &= r_{xx}[n-1-n+k+2]= r[k+1] \\
 E x[n-2]x[n-k-1] &= r_{xx}[n-2-n+k+1]= r[k-1] \\
 E x[n-2]x[n-k-2] &= r_{xx}[n-2-n+k+2]= r[k]
\end{align}
Replacing all results in equation $(1)$, we get
\begin{equation}
\begin{split}
  r[k] &= \delta(k) +0.75\delta(k-1) + 0.75^2\delta(k-1) -0.25\delta(k-2)\\
   &+0.75^2 r[k] +0.75(0.25)r[k+1]-0.75(0.25)r[k-1] + 0.25^2 r[k]
\end{split}
\end{equation}
which suggests that
\begin{equation}
 0.375r[k] = \delta(k) + 1.3125\delta(k-1) - 0.25\delta(k-2) +0.1875(r[k+1] - r[k-1])
\end{equation}
For $k = 0$, we get
\begin{equation}
 0.375r[0] = \delta(0) + 1.3125\delta(-1) - 0.25\delta(-2) +0.1875(r[+1] - r[-1])
\end{equation}
we know that $\delta(0) = 1$ and $\delta(-1) = \delta(-2) = 0$. Also $r[1] = r[-1]$ due to the fact that $r[k]$ is even
\begin{equation}
 r[0] =\frac{1}{0.375} \delta(0) = \frac{8}{3}
\end{equation}
For $k = 1$, we get
\begin{equation}
 0.375r[1] = 1.3125 + 0.1875(r[2] - r[0])
\end{equation}
which is (again because $r[k]$ is an even function)
\begin{equation}
 0.375r[1] = 1.3125 + 0.1875(r[2] - \frac{8}{3}) \tag{*}
\end{equation}
and at $k = -1$
\begin{equation}
 0.375r[-1] = 0.1875(r[0] - r[-2])
\end{equation}
which is (again because $r[k]$ is an even function)
\begin{equation}
 0.375r[1] = 0.1875(\frac{8}{3} - r[2])\tag{**}
\end{equation}
Now solve the system of two equations in two unknowns $(*),(**)$ to get $r[1],r[2]$, i.e.
\begin{equation}
 \begin{bmatrix}
    0.3750 &  -0.1875\\
    0.3750  &  0.1875
 \end{bmatrix}
 \begin{bmatrix}
  r[1]\\
  r[2]
 \end{bmatrix}
 =
 \begin{bmatrix}
    0.8125\\
    0.5000
 \end{bmatrix}
\end{equation}


Is there a clever way to implement cascaded moving average filters? I'd like to effeciently approximate a gaussian filter's step response using cascaded moving average filters.

I know about recursive moving average, but is there some clever algorithm to cascade them apart from just putting them in series? CIC filters seem related but are mostly used for decimation, I'm interested in pulse shaping.

[EDIT:20181225, added a couple of references] This topic is coming back. P. Getreuer provides C code, and lists an handful of them in A survey of Gaussian convolution algorithms,  Image Processing On Line, 2013:

  Gaussian convolution is a common operation and building block for
  algorithms in signal and image processing. Consequently, its efficient
  computation is important, and many fast approximations have been
  proposed. In this survey, we discuss approximate Gaussian convolution
  based on finite impulse response filters, DFT and DCT based
  convolution, box filters, and several recursive filters. Since
  boundary handling is sometimes overlooked in the original works, we
  pay particular attention to develop it here. We perform numerical
  experiments to compare the speed and quality of the algorithms.

Here are a couple references:

Robinson, J. A., Efficient Gaussian filtering using Cascaded Prefix Sums, 2012

  This paper introduces a method for multidimensional Gaussian filtering
  using an efficient one-pass cascade of overlapping local average
  windows driven by prefix sums. Each local-average filter is
  implemented in n dimensions, with non-integer lengths, allowing
  accurate approximation of Gaussians of any variance. In axis oriented
  form the method has a scan-rate hardware realization and fast software
  implementation using minimal extra memory. In this latter case the new
  method consistently outperforms the fastest alternative Gaussian
  filtering method both in accuracy and speed.

Bernardo Rodrigues Pires, Karanhaar Singh, José M. F. Moura, Approximating image filters with box filters, 2011

  Box filters have been used to speed up many computation-intensive
  operations in Image Processing and Computer Vision. They have the
  advantage of being fast to compute, but their adoption has been
  hampered by the fact that they present serious restrictions to filter
  construction. This paper relaxes these restrictions by presenting a
  method for automatically approximating an arbitrary 2-D filter by a
  box filter. To develop our method, we first formulate the
  approximation as a minimization problem and show that it is possible
  to find a closed form solution to a subset of the parameters of the
  box filter. To solve for the remaining parameters of the
  approximation, we develop two algorithms: Exhaustive Search for small
  filters and Directed Search for large filters. Experimental results
  show the validity of the proposed method.

Elboher, Elhanan  and Werman, Michael, Efficient and Accurate Gaussian Image Filtering Using Running Sums, 2011

  This paper presents a simple and efficient method to convolve an image
  with a Gaussian kernel. The computation is performed in a constant
  number of operations per pixel using running sums along the image rows
  and columns. We investigate the error function used for kernel
  approximation and its relation to the properties of the input signal.
  Based on natural image statistics we propose a quadratic form kernel
  error function so that the output image l2 error is minimized. We
  apply the proposed approach to approximate the Gaussian kernel by
  linear combination of constant functions. This results in very
  efficient Gaussian filtering method. Our experiments show that the
  proposed technique is faster than state of the art methods while
  preserving a similar accuracy.

Kovesi, Peter, Fast Almost-Gaussian Filtering, 2010

  Image averaging can be performed very efficiently using either
  separable moving average filters or by using summed area tables, also
  known as integral images. Both these methods allow averaging to be
  performed at a small fixed cost per pixel, independent of the
  averaging filter size. Repeated filtering with averaging filters can
  be used to approximate Gaussian filtering. Thus a good approximation
  to Gaussian filtering can be achieved at a fixed cost per pixel
  independent of filter size. This paper describes how to determine the
  averaging filters that one needs to approximate a Gaussian with a
  specified standard deviation. The design of bandpass filters from the
  difference of Gaussians is also analysed. It is shown that difference
  of Gaussian bandpass filters share some of the attributes of log-Gabor
  filters in that they have a relatively symmetric transfer function
  when viewed on a logarithmic frequency scale and can be constructed
  with large bandwidths.

Wells, William M., Efficient Synthesis of Gaussian Filters by Cascaded Uniform Filters, 1986

  Gaussian filtering is an important tool in image processing and
  computer vision. In this paper we discuss the background of Gaussian
  filtering and look at some methods for implementing it. Consideration
  of the central limit theorem suggests using a cascade of simple''
  filters as a means of computing Gaussian filters. Amongsimple''
  filters, uniform-coefficient finite-impulse-response digital filters
  are especially economical to implement. The idea of cascaded uniform
  filters has been around for a while [13], [16]. We show that this
  method is economical to implement, has good filtering characteristics,
  and is appropriate for hardware implementation. We point out an
  equivalence to one of Burt's methods 1, 3 under certain
  circumstances. As an extension, we describe an approach to
  implementing a Gaussian Pyramid which requires approximately two
  addition operations per pixel, per level, per dimension. We examine
  tradeoffs in choosing an algorithm for Gaussian filtering, and finally
  discuss an implementation.


Filter message from noisy voice signal I am trying to decipher a message hidden within a very noisy audio file(.wav)(I think it is white noise with an additional low drone).  The message is a six digit number. I have not have any further details about the noise.

I attempted to use a low-pass filter in the hopes that eliminating most of the higher frequencies would allow me to hear the numbers but, I seem to be unable to also get rid of enough of the low drone to hear the voice well enough. My attempt was as follows (the employed function freq_space_low_pass_filter is included at the end):

[data, SampleRate, NbitsPerSample]=wavread('noisy_msg6.wav');

y=data(:,1); % we will work only with one channel in this demo
N=length(y); %number of sample points
t=( (1:N)*1/SampleRate ).'; % time spacing is 1/SampleRate and we want column vector

Y=fft(y);

spectrum_freq=fourier_frequencies(SampleRate, N);

Freq3db=100;
[spectrum_filtered,g_vs_freq]=freq_space_low_pass_filter(Y, SampleRate, Freq3db);

y_filtered=ifft(spectrum_filtered);

y_filtered=real(y_filtered);

wavwrite(y_filtered/(0.1+max(y_filtered)), SampleRate, NbitsPerSample,
         'noisy_msg6_filtered.wav');

%%%%%%%%down sampling%%%%%%%%

indexes=(abs(spectrum_freq) < 10*Freq3db);
spectrum_freq_down_sampled = spectrum_freq(indexes);
spectrum_down_sampled = spectrum_filtered(indexes);
N_down_sampled = length(spectrum_down_sampled);

spectrum_down_sampled=spectrum_down_sampled*N_down_sampled/N;

SampleRate_down_sampled=SampleRate*N_down_sampled/N;

y_down_sampled=real(ifft(spectrum_down_sampled));
t_down_sampled = ( (1:N_down_sampled)*1/SampleRate_down_sampled ).';

sound(y_down_sampled, SampleRate_down_sampled)

function [spectrum_filtered,g]=freq_space_low_pass_filter(spectrum, SampleRate, Freq3db)
%% applies low pass filter in the frequency domain
% spectrum - result of fft on time series data (column vector is expected)
% SampleRate - measured in Hz, 1/dt where dt spacing of the points in time domain
% Freq3db - desired 3db roll off point in Hz

N=length(spectrum);

function G=filter_gain(freq, Freq3db)
    G=1./(1+1i*freq/Freq3db); % this corresponds to low pass RC filter
end

spectrum_freq=fourier_frequencies(SampleRate, N);

% calculate filter gain for each spectrum frequency
g=filter_gain(spectrum_freq, Freq3db);
spectrum_filtered=spectrum.*g;
end

Plot of the signal spectrum:

A few points:

Filtering in the frequency domain is complicated and needs a real algorithm like overlap add to avoid time domain aliasing. It's much easier to filter directly in the time domain: [b,a] = butter(1,100/(SampleRate/2); y_filtered = filter(b,a,y); is much better
You probably want to put notch filter on the drone frequency
For broad-band stationary noise a good method is Wiener filtering or spectral subtraction. Many papers have been published on this.


What is the general feature of a time domain signal that gives a step like spectrum? I'm trying to study a nonlinear system by sending a pulsed signal into the system and look at the response spectrum. The signal I send to the system is something like

and the response spectrum from I get from the system is something like (note the spectrum is plotted log scale, ie, (frequency, log of absolute value of Fourier transform at that frequency))

where we can see that there is a stair like structure in the spectrum. So what is the general feature the time domain response signal should have in order to give a stair like spectrum like this? And have you ever encountered a spectrum like this, and is there something in general we can say about the system just from the stair like spectrum?

Could you give a model of this behavior that generate stair like spectrum?

If there is a time invariant linear response region in your non-linear system, then a time domain input resembling a Sinc function $\frac{\sin(x)}{x}$ in that linear region, will produce an approximation to rectangular frequency response in the frequency domain.


Early - late gate spread spectrum For situation that received signal is complex, how Early and late are calculated in spread spectrum tracking?
How I can consider both real and imaginary in error calculation?

This is the common case for receiving spread spectrum signals since we do not know carrier frequency and phase offset but can determine that from the complex output of our correlators.

Each correlator (Early, Prompt and Late) would be complex; the complex Prompt correlator output is used for Carrier Recovery (measure rate of rotation or change in phase from one complex symbol to the next), and the magnitudes of the Early and Late correlators are used for Timing Recovery.  If you are using spread spectrum that once carrier offset is removed only resides on the real axis (BPSK modulation), then the Early and Late Correlators can be simplified to be real only as long as the carrier offset is removed prior to timing recovery and carrier tracking is implemented in such a way to rotate the signal to the real axis. This is not necessary as complex Early and Late correlators can optionally be used in which case you would use the magnitudes to derive timing.

Here is a link to a decent online tutorial showing the use of the complex (I and Q) correlator outputs and further details timing recovery using Early-Late.

http://what-when-how.com/a-software-defined-gps-and-galileo-receiver/carrier-and-code-tracking-gps-and-galileo-receiver-part-3/

From which I clipped this diagram:

Additional methods to determine the error that I didn't see mentioned in the tutorial (after a quick review, perhaps they are detailed there as well):

$err = Re((P)^*(E-L))$

Where $(P)^*$ is the complex conjugate of the complex Prompt correlator output ($I_P-jQ_P$) sampled at the estimated symbol timing location, and E and L are the complex correlator outputs 1/2 a chip earlier ($I_E+jQ_E$) and 1/2 a chip later ($I_L+jQ_L$) than the estimated symbol timing location.  The operation $Re()$ is to take the real part of the result: if there is no carrier offset, the phase of P should be the same as the phase of E-L on average such that the complex conjugate multiplication will rotate the result to the real axis, however due to noise contributions or possible carrier offset the result will have an imaginary component so the Re() operation is to remove that.

Multiplying out this form of the error discriminator based on I and Q outputs of each correlator results in:

$$Re((I_P-jQ_P)(I_E-I_L-jQ_E+jQ_L))$$

$$ = I_P I_E - I_P I_L - Q_P Q_E +Q_P Q_L$$

If the signal is rotated to the real axis first using carrier recovery then this can be $P(E-L)$ were P, E and L are the output of the I (Real) correlator only in this case.

By multiplying by P as in both approaches above we both normalize the result and strip the data modulation from effecting the error result.


Proof from wave equation of upper limit of active noise cancellation Anybody know of a proof of the hard upper frequency limit of active noise cancellation using the ear canal volume as a boundary condition to solve for the upper frequency limit? I keep on encountering 1kHz as the upper limit in practice.

I don't have a mathematical proof for you, but here are 2 issues to think about.

The impulse response from the actuator (speaker) to the error measurement (microphone) is non-minimum phase. This is unavoidable because it takes a non-zero amount of time for an acoustic wave to travel from the speaker to the microphone. In a digital system there is additional latency introduced because you cant instantaneously sample, process, and output a signal. This latency limits your cancellation bandwidth.

The other issue is that even if you had a super low latency system and could get cancellation up to 10kHz it would still only be cancelling noise at the error microphone. At 10kHz your wavelength is about 3.4cm. The distance from the error microphone to the eardrum is a sizable portion of this wavelength, which means that whatever cancellation you thought you had at 10kHz is gone, and you might even be boosting now. So your perceived cancellation is limited by how accurately your microphone measurement approximates the sound field at the eardrum.


Fixed SNR with unitary noise variance I've just one question :

How can I write a model like y = x + w,( with w a white gaussian noise) with a fixed SNR and a noise variance equal to 1. What coefficient may I have before x ?

Thanks !

If you define the SNR as the ratio of the signal power and the noise power in dB, you have

$$SNR_{dB}=10\log \left(\frac{P_s}{P_w}\right)\tag{1}$$

where $P_s$ is the power of the desired signal and $P_w$ is the noise poiwer. If the noise $w$ has a mean of zero, then $P_w=\sigma^2_w=1$. From (1) (with $P_w=1)$ you get the desired value of $P_s$ for a given value of $SNR_{dB}$:

$$P_s=10^{SNR_{dB}/10}$$

In order to normalize the signal $x$ such that it has the desired power $P_s$ you first need to know its power $P_x$. Dividing $x$ by $\sqrt{P_x}$ will give you a unity power signal, which can then be multiplied by $\sqrt{P_s}$ in order to obtain the desired SNR:

$$s=\sqrt{P_s}\frac{x}{\sqrt{P_x}}=10^{SNR_{dB}/20}\frac{x}{\sqrt{P_x}}$$


Shifting signal smaller than discrete step I have a image that I need to shift with less than a pixel.

My plan was to do a Fourier transform and multiply the signal with $e^{-aiu-biv}$ where $a, b$ are the shifts in x and y direction. This might make sense in theory but in practice I don't really know how to represent $e^{-aiu-biv}$.

I am using numpy to do a 2d FFT of the image which, ofcourse, is represented as a matrix of complex values. From this point of I'm kinda lost, how do I multiply this with $e^{-aiu-biv}$?

It seems like fundamentally something is wrong as I need to temporarily represent the image in a higher resolution, shift is slightly and then collapse it back to the original resolution.

I don't do much image processing but maybe the simpler 1D case will help:

A delay in the time domain is a phase shift in the frequency domain, this notion can be used to shift your FFT in the frequency domain by a rotation much smaller than a sample

$$x(n-n_o) = W_N^{kn_o}X(k)$$

$$W_N^{kn_o} = e^{\frac{-j 2 \pi k n_o}{N}}$$

You need to apply this principle row wise and column wise. I tried giving it a go by building a 2D matrix of the x y shift and then applying it to the image. For integer shift this works, when you get into fractional shift (0-.99) for $n_o$ it gets weird.

Non-fractional shift of a "cat eating fancy ice cream" to prove that it does what I think it does:

Matlab code:

close all
clear all

%load image
I = imread('cat_eating_fancy_ice_cream.jpg');

I = im2double(I);

%make phase shift matrix for frequency domain shifting
sz = size(I)
y = e.^((-j*2*pi*[0:sz(1)-1]*100)/sz(1)).';
x = e.^((-j*2*pi*[0:sz(2)-1]*200)/sz(2));
shift = y*x;

%convert image to grayscale for easier manipulation
pF1 =(fft2(rgb2gray(I)));

%phase shift frequency domain info
pF = pF1.*shift;

%plots
figure
subplot(121)
imshow(rgb2gray(I))
subplot(122)
imshow(ifft2(pF))


How to convert dB back to manitude and then to rectangular format I am experimenting with audio (wav files) using Short-Time Fourier Transform (STFT) in Python using scipy.signal.stft.
As I understand we get data for freq, time and magnitude + phase in rectangular format.
for e.g. (shortened for illustration)

freq, time, Zxx = scipy.signal.stft(signal, fs=fs , .. and some other parms)

If I take a single value from Zxx , e.g.

x = -0.01224990-0.01877387j ,

I would like to convert this value to dB, do some modification such as

add/reduce few dB's
Also change angle or phase
and reconstruct the value back to x

So far I have,

#Magnitude
m = np.abs(-0.01224990-0.01877387j)
0.022416918717497732

#angle in degrees
ad =np.angle(-0.01224990-0.01877387j, deg=True)
-123.12434072219258

#phase
cmath.phase(-0.01224990-0.01877387j)
-2.1489251349495935

#magnitude in dB
dB = 20*np.log10(abs(p))
-32.988481658150704

#modified magnitude
dB = dB + 10.0

#Modified angle
ad = -123.12434072219258 + 50.0

Now, How do I convert it back to modified x. using above values?

I researched but couldn't find a straightforward solution. Before reinventing the wheel of constructing the complex mathematical functions, i thought I would ask the experts.

The only reference I found was mag2db and db2mag at matlab, but not sure how to access or implement these in python.

EDIT1
Please note: I am not trying to reconstruct original signal (or x value) back. I need to modify Zxx values and hence I know I will get a different signal back using scipy.signal.istft.

Although I can directly manipulate x+jy which affect both magnitude and phase, I would like to modify them individually and Hence this question.
Thanks in advance.

If $X_\mathrm{db} = 20 \cdot \mathrm{log}_{10}(|X|)$, then $|X| = 10^{X_\mathrm{db} / 20}$, of course. Furthermore, you can compute the cartesian representation of a complex number $c = a + ib$ from its magnitude $|c|$ and phase $\varphi$:  $c = |c| \cdot \mathrm{exp}(i \varphi)$

So you are probably looking for something like this:

import numpy as np
from matplotlib import pyplot as plt

fs = 100
t = np.arange(100)
x = np.sin(2 * np.pi * 3 * t / fs)

X = np.fft.fft(x)
X_mag = np.abs(X)
X_db = 20 * np.log10(X_mag)
X_phase = np.angle(X)

# ... do modification on X_db and X_phase ...

X_mag_rec = 10 ** (X_db / 20)
X_rec = X_mag_rec * np.exp(1j * X_phase)
x_rec = np.real(np.fft.ifft(X_rec))

plt.plot(t, x, label='original')
plt.plot(t, x_rec, label='reconstructed')
plt.legend()
plt.show()

If you are doing the modification in the STFT domain, note there may be no valid inverse STFT. For details on that look into:

D. Griffin, J. Lim: Signal estimation from modified short-time Fourier transform, in: IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 32, no. 2, pp. 236–243, Apr 1984.


The expectation in power spectral density I'm a bit confused with the definition of the power spectral density (PSD). From Wiki https://en.wikipedia.org/wiki/Spectral_density , I found the definition is:

$$
S_{xx}(\omega) = \lim_{T\rightarrow \infty}\mathbb{E}[|x(\omega)|^2],
$$
where $x(\omega)$ is the Fourier transform of the process $x(t)$.

I'm really confuses with that expectation $\mathbb{E}$ in the above equation. The expectation is taken w.r.t. $x(t)$, but $|x(\omega)|^2$ is a function of angular $\omega$. There is no $x(t)$ in the integrand at all, thus
$$
\mathbb{E}[|x(\omega)|^2] = |x(\omega)|^2.
$$
What is the point of this expectation?

However, in the other hand,
$$
\mathbb{E} \left[ \left  | x(\omega) \right |^2 \right] = \mathbb{E} \left[ \frac{1}{T} \int_0^T x^*(t) e^{i\omega t}\, dt \int_0^T x(t') e^{-i\omega t'}\, dt' \right] = \frac{1}{T} \int_0^T \int_0^T \mathbb{E}\left[x^*(t) x(t')\right] e^{i\omega (t-t')}\, dt\, dt' \neq |x(\omega)|^2
$$
the PSD is the Fourier transform of the cross-covariance (auto-correlation) of the process. This exectation is indeed needed.

Where am I wrong? Can't those integrals switch orders?

I think a better definition of the power spectrum is the following:

  The power spectrum of $x(t)$ is the Fourier transform of the autocorrelation function of $x(t)$, where $x(t)$ can be either a deterministic power signal, or a wide-sense stationary (WSS) random process. The definition of the autocorrelation function depends on the model for $x(t)$.

If $x(t)$ is modeled as a WSS random process, then the autocorrelation function is defined by

$$R_x(\tau)=E\big\{x^*(t)x(t+\tau)\big\}\tag{1}$$

For deterministic power signals, the autocorrelation function is given by

$$R_x(\tau)=\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{T/2}x^*(t)x(t+\tau)dt\tag{2}$$

In this answer it is shown that the following definition of the power spectrum for a WSS random process $x(t)$

$$S_x(\omega)=\lim_{T\rightarrow\infty}E\left\{ \frac{1}{T}\left| \int_{-T/2}^{T/2}x(t)e^{-j\omega t}dt \right|^2 \right\}\tag{3}$$

is equivalent to the definition of the power spectrum as the Fourier transform of $(1)$.

For deterministic power signals, the corresponding definition of the power spectrum is

$$S_x(\omega)=\lim_{T\rightarrow\infty}\frac{1}{T}\left| \int_{-T/2}^{T/2}x(t)e^{-j\omega t}dt \right|^2 \tag{4}$$

which can also be shown to be equivalent to the Fourier transform of $(2)$.

The definitions of autocorrelation and power spectrum of deterministic power signals is described in Chapter 12 of

Papoulis, A., The Fourier Integral and its Applications, McGraw Hill, 1962.

A good reference on random processes and the corresponding definitions of autocorrelation and power spectra is

Papoulis, A. and S.U. Pillai, Probability, random
variables, and stochastic processes, Boston: McGraw-Hill, 2002.


analog-to-discrete-to-analog system sampling problem This is not a homework problem, I am solving practice problems for my exam.

Consider the analog-to-discrete-to-analog system shown in figure 1. The CT signal $x_a(t)$
is sampled at a frequency of $F_s = 2000$ Hz ($T_s = 0.5$ msec). The resulting impulse train
is then converted to a discrete time sequence $x_d[n]$. The Lowpass DT filter $H_d(e^{j\Omega}$ is
subsequently used to filter $x_d[n]$ giving $y_d[n]$. Finally, a CT version of the output $y_a(t)$ is
created, using an ideal DT-to-CT converter (at the same sampling frequency $F_s = 2000$ Hz).

Note : $H_d(e^{j\Omega}$, which is obviously periodic, is shown for only one period.

a) For the CTFT of $x_a(t)$ given by $X_a(\omega)$ in the figure with $B = 2000\pi$ rad/sec, sketch the
$X_d(e^{j\Omega})$, the DTFT of the DT sequence $x_d[n]$.

b) Sketch $Y_a(\omega)$, the CTFT of the CT signal $y_a(t)$. Again, assume that we are using the
same frequency as that of sampling $F_s = 2000$ Hz.

I followed the following steps:

Step 1: Multiplied $x_a(t)$ by an infinite pulse train and then switched to frequency domain so I can have an expression for $X_d(\omega)$ in CTFT.

Step 2: After conversion, the "envelope" (which in this case is the triangular shapes) are not drawn but instead we have pulses, something like this:

(Excuse my very bad drawing skills) However, only the red lines exist in this case since we are now in discrete time.

Step 3: I believe in step 2 I have found how $X_d(e^{j\Omega})$ should look like. But then here I get stuck about how should I filter the signal. I get that the filter is periodic but how do I use this in here? Isn't the width of the filter much smaller than the width of the signal?

Your first step is correct. However, a signal which is discrete in time domain does not have a discrete frequency domain representation per se, as in theory the Discrete-Time Fourier Transform may be calculated. The signal spectrum would be discrete in combination with the Discrete Fourier Transform (DFT) only if a finite representation was required.

The main consequence of time-domain discretisation is periodicity in the frequency domain. Your solution for exercise 2 should be, as well. For this, you may need to check whether the spectrum will overlap at $\frac{F_s}{2}$.

Your assumption about the discrete filter transfer function being periodic in the frequency domain is correct. Since the filter is applied in frequency domain, you should multiply it with the signal spectrum to get the output spectrum. Should the filter be smaller than the signal, the signal components outside the passband are attenuated (or set to zero, as this filter is ideal). The computation looks like this:

$$y_d[n]=x_d[n] * h_d[n] = \mathcal{F}_\text{DTFT}^{-1}\{X(e^{j\Omega}) \cdot H(e^{j\Omega})\}$$


Accurately measuring relative distance between a set of fiducials (Augmented reality application) Let's say I have a set of 5 markers. I am trying to find the relative distances between each marker using an augmented reality framework such as ARToolkit. In my camera feed thee first 20 frames show me the first 2 markers only so I can work out the transformation between the 2 markers. The second 20 frames show me the 2nd and 3rd markers only and so on. The last 20 frames show me the 5th and 1st markers. I want to build up a 3D map of the marker positions of all 5 markers.

My question is, knowing that there will be inaccuracies with the distances due to low quality of the video feed, how do I minimise the inaccuracies given all the information I have gathered?

My naive approach would be to use the first marker as a base point, from the first 20 frames take the mean of the transformations and place the 2nd marker and so forth for the 3rd and 4th. For the 5th marker place it inbetween the 4th and 1st by placing it in the middle of the mean of the transformations between the 5th and 1st and the 4th and 5th. This approach I feel has a bias towards the first marker placement though and doesn't take into account the camera seeing more than 2 markers per frame.

Ultimately I want my system to be able to work out the map of x number of markers. In any given frame up to x markers can appear and there are non-systemic errors due to the image quality.

Any help regarding the correct approach to this problem would be greatly appreciated.

You can use structure from motion type of algorithm to estimate the camera pose from the environment, not from the markers and then fuse this camera pose with the marker poses in order to detect the locations of the markers accurately. Knowing the extrinsic pose of your camera (by SFM), you can triangulate all the 3D positions.

For pose estimation, 5-point methods are usually more accurate than 8-point algorithm.

Presumably, you should do further bundle adjustment so that the overall accuracy increases.


Can you downsample unevenly with even distribution? Could I downsample a series of length 10 to length 7, with even distribution?
i.e. Numbers get split across two (neighbouring) resultant numbers, as opposed to one result being a combination of 1 number, and another result being a combination of 2.
i.e. All numbers be a combination of 10/7 numbers.
For this example, with my current - naive - approach, some of the 7 "buckets" would average 1 number, and some other buckets would average 2 numbers: avg(1), avg(2, 3), avg(4), avg(5, 6), avg(7), avg(8, 9), avg(10)
Would the solution be something like: weighted_avg((1, 1), (2, 10/7 - 1)), weighted_avg((2, 1), (3, 10/7 - 1)), ...
where weighted_avg = ((value_x_1, weighting_1), (value_x_2, weighting_2)) -> weighted average

Absolutely positively -- uh -- maybe.
It depends entirely on the underlying process that generated the original series, what the series "means", how much unique information is actually present in those ten samples, how much you're willing to throw away, and what you know about how that information is structured.
At the extreme "no" end would be if those ten data points were essentially completely random -- i.e., if you recruited ten people on the street and asked them their favorite number, then there's no way of relating the 1st sample to the second.
The trivial extreme "yes" would be if those ten data points were super highly correlated -- the altitudes of ten points a meter apart on a nice smooth sidewalk, for instance.  Then your proposed simple interpolation scheme would work quite nicely.
What you're trying to do is a generalization of downsampling.  If you had an infinite sequence of points that it made sense to describe in the frequency domain, then as long as the sequence was band-limited to less than 7/20th of the original sampling rate, you could upsample by a factor of 7, filter the resultant signal, then downsample it by a factor of 10.
You can't do that, because filters suffer from end effects, and from the point of view of the necessary filters, your data sequence is all end and no middle.
The absolute best formal way to do the job would be to understand whatever the underlying process that generated those ten points.  From that, figure out how to work backward to the most likely location or the seven points, and voila, you either have an answer or your loved ones retrieve your body from under a mountain of scribbled on paper and give you a decent burial.
A pretty good ad-hoc way to do this, if you have reason to believe that those ten points represent some sort of a bandlimited process, is to just use a set of splines to make a smoothed line between the ten points, and then sample that line at the seven points.  This really only works well if the underlying process happens to match up well with a spline -- but if it's bandlimited it does -- and if it's something like the positions of a beam at those ten points then it definitely does, because the mathematical sort of spline is a pretty good simulation of the draftsman's spline, which is a beam that's pinned down at certain points to make a sparse set of points into a nice smooth curve, for boatbuilding.


Smoothing a Signal which contain White Noise I'm trying to smooth a signal in order to reduce white noise as a result of changes in amplitude across sequential frames of an audio signal.

I've been investigating using a moving average or boxcar filter to help minimise the issue.  In my investigations it would appear normal implementations of a boxcar filter only smooth the contents of a single frame/buffer of audio.

Can a boxcar filter be implemented in order to smooth the transition across frames as well? or is there a better approach?

A boxcar filter is a FIR filter, suppose your frame has length N and your boxcar has a length of M, you will have M+N-1 filtered values.There are 2 properties, time invariance and superposition that you can use.

Assume your boxcar has length 5, and examine 2 sequential frames,data is $x$

$$
\begin{align}xxxxxxxxxxxxxxxxxxxx&0000  \quad \text{Frame 1 with 4 zeros post appended} \\
                0000&xxxxxxxxxxxxxxxxxxxxxxx \quad \text{Frame 2, 4 zeros pre appended}
\end{align}
$$

If we could time align the 2 frames, and added them together, we would have continuity across the frames.  This is a consequence of superposition.

If we run the boxcar over each frame,including the zeros
$$
\begin{align}yyyyyyyyyyyyyyyyyyyyy&yyyy \quad \text{Filtered Frame 1 } \\
                yyyy&yyyyyyyyyyyyyyyyyyyyy \quad \text{Filtered Frame 2 }
\end{align}
$$

Adding the 2 time aligned filtered frames, is the same as filtering a gapless buffer.  This is superposion, and time invariance of the filter.


Discrete Fourier Transform by hand I have an assignment where I'm given the DFT of a sequence $x[n]$ as $X[k]=\{4,3,2,1,0,1,2,3\}$ and also
$$y[n] = \left\{ \begin{array}[cc]
xx[n/2] & \text{if n is even} \\
0 & \text{otherwise}
\end{array} \right. $$

and I'm supposed to find and sketch the DFT of $y[n]$.

So $y[n] = \{x[0], 0, x[1], 0 ... x[7], 0\}$ and it's not complicated to find $Y[k]$ if we know $x[n]$

I know how to use the definition of the DFT and IDFT to calculate $x[n]$ but it's a tedious task to do by hand, especially when the sequence is longer than a few items. Is there a quicker way to calculate the DFT and IDFT by hand without using a program like Matlab?


Signal processing using panda board I'm an absolute beginner.
I'm doing a project that involves estimating the time delays between the sound signal arriving at different microphones in an array (the project is a bit complicated but I'll find my way through). I'm planning to use a panda board and connect the microphone array to it. But before going about implementing the array, I thought I should first learn how to do some DSP on an analog input. Basically, how would I connect ONE microphone to the panda board and do, for eg, an FFT on the signal?

The panda board has an audio in port but that wouldn't help much when I have to use multiple microphones. Can I connect the condenser mic (the really cheap one ) using 2 wires to a general input pin? How would I read input from it? Will I be able to record sound using such a setup? And I understand that I have to do the Analog to digital conversion (In fact I WANT to do it by myself, so that I would better understand how stuff work).

And now the million dollar question : If I connect a few more microphones to different input pins of the panda board and give them the same sampling rate (I have no idea how but I'll learn it) will that make a microphone array?

To connect one microphone to the Panda board, use its audio input. The homepage says it's a stereo input, so yo can actually connect two microphones assuming your microphones are mono. To apply any DSP on the analog signal it must first be converted to a digital signal, of course. For that you need an analog-to-digital converter (ADC). It's an electronic device. You cannot do analog-to-digital conversion by yourself. I suspect the Panda board audio input to have a dual-channel ADC, otherwise it would be quite useless.

I haven't dived into the technical specs of the Panda board but I'd be very surprised if its general input pins were made for analog signals. So you cannot use those for recording audio, because they're lacking ADCs (unless you manage to use the ADCs of the RF interfaces present on the OMAP platform). For each microphone you need an ADC and an interface to the Panda board. Any selfmade solution will be a real hassle here, I would not recommend it.

A possible alternative solution is to use a professional audio interface that has multiple microphone inputs and a USB output. Then connect this audio interface to the Panda board via USB. The challenge here is to get a driver running on the Panda board. But I think it should be feasible. Keep Linux support in mind when choosing an audio interface.

If you've managed that, and you arrange your microphones in a reasonable way you have a microphone array.


Time shifts and change in phase My question here is in regards to implicit information/terminology.

Basically, whenever we say that "a shift in time means a change of phase in fourier space", is it implicit that the shift is always circular? Are there any exceptions to this? Caveats? etc?

For a continuous-time signal $x(t)$ with Fourier Transform $X(f)$, the
Fourier transform of $x(t-t_0)$ is $e^{-j2\pi ft_0}X(f)$, that is, the value
of $X(3)$, say, gets changed to $e^{-j2\pi 3t_0}X(3)$ which has the same
magnitude as $X(3)$ but a different phase.

For a discrete-time continuous-amplitude signal $x[n]$ with DFT $X[k]$
of length $N$, you have essentially discarded every piece of information
conveyed to you by $x[-1], x[-2], \ldots $ and by $x[N], x[N+1], x[N+2], \ldots$
and have retained only $x[0], x[1], \ldots, x[N-1]$. If you phase-shift
the DFT appropriately, then yes, you will get a circular shift of the $x[n]$
values. The result will be the linear shift of the $x[n]$ values only
in the fortuitous circumstance that the sequence $x[n]$ is itself periodic
with period $N$, and so the untrained observer such as myself cannot tell
if that symbol that just shifted into my DFT wndow is $x[-1]$ coming from
the left (linear shift) or $x[N-1]$ from the right (circular shift)
because by golly they look exactly the same to me.  Savvier denizens
of this forum will give you more elaborate answers explaining how you
can tell the difference between the two cases.


How do I explain a complex exponential intuitively? What is a complex exponential, explained intuitively?
How do I explain to an adolescent a complex exponential function?

From a graphical point-of view, it is  an infinite spring, whose distance between adjacent coils reflects the frequency of the complex exponential:

If you have a 1D time x-axis, you may be used to draw functions along a single 2nd y-axis  dimension: sines, cosines, etc. If you want to plot a complex function,  you need one x-axis, and 2 y-axes for the real and imaginary parts. You can plot this into a 3D visualization, see above, and you see a spring (although it is some hot summer here). Changing the frequency dilates or contracts the spring. It is known as the Heyser corkscrew, or Heyser spiral, in DSP.

A little more details can be found there:

3D wiggle plot for an analytic signal: Heyser corkscrew/spiral
DSP Concepts Visually Explained

The Fourier transform tells you that any wire (a function) can be reproduced by a superposition of scaled and shifted springs.

After the what, the why. Complex exponentials (or cisoids) are special in  that if one is filtered (with a moving average) it keeps the same shape. So, they are invariant, under  Linear-Time-Invariant (LTI) systems. Invariant vectors/functions are often an appropriate way to study systems or transformations. Morever, as they form an orthogonal basis, they form a basis of choice to decompose arbitrary vectors, and to study how the latter are affected by LTI systems.

Last, the complex exponential is itself invariant under differentiation ($(e^{z})' = e^{z}$), a specific linear and invariant operator), which makes it quite unique, with interesting properties.


Is it possible to know the order of the filter, just looking on the pole zero plot? is it possible to know the order of the filter, just looking on the pole zero plot.

I know how to get the order of the filter using calculations(highest order), but I want to know is it possible to see that in pole zero plot.

So far my idea is that the higher number of zeros/poles (pairs) is the order, for example if I have 3 poles and 1 zero in a zero ploe plot, than is the filter 3 order.

My knowledge over this is still small, so I am sorry if it is a stupid question.

Thanks

Count the number of poles and zeros; if the numbers aren't the same then the larger number is the filter order. Don't forget to count multiple poles or zeros with their multiplicity. Actually, you always get the same number of poles and zeros but on a pole-zero plot you might not see them all because some could be at infinity or at the origin.

E.g., in the following pole-zero plot you see 12 zeros, 6 poles away from the origin, and a pole at the origin with multiplicity 6. So the total filter order is 12. (This is an IIR low-pass filter with approximately linear phase in the passband):


DTFT of $(-1)^{n}\frac{\sin((7\pi / 8 )n )}{\pi n}$ The DTFT without $(-1)^{n}$ is relatively easy but with $(-1)^{n}$ I am not sure how it would be solved. I hope anyone can help

HINT:

$(-1)^n=e^{jn\pi}$. What is its DTFT?
Multiplication with $(-1)^n$ corresponds to convolution in the frequency domain with the DTFT you obtained in step 1.


How to Calculate Total Variation (TV) of an Image? TV is L1 norm of gradient of an image. So we've to find gradient of the image (which is still matrix, right?). Then take the sum of absolute values of the gradient matrix (So now it must be a scalar, right?).

Can I simply use imgradient function in matlab for first step? It gives two values magnitude and direction, so which one should I consider for next step?

After getting gradient matrix say Igrad, will just summing absolute values (sum(abs(Igrad(:))) give me TV? And it must be just a single scalar value right?

Thank you for reading. Keenly waiting for your reply.

The Total Variation of an image $I$ can be computed using two formulas:

$TV(I) = \sum_{x} \| \nabla I (x) \|_1 $ (anisotropic TV);
$TV(I) = \sum_{x} \sqrt{ \| \nabla I (x) \|_2^2 }$ (isotropic TV).

In practice, both formulas yield almost the same result.

Using matlab, what you propose implements the first formula, while the second one can be obtained by sqrt(gradI_x.^2 + gradI_y.^2).


Write a 440 Hz Sine Wave to .wav-File using Python and SciPy This is my first post on DSP.SE! I'm a student in Computer Sciences and am just getting started with Signals and Systems, thus getting in touch with discrete signals and analyzing them.

I wanted to explore this world a little further and so my first goal is to create a Sine Wave that has a Frequency of 440 Hz and write it to a .wav-File. Of course i have all the calculus in mind, but i can't seem to figure out the problem myself!

So far, i have this:

#!/usr/bin/env python3

import numpy as np
from scipy.io import wavfile

sampleRate = 44100

t = np.linspace(0, 20, sampleRate * 5)  #  Produces a 5 second Audio-File
y = np.sin(440 * t)  #  Should have frequency of 440Hz

wavfile.write('Sine.wav', sampleRate, y)

I understand why this code always creates a File, that is 5 seconds long, but i can't figure out, how to get the frequency right!

Have i even gotten the rest of the code right, Dsp-Guy-wise?

After some experimenting i came up with the following:

#!/usr/bin/env python3

import numpy as np
from scipy.io import wavfile

sampleRate = 44100
frequency = 440
length = 5

t = np.linspace(0, length, sampleRate * length)  #  Produces a 5 second Audio-File
y = np.sin(frequency * 2 * np.pi * t)  #  Has frequency of 440Hz

wavfile.write('Sine.wav', sampleRate, y)

I pulled out one of my math-books and figured out how to write the sine-function:

Since sin(x) has a period of 2*pi, sin(2*pi*x) has a period of 1/Frequency of 1Hz. Thus, sin(440*2*pi*x) has a Frequency of 440Hz.

The last thing that was unclear to me was how to get the array-length and element-count right!

Opting for a File that has a sample rate of 5 seconds at 44.1 kHz Sample Rate, i need 5 * 44.1k samples, so that was already ok. But since sin(440*2*pi*x) performs 440 periods until t == 1, the array simply has to have a duration of 5 (5 seconds in my case to be precise). enlarging the array boundaries would have squished the array-values together, resulting in a higher pitched sound.

Sorry to anyone that has already scratched their head, i probably should have thought for a little bit longer!


Can you simulate the "response time" of a lower window size FFT on a higher window size? Can you simulate the "response time" of a lower window size FFT on a higher window size?

That is, is it possible to calculate the equivalent lower window sized FFTs from a higher one in order to iterate over the equivalent e.g. 128 window sized FFTs, but plot the higher one e.g. 8192 samples window size?

This is a problem in e.g. when one is sampling at a higher window size for the FFT, but some signal parts (e.g. strong amplitude peaks in audio) appear at a faster pace and are thus "missed" by the FFT. Or not necessarily the FFT, but any processing that uses that higher window size FFT.

What I therefore need is the visible resolution of a higher window size and the faster response of a lower window size.

My application (in audio) requires me to be able to track frequencies that follow closely to the beats of the audio. If e.g. a drum loop plays at a high bpm, then a 8192 samples FFT might be too slow to be done reading form since it changes slower than the drum loop's frequencies.

No. The time of accumulation will determine the resolution of frequencies - how far are measures separated.