Decimating a discrete time cosine function Let $x[n] = \cos(wn)\ \ $ be our signal. Let the other signal be $x_2[n] = \cos(2wn)\ \ $ . The spectrum of $x[n]$ and $x_2[n]$ are well defined and consist of impulses at the relevant frequencies with weight $\pi$. However, if we interpret $x_2[n]$ as the decimated version of $x[n]$, the spectrum of $x_2[n]$ must have impulses with weight $\frac{\pi}{2}$. I have the explanation but I would like to hear the alternatives, what is the best way to resolve  these seemingly contradicting facts?

The following correction of Jason R's answer completely resolves the problem: The frequency spectrum of $X_2(ω)$  consists of two components, which are scaled versions of $X(ω)$. Those are $4\pi$-periodic instead of being $2\pi$-periodic. When those $4\pi$-periodic versions are added up, the sum becomes $2π$-periodic again. THERE IS NO ALIASING or IMPULSES SUPERIMPOSED ON EACH OTHER. The derivation from
$$
X_2(\omega) = \frac{1}{2} \pi\left(\delta\left(\frac{\omega}{2}-\omega_0\right) + \delta\left(\frac{\omega}{2}+\omega_0\right) + \delta\left(\frac{\omega}{2}-\omega_0-\pi\right) + \delta\left(\frac{\omega}{2}+\omega_0-\pi\right)\right)
$$
to
$$
X_2(\omega) = \frac{1}{2} \pi\left(\delta\left(\omega-2\omega_0\right) + \delta\left(\omega+2\omega_0\right) + \delta\left(\omega-2\omega_0-2\pi\right) + \delta\left(\omega+2\omega_0-2\pi\right)\right)
$$
is not correct. Instead it must read
$$
X_2(ω)=π(δ(ω−2ω_0 )+δ(ω+2ω_0)+δ(ω−2ω_0−2π)+δ(ω+2ω_0−2π)) \qquad (∗)
$$
This is due to the scaling property of impulses
$$
δ(aω)=\frac{δ(ω)}{|a|}
$$
Then, observing that (*) is $2π$-periodic, two of the impulses are redundant and hence the DTFT of $x_2[n]=x[2n]$  becomes
$$
X_2(ω)=π(δ(ω−2ω_0)+δ(ω+2ω_0))
$$


Housekeeping operations in filters I have to implement a second order filter in DF1 or DF2. I should avoid housekeeping operations by suitable addressing strategies.

I don't understand what "housekeeping operations" means. Anyone could help me or propose a source i can find more information?

EDIT :

I use DF2 for the second-order filter:

As pichenettes proposes with Housekeeping the code in matlab is something like:

%x=input sample
%state=containing old input values
%compute output value
y=b(1)*x+state(2)

%Update states
state(2)=b(2)*x-a(2)*y+state(1)
state(1)=b(3)*x-a(3)*y

Without Housekeeping is going to be like:

%compute output value
y=b(1)*x+state(pointer+1);
pointer=rem(pointer+1,N-1); %increment pointer in modulo form, N number of taps

%Update states
state(pointer+1)=b(2)*x-a(2)*y+state(pointer+1);
pointer=rem(pointer+1,N-1);

%Overwrite oldest sum with b(N-1).x
state(pointer+1)=b(3)*x-a(3)*y;
pointer=rem(pointer+1,N-1);%Increment pointer modulo-(N-1)

Is that right?Any mistakes?

First, maybe you could give us more context (is that homework?).

I think that what is meant by "housekeeping operations" are the data transfers between registers / state variables ("delay propagation") that naturally arise when naively implementing filters. For example, a naive implementation of a 4-tap FIR would look like this:

output = h[0] * input + h[1] * x[0] + h[2] * x[1] + h[3] * x[2]
// "Housekeeping" operations
x[2] = x[1]
x[1] = x[0]
x[0] = input

The data transfers operations can be avoided on some hardware architectures (especially DSP) by using modulo addressing. I think this is what is meant by "suitable addressing strategies".


Where did Arnold Tustin first introduce the bilinear transform? It's well known that the bilinear transform is also known as Tustin's Method.  As far as I know, Arnold Tustin really did introduce the idea into the control systems literature, so the name isn't just a case of Stigler's Law.  For example, I managed to find the following reference:

  Tustin in the UK developed the bilinear transformation for time series
  models, while Oldenbourg and Sartorius also used difference equations
  to model such systems. $[1]$

What isn't clear is where he first introduced the idea--even when browsing the titles of his publications.  I'm guessing that it only became known as the bilinear transform later on, so he probably didn't use that terminology.  I'd like to read his exposition of the technique. Does anyone know where he first published it?

Bissel, C.C. A History of Automatic Control. link.

Even though I couldn't find a decent table of contents, I took a shot on ordering a source book of papers edited by Tustin (i.e. Automatic and Manual Control), thinking he may have included something relevant.  The direct source wasn't in there, but in a paper by B. M. Brown, Application of Finite Difference Operators to Linear Systems, the following citation is given:

  The simplest of these [approximations in the form of rational
  fractions] has been derived ad hoc by $Tustin^6$.  It is the time
  series equivalent of the operator   $$ {2\over{\delta}} \space {1 - E^{-1} \over {1 + E^{-1}}} $$

  This can be written as

  $$ {2\over{\delta}} {E - 1 \over{ E + 1 }} \space = \space {1\over{\delta}} {\Delta \over { 1 + 1/2 \Delta} } \space = \space {1\over{\delta}} (\Delta - 1/2 \Delta^2 + 1/4 \Delta^3  - \space \ldots \space) $$

That looks like the bilinear transform. The actual reference is a little vague,

$$
^{6} \space Tustin, A. \space \space J. Inst. Elect. Engrs. \space
94 \space (1947) \space 130.
$$

but I'm pretty sure it corresponds to A method of analysing the behaviour of linear systems in terms of time series.  That paper has a pretty detailed exposition, and Tustin doesn't seem to cite any of his other works as references, so it's reasonable to assume that's where he first introduced it. What's really interesting is the notation he used,
$$
{2\over{\delta}} { [1,-1] \over{ [1, 1] } }.
$$


How does chessboard configuration affect camera calibration? When calibrating a camera with a chessboard, are more squares better or worse? When would you use more or less? And how does the size of the blocks effect the calibration results?

Generally speaking, when you are estimating a model from data, the more (independent) data points, the better the estimation.

In this particular case, the goal is to get as many 3d -> 2d point correspondences as possible, filling as much as possible of the volume of space of interest that is observed by the camera. The size of the squares directly affects how many data points you get per image (larger squares implies less visible corners/points).

It is also important to ensure that the camera parameters you want to calibrate are observable given the data. For example, in order to calibrate the focal length (or the field of view, equivalently) you need to observe a significant perspective foreshortening of the (planar) checkerboard in a few images. In order to accurately calibrate the nonlinear (e.g. barrel) lens distortion, you'll need to observe plenty of aligned data points near the edge of the image, where the distortion is larger.

There is a longer response on StackOverflow that you may want to look up.


How to implement Viola-Jones algorithm in Matlab? I am new in Image Processing and I want to implement Viola-Jones algorithm along with Haar Features for face detection in Matlab. Is there any proper documentation for this?

Somebody did it for you:

Viola-Jones Face Detection for Matlab
Viola Jones Object Detection

and I guess you can find a few others. A nice description, in pseudo-code, can be found in An Analysis of the Viola-Jones Face Detection Algorithm, IPOL, 2014, which you can follow to code your own.


Why is there a time-frequency precision tradeoff I understand the science - I think. Looking at a spectrogram, the lower frequencies are fuzzier due to a lower ratio of: frequency/time.
My main question is: why can I as a human hear lower frequencies sharply, but the FT cant?
Is this a psychological thing - perhaps my ears suffer from the same issue, but my brain is filling in the missing information?

Notice how the band gets wider, demonstrating the increase in innacuracy.

why can I as a human hear lower frequencies sharply, but the FT cant?

No you can't. Human pitch discrimination at low frequencies is actually not particularly good. The "sweet spot" for human pitch perception is between 250 Hz and 4kHz and it gets rapidly worse above and below that. When tuning a bass instrument you typically try to reduce the fundamental and the lower harmonics as much as possible and emphasize the higher harmonics a lot since this makes it a lot easier to determine the correct pitch. A violin slightly out of tune makes everyone cringe: a double bass slightly out of tune is hardly noticeable.
So in this regard the human auditory system is no better than the Fourier Transform.


Removing noise from F2F signal (This question relates to Extracting Binary Magnetic-Strip Card Data from raw WAV)

I am extracting the binary sequence from the magnetic strip on a credit card.

As you can see, the signal is degraded clearly in one place.   also there is a minor degradation right on the left of the image.

Just using IIR (i.e. $X_{\rm out} = 0.9X_{\rm out_{last}} + 0.1X_{\rm in}$) smooths it, but the resulting signal is not mathematically smooth;  if I differentiate the signal a couple of times the noise comes back with a vengeance:

My question is: can I remove the noise in such a way that the derivatives come out clean?

If so, how?

EDIT: Here is a close-up of some damaged waves:

EDIT(2):
A couple of approaches I am considering:

Firstly I could make a taylor approximation of the signal either side of the damaged sector, and blend the approximations together.
Secondly I could FFT, remove high-frequency components and reverse FFT. I'm going to try that second approach now...

To remove high-frequency "noise" without removing sharp transitions, you might have to try some sort of non-linear filtering process.

A random example might be applying (cross fading to) a median filter only when the distance to the nearest N local maxima (or minima, and beyond some noise threshold) goes under some distance threshold (where this distance, N and the median filter width, are calibrated to the period of the suspected "noise").


Having problems with understanding the meaning of PSD for colored Gaussian noise After reading two articles on signal processing stack exchange:
On coloured Gaussian noise
How the white and colored noise differ in time domain
I do understand that variance do not change over frequency (or time) for colored Gaussian noise.
Let's say I have a random process $y_n$, which is defined as:
$$ y(n) = x(n-1) - x(n), \quad\text{where }x_n\text{~} N(0, \sigma_n^2)
$$
For the clarification, $x_n$ is white Gaussian, and hence $E[x_n\cdot x_m]=0 \text{ for } n\neq m$.
Here, we can see that PSD of $y_n$ has high-pass form since
$R_{yy}(\tau)=2R_{xx}(\tau)-R_{xx}(\tau-1)-R_{xx}(\tau+1).$
we can also easily see that its variance is $2\sigma_n^2$.
Next, we also introduce another random process which has a 2x higher variance than $x_n$:   $w_n\text{~}N(0,2\sigma_n^2)$
Now, I would like to compare the random process $y_n$ to $w_n$.
First, we know that their variance are identical to $2\sigma_n^2$.
Second, although $y_n$ is colored Gaussian and $w_n$ is AWGN, their PDF must look the same, because they are still Gaussian process.
What I do not understand clearly is the meaning of PSD of $y_n$. As shown below, plot with a red color is a PSD of $y_n$. It shows high-pass shape, start having higher power magnitude after passing $0.5 f_s$. But what does it mean for a Gaussian noise having higher power magnitude at high frequency and having lower power magnitude at low frequency?
It cannot be that $y_n$ has higher "variance" at high frequency and lower "variance" at low frequency because I mentioned earlier that variance should not vary over frequency or time.
The next thing I thought about is the difference between certain sequences of $y_n$ over time. For instance, if $y_n$ has a high frequency pattern (i.e., transition density of 1) for $n=1\cdots 100$, then it may temporarily have a higher variance for this particular high frequency pattern. But soon I realized that this argument does not make sense neither since Fourier transform does not care about behavior of $y_n$ for some local pattern.
If anyone can give me any good insights about the meaning of PSD for the example above, I would be greatly appreciated!
Please let me know if anything is not clearly in the above. English is not my first language so there is a very high chance I did not explain things very clearly.
$Y(f)$ vs $W(f)$" />

Think in terms of noise realizations. When your head is full to the brim of covariances and multivariate Gaussian variables, generate a Gaussian white noise realization of length N, and create, using this sequence, a sequence $y(n) = x(n) - x(n-1)$. Let us call the latter a blue Gaussian noise realization. I plot histograms of these sequencies which show how many values fall into each bin, those bins uniformly partitioning a range of $y, x$ values.

In compliance with your convention, a blue colored graph represents the white noise data, a red colored, the blue noise data. This histograms visually agree with your observation that

$w_n \sim N(0,2σ_n^2)$, 2x higher variance than $x_n$

Each histogram plotted here is an approximate representation of the distribution of noise signal data, id est, the Probability Density Function of noise signal. The realization sequence's length N is 220 samples. The longer sequence length gets and the greater bin size gets, the closer these numerical PDF graphs resemble theoretical curves.
Go on to Power Spectral Density. Here are discrete Fourier transforms (actually, FFTs) of our white-noise and blue-noise realizations, squared and plotted. Color scheme according to your convention (white noise in blue, blue noise in red):

A white-noise PSD (blue graph) is near-constant, a blue-noise PSD graph (red graph) is two symmetrical parabolic segments sewn at 1/4 and 3/4 of the entire range (the shapes are revealing $\omega$ squared function and FFT mirroring). The high-frequency end is in the middle of horizontal axis.
No relationship between variance and PSD is revealed in our calculations and pictures. It appears that these concepts are things that can only be related indirectly, using formulas that have yet to be derived at this stage of learning.
Variance 'measures how far a set of numbers is spread out from their average value' (Wikipedia). Variance is defined in the time/space domain, and this definition never mentions frequencies or wavenumbers or Fourier transform. Although variance can be measured for any "set of numbers", in many signal processing scenarios said set of numbers is almost synonymous with sampled noise values.
A PSD thing, on the contrary, belongs in the frequency domain. The article in Wikipedia states from the beginning

According to Fourier analysis, any physical signal can be decomposed
into a number of discrete frequencies

and continues

PSD of the signal describes the power present in the signal as a
function of frequency. PSD of noises is a specific case of PSD of
signals in general.

So variance and PSD are sort of "orthogonal" concepts; whatever, variance cannot be a function of frequency in any meaningful way, while PSD is.
Compare also the histograms of the white-noise realization PSD (top) and the blue-noise realization PSD (bottom). The white-noise PSD shape replicates, as it should, a normal distribution of source signal values. A peak of the blue-noise PSD shape shows greater share of the high-frequency components in the blue noise signal.

PS There is a number of Guassian Process tutorials on github and elsewhere, mostly Jupyter notebooks. See, for example, https://peterroelants.github.io/posts/gaussian-process-kernel-fitting/.


Convolution a 3-D signal with 2-D signal I have a 3D signal $M(x,y,t)$ and I want to do convolution with signal $G(x,y)$ and $Et(t)$,

I mean, how to calculate this Convolution below in Matlab?

$$A(x,y,t)=M(x,y,t) \star G(x,y) \star Et(t)$$

In my case, $M(x,y,t)$ is a movie with $k$ frame; $x,y$ are position in image and $t$ depend on frames number and frame rate.
$G(x,y)$ is a gaussian shape with $x,y$ and $Et(t)$ is a simple exponential in time.

I propose this is an abuse of notations, which can be tackled as follows: the convolution in the time dimension $t$ is independent of the one in the space dimensions $(x,y)$, in a separable way. In practice, you can convolve each 2D data $A(x,y,.)$ with $G(x,y)$, then each resulting 1D data $A'(.,.,t)$ with $Et(t)$, or in the opposite order, which is allowed due to linearity.

A cleaner version would construct the 3D operator by a tensor product operation between $G$ and $Et$, resulting in:

$$ M\ast (G\otimes Et)\,.$$

This is sometimes called a (separable) space-time filter. It is simpler to implement, supposing that the perception rules differ in space and time. However, it reached its limits in following fast motions between time frame (sometimes, non-separable filters or full 3D filters are used in video).

Basically, it combines 2D and 1D filters into a genuine 3D filter $H_{G,Et}(x,y,t)$, so that, for each $(x,y,t)$:

$$H_{G,Et}(x,y,t) = G(x,y)\cdot Et(t)\,.$$

In practice, using separability is more efficient. In Matlab, if you can load all the data, if you can for instance first 2D filter each 2D frame with $G$, then "transpose" to get the time in the first dimension, and apply:

  y = filter(b,a,x)

which is vectorized, and revert to the initial data order.

Outside this natural interpretation, I consider the above  abuse very abusive.

Additional reading on separability:

How does a separable filter work?
Why does the separable filter reduce the cost of computing the operator?
How to prove a filter is separable?
Filtering $n×n$ images by separable $m×m$ filters. Computation time for filtering using FFT, 2-D convolution, and two 1-D convolutions


Reproducing MATLAB's randn Gentlemen! I have to reproduce MATLAB random normal generator (randn) with some fixed seed.

Also, I have to reproduce it in an external script, without MATLAB itself. Has anyone done this? What random normal generator implementation is better to use?

P.S. Preferred languages are Python and C/C++, but other languages may be useful too.

Its hard to say without knowing exactly how Matlab's randn is processed from the RNG you're using and how Matlab uses the seed.

The most reasonable way is  to save pre-calculated noise data or just deal with a different random number generator (Boost among others has a normal random number generator, and you can pull it into Matlab to use via MEX if necessary).


How to derive the transmitted signal from a constellation diagram Let us assume we have $16$ symbols to transmit.  We can represent these $16$ symbols by $16$ unique signals and transmit.  If at the receiver we can identify these $16$ signals correctly, we have identified the symbols transmitted correctly.  I have used qammod([0:15],16,0) in MATLAB and got the output that is attached here.

>> a = qammod([0:15],16,0); a'

ans =

  -3.0000 - 3.0000i
  -3.0000 - 1.0000i
  -3.0000 + 1.0000i
  -3.0000 + 3.0000i
  -1.0000 - 3.0000i
  -1.0000 - 1.0000i
  -1.0000 + 1.0000i
  -1.0000 + 3.0000i
   1.0000 - 3.0000i
   1.0000 - 1.0000i
   1.0000 + 1.0000i
   1.0000 + 3.0000i
   3.0000 - 3.0000i
   3.0000 - 1.0000i
   3.0000 + 1.0000i
   3.0000 + 3.0000i

$0,1,\ldots,15$ represent the $16$ symbols. $16$ represents the number of symbols or the number of unique signals - sinusoids with different amplitudes and phases - to be transmitted.  $0$ represents the offset phase.

Now I want to interpret the output.  I actually got $16$ complex numbers.

What do they represent?
What are the magnitude, phase, real part and imaginary part of these complex numbers that they represent?  I guess they represent the parameters of sinusoidal signals.
In both the real and imaginary parts we see $-1, -3 , 1$ and $3$.  We do not see $-2, 0$  and $2$.  Why?

This is a basic question about how passband pulse-amplitude-modulation (PAM) works. It has nothing to do specifically with QAM, but it applies to any type of passband PAM (such as PSK, or any other choice of a constellation).

The question is how to obtain a transmit signal from a given point in a constellation diagram. Let's say that at time $t=kT$ we want to transmit a signal corresponding to the point $a_k$ in the constellation diagram, where $a_k$ is a complex number:

$$a_k=\text{Re}(a_k)+j\text{Im}(a_k)=|a_k|e^{j\phi_k}$$

The complex number is multiplied by a (real-valued) transmit pulse $g(t)$ and by a complex carrier, and the transmitted signal is the real part of that complex-valued signal:

$$\begin{align}s_k(t)&=\text{Re}\left(a_kg(t-kT)e^{j\omega_0t}\right)\\&=\text{Re}(a_k)g(t-kT)\cos(\omega_0t)-\text{Im}(a_k)g(t-kT)\sin(\omega_0t)\\&=|a_k|g(t-kT)\cos(\omega_0t+\phi_k)\end{align}\tag{1}$$

In practice you transmit a sequence of symbols, and the corresponding transmitted signal is given by

$$\begin{align}s(t)&=\sum_ks_k(t)\\&=\cos(\omega_0t)\sum_k\text{Re}(a_k)g(t-kT)-\sin(\omega_0t)\sum_k\text{Im}(a_k)g(t-kT)\\&=\sum_k|a_k|g(t-kT)\cos(\omega_ot+\phi_k)\end{align}\tag{2}$$

As you can see from $(1)$ and $(2)$, the transmitted signal can be represented either as two amplitude-modulated orthogonal carriers (sine and cosine), or as an amplitude and phase modulated carrier. In the first case the real and imaginary parts of $a_k$ determine the amplitudes of the two carriers, and in the second case the magnitude $|a_k|$ determines the amplitude of the carrier, and $\phi_k=\arg\{a_k\}$ is the phase of the carrier.


Is the below filter linear phase $$h(t) = \frac{1}{1 + t^2}$$ and is it IIR or FIR filter. I tried finding the Laplace transform of this filter to get the data flow diagram with 5 taps and T=2s, however, I am unable to solve this. Probably the solution is simpler than I tried, however, being novice in DSP, I am unable to solve the same.

Thank You.

If

$$h(t)=\frac{1}{1+t^2}\tag{1}$$

denotes the filter's impulse response, then the filter is an infinite impulse response filter, simply because $h(t)$ has infinite support, i.e. it doesn't vanish anywhere (other than for $t\rightarrow \infty$). Note that if the definition $(1)$ of $h(t)$ is valid for all $t$, then the filter is not causal because $h(t)=0$ for $t<0$ doesn't hold.

Since the impulse response is real-valued and it satisfies $h(t)=h(-t)$, its Fourier transform must be real-valued (and symmetrical). Since the Fourier transform is real-valued, the filter's phase is zero.

If you sample that impulse response using $5$ taps at times $t_{-2}=-2T$, $t_{-1}=-T$, $t_0=0$, $t_1=T$, $t_2=2T$, you get an FIR filter, but the zero-phase property remains unchanged due to the symmetry of the filter coefficients. For implementation, the filter impulse response needs to be shifted by two samples to make it causal. This causes a delay of two samples, which changes the filter's phase $\phi(\omega)$ from zero to a linear phase $\phi(\omega)=-2\omega$.


Unable to understand the time-shifting property of CTFS The CTFS of $x(t)$ is $c_{k}$ the Fourier series coefficients. Furthermore, $x(t-t_{0})$ is known to be $e^{-j\omega t_{0}}c_{k}$, the proof is given as follow :
$$
\begin{aligned}
\mathscr{F}\left(f\left(t-t_{0}\right)\right) &=\forall n, n \in \mathbb{Z}:\left(\frac{1}{T} \int_{0}^{T} f\left(t-t_{0}\right) e^{-\left(i \omega_{0} n t\right)} \mathrm{d} t\right) \\
&=\forall n, n \in \mathbb{Z}:\left(\frac{1}{T} \int_{-t_{0}}^{T-t_{0}} f\left(t-t_{0}\right) e^{-\left(i \omega_{0} n\left(t-t_{0}\right)\right)} e^{-\left(i \omega_{0} n t_{0}\right)} \mathrm{d} t\right) \\
&=\forall n, n \in \mathbb{Z}:\left(\frac{1}{T} \int_{-t_{0}}^{T-t_{0}} f(\bar{t}) e^{-\left(i \omega_{0} n \bar{t}\right)} e^{-\left(i \omega_{0} n t_{0}\right)} \mathrm{d} t\right) \\
&=\forall n, n \in \mathbb{Z}:\left(e^{-\left(i \omega_{\mathrm{o}} n \bar{t}\right)} c_{n}\right)
\end{aligned}
$$
I did not understand what happened in the second step. Furthermore, in the third step why did we replace $t-t_{0}$ by $\overline{t}$ without modifying the boundary of integration and replacing $dt$ by $d\overline{t}$? Any help is much appreciated and thank you.

Even though Latex is highly encouraged instead of just copying stuff, in this case it would have been interesting to see the original "proof" because it's hard to know which errors come from the original proof and which are just typos introduced by entering the Latex code.
Of course $dt$ should have been replaced, that's just a simple typo. And the integration bounds could have been changed, but that's not necessary in this case because the integrand is $T$-periodic, and hence, as long as you integrate over one period the exact upper and lower bounds don't matter:
$$\int_{0}^Tg(t)dt=\int_{t_0}^{T+t_0}g(t)dt\qquad\text{if}\quad g(t)=g(t+T)$$


bandwidth limited pulse,what window function is it? I want to know what kind of shape a multiple cycle sine wave pulse ( or maybe the correct term is burst ) should have to have the narrowest bandwidth possible at given pulse width. In the laser science field, modelocked lasers producing ultrashort pulses are said to make transform limited pulse. Transform limited and bandwidth limited are the same thing: it's a pulse with minimum spectral width possible at the specific length of that pulse. When I was looking at the amplitude envelope shape of these pulses, it reminded me of the window function like it's used for example in fft, and it looked kind of like gaussian (the sort of gaussian that starts and ends with zero):

My question is what amplitude envelope shape should that sinewave multiple cycle pulse have to have minimum spectral bandwidth for its duration? At first I thought, that's easy, it must be gaussian (the "confined gaussian" in wikipedia, zero at edges), but when I thought about it more I am not so sure anymore. Wouldn't for example a pulse with Nuttall shape be narrower spectraly? What about Hanning and Blackman? I know all these look almost the same but hey, the true transform/bandwidth limited pulse can be only one shape, so which is it? https://en.m.wikipedia.org/wiki/Bandwidth-limited_pulse

From the great many types of window functions described on Wikipedia, the familiar rectangular window stand out has having a pretty small bandwidth. It's unfortunately also one with about the worst spectral leakage.

Now for a given maximum allowable spectral leakage level, you might want to have a look at the Dolph–Chebyshev window, which minimize the norm of the side-lobes for a given main lobe width (or equivalently minimizes the main lobe width for a given maximum side-lobe level). Note however that this window has impulse like discontinuities at the edges in the time-domain due to the constant side-lobe level in the frequency-domain. These discontinuities are more significant for larger side-lobe level as can be seen in the following graph:


Differential Equation of a Band Pass Filter I'm struggling with a signal analysis electrical engineering college assignment. We need to a derive differential equation for a low-pass filter, high-pass filter and a band pass filter (made by connecting the output of a low pass filter to the input of a high pass filter). Here are diagrams for reference:

Low pass filter

High pass filter

$$$$
Band pass filter (buffer can be removed)

The equations I derived are:

low pass filter:
      $$V_o=V_i-RC \frac{dV_i}{dt}$$
high pass filter:
      $$V_o=RC \frac{dV_i}{dt}$$
band pass filter:
      $$V_o = R_2C_2 ( (1-C_1-R_1)- R_1C_1  \frac{d^2V_i}{dt^2} )$$

The equation for the band pass filter I found by making the input of the high pass filter the output of the low pass filter. So I derived the equation of the low pass filter with respect to time and got:
$$
\frac{dV_o}{dt} = \frac{dV_i}{dt} - R_1\frac{dVi}{dt} - C_1\frac{dV_i}{dt} - R_1C_1\frac{d^2V_i}{dt^2}
$$

I substituted in this rate of change of voltage into the equation of the high pass filter to get the equation I derived for the band pass filter.

Could anyone please tell me if what I have done is valid and if not, please explain to me what the correct strategy is. Thank you kindly.

I'm afraid your results are not correct. Let's have a look at the low-pass filter. The current $i(t)$ through $R$ and $C$ must be

$$i(t)=\frac{v_i(t)-v_o(t)}{R}=C\frac{dv_o(t)}{dt}$$

So we get

$$v_i(t)=v_o(t)+RC\frac{dv_o(t)}{dt}\tag{1}$$

For the high-pass filter we have

$$i(t)=\frac{v_o(t)}{R}=C\frac{d(v_i(t)-v_o(t))}{dt}$$
which gives

$$RC\frac{dv_i(t)}{dt}=v_o(t)+RC\frac{dv_o(t)}{dt}\tag{2}$$

Let's call the time constants of the low-pass and high-pass filters $\tau_L=R_LC_L$ and $\tau_H=R_HC_H$, respectively. For the band-pass filter we need a relation between $v_i(t)$ of the low-pass filter and $v_o(t)$ of the high-pass filter. If we use $v_{oL}(t)$ to denote the output of the low-pass filter, which equals the input of the high-pass filter, we get from (1) by taking the derivative

$$\frac{dv_i(t)}{dt}=\frac{dv_{oL}(t)}{dt}+\tau_L\frac{d^2v_{oL}(t)}{dt^2}\tag{3}$$

From (2) we have

$$\tau_H\frac{dv_{oL}(t)}{dt}=v_o(t)+\tau_H\frac{dv_o(t)}{dt}\tag{4}$$

and (by taking the derivative)

$$\tau_H\frac{d^2v_{oL}(t)}{dt^2}=\frac{dv_o(t)}{dt}+\tau_H\frac{d^2v_o(t)}{dt^2}\tag{5}$$

Plugging (4) and (5) into (3) we finally get for the band-pass filter

$$\tau_H\frac{dv_i(t)}{dt}=v_o(t)+(\tau_L+\tau_H)\frac{dv_o(t)}{dt}+\tau_L\tau_H\frac{d^2v_o(t)}{dt^2}$$

By the way, if you were allowed to use the Laplace or Fourier transform, everything would be a lot easier.


Impulse vs Step response What is the difference between an impulse and a step response diagram?

If the impulse response is $h[n], n\geq 0$, then the step response $s[m]$ at any time $m≥0$ is
$$s[m]=h[0]+h[1]+h[2]+\cdots +h[m].$$
Note that if the filter is a finite impulse response (FIR) filter so that
$h[m] = 0$ for $m > M$, then the step response at $M$ is
$$s[M] = h[0]+h[1]+h[2]+\cdots +h[M]$$
and remains at this value forever afterwards, that is, $s[m]=s[M]$ for all $m \geq M$.

What other information do you want to get? The transfer function of the filter? That's easy:
$$H(z)=h[0]+h[1]z^{−1}+h[2]z^{−2}\cdots $$
which works out to be polynomial of degree $M$ in $z^{-1}$ for a FIR filter.


Determine Parameters of gaussian distribution to fit a certain Image-Filter I'd like to know what the parameters for the two-dimensional gaussian distributions are, which create this (I think) quite common image-processing filter:
$$
h_\sigma(x,y)=
\left[
\matrix
{
1&2&1\\
2&4&2\\
1&2&1
}
\right]
$$
The two-dimensional gaussian equation looks like:
$h_\sigma(x,y)= A\cdot e^{-\frac{x^2+y^2}{2\sigma^2}}$, where $A = \frac{1}{{2\pi}\sigma^2}$, but I think $A$ is not important as it is only a factor, isn't it?

With $\sigma=1$ I get:
$$h_{\sigma=1}(x,y)=
\left[
\matrix
{
h(-1,-1)=0.3679&h(0,-1)=0.6065&h(1,-1)=0.3679\\
h(-1,0)=0.6065&h(0,0)=1.0&h(1,0)=0.6065\\
h(-1,1)=0.3679&h(0,1)=0.6065&h(1,1)=0.3679
}
\right]
$$
As you see, I think the entries in the filter-matrix correspond to sampling values of the gaussian function from [-1|-1] to [1|1]. I hope that is at least right.

Multiplied $h_{\sigma=1}(x,y)$ with $4$, and then rounded gives me my desired filter, but I don't think that is the correct way, is it?

$$\mbox{round}(4\cdot h_{\sigma=1}(x,y))=
\mbox{round}(4\cdot
\left[
\matrix
{
0.3679&0.6065&0.3679\\
0.6065&1.0&0.6065\\
0.3679&0.6065&0.3679
}
\right]
)=
$$
$$
h_{\sigma=1}(x,y)=
\mbox{round}(
\left[
\matrix
{
1.4716&2.4260&1.4716\\
2.4260&4.0&2.4260\\
1.4716&2.4260&1.4716
}
\right]
)=
$$
$$
\left[
\matrix
{
1&2&1\\
2&4&2\\
1&2&1
}
\right]
=h_\sigma(x,y)
$$

There must be a better solution, or mathematical explanation for creating that matrix? I especially dislike the rounding operation, as the values are all close to $.5$ and not to $.0$ which should results in a huge error...

My question in other words: How do I determine the filter from scratch...

The filter you referenced is known as a Binomial filter. It is an approximation to a Gaussian, but for smaller filters it's a very crude one. The design is oriented more toward efficiency than accuracy. Rather than sampling a Gaussian directly, the idea behind the approximation is based on the Central Limit Theorem. In this case, it means that a small moving average filter $[ 1 \space  1 ]$ convolved with itself over and over will become more and more Gaussian. It also happens to yield the rows of Pascal's Triangle. The coefficients for larger filters are given by the following formula:

$$
{N \choose k} = \frac{N!}{(N - n)!\space N!}.
$$

For example, for $N = 5$, taking the outer product of the corresponding row of Pascal's triangle with itself yields

$$
\begin{bmatrix} 1 \\ 4 \\ 6 \\ 4 \\ 1 \end{bmatrix}\begin{bmatrix} 1 & 4 & 6 & 4 & 1 \end{bmatrix} =
\begin{bmatrix}
1 &   4 &   6  &  4  &  1 \\
4 &  16 &  24 &  16  &  4 \\
6 &  24 &  36 &  24  &  6 \\
4 &  16 &  24 &  16  &  4 \\
1 &   4 &   6 &   4  &  1
\end{bmatrix}.
$$

The paper I referenced in the first link has the following diagram which shows how the quality of the approximation increases with increasing filter order:


Hand written DFT vs MATLAB FFT has flipped conjugates? I had to compute the DFT of a signal with 5 points with a hand-written method and using MATLAB's FFT just to see if they are equal. Taking the absolute value of both gives the same spectrum, however the conjugates are in a different order for both. Why is this?

values using handwritten DFT, negative conjugates first

values using MATLAB FFT, positive conjugates first

The code I used is shown below

Now we compute the DFT of the finite sequence to show that DFT is
actually a sampled version of the DTFT.

N=5;
W=exp(-1i*2*pi/N); %the complex term W_N defined in lecture 2;
DFTmatrix = W.^((0:N-1)'*(0:N-1)); %the DFT matrix of size N
%Multiply the DFT matrix with the time-domain samples to get DFT
DFTofTruncatedSampledSignal = DFTmatrix*sampledsignal.';

Task 5:
FFTSpectrum =fft(sampledsignal,N)'; %MATLAB fft function to compute DFT

I think your problem might be in this line:

FFTSpectrum = fft(sampledsignal,N)';

Note that in MATLAB, ' is the complex conjugate transpose. Use .' to transpose a matrix without applying the complex conjugation.


How to interpret the sign of a cepstral coefficient? I feel like I've got a decent intuition for what the magnitude of a cepstral coefficient represents (i.e. how much of the signal's energy is contained in harmonics of a given wavelength), but I don't really understand what the sign of the cepstrum represents. I feel like it's important to know since some speech processing applications take the max of the real part of the cepstrum (see here), and I'm wondering if there's a reason why the max is in any way more informative than the min, or if I should instead be taking the max of the absolute value of the cepstrum?

For the zero-th cepstral coefficient, the sign is probably always positive since it represents the average energy of the entire spectrum.

For higher-order coefficients, seeing that the question was first posted in 2017, does the OP maybe already have an answer? Would be interested to learn.


Kalman filtering in image processing, resources? I'm looking for a good resource (book, tutorial, lesson etc.) that explains the usage of Kalman filtering in image processing applications.

I'm aware of the fact that Kalman filtering is an optimal estimator which is a tool not a concept in image processing because of that I didn't find a book that explains the usage of Kalman filtering for real applications (from a technical perspective).

I use OpenCV (Matlab will be fine), so I prefer resources that handle image processing in OpenCV.

Here are a couple of examples from the computer vision system toolbox:
Using kalman filter for object tracking

For a more in-depth explanation the best book is Multiple Target Tracking with Radar Applications by Samuel Blackman.


How can I "smooth" square wave/pulse data to remove/reduce harmonics? First, I have reviewed close-match questions, but can't find something that seems to match what I'm looking for.  Also, it's been 30 years since college, so forgive me for forgetting stuff.. :)

I am analyzing IP Packet data.  I am creating bins (of varying window sizes.. The window seems not to be sufficient to remove/smooth the harmonics) to count the frequency of certain attributes/behaviors within the data.  From this, I am trying to determine fundamental frequencies that activities are occurring.  I believe that the artifacts that I am seeing are harmonics since the data tends to be pulse/bursty, resulting in what amounts to a square wave.  (see below) . In this sample data, as you can see, there is a strong signal occurring every hour, over the course of 47 hours, yet it is almost impossible to "see" because of the harmonics. (at least, that's my guess/intuition)

Is there something that I can do to smooth things out so that the harmonics are suppressed (or largely suppressed) and I only see the fundamental frequencies?  I feel as though I've forgotten something basic.

If you are interested in finding the periods, have a look at autocorrelation lag.
Maybe try some xcorr options too.
MATLAB example code:

nh = 47;  % 47 hours.
% Pulse period of 1 hr and make sequence length not a multiple of period
n = 3600*nh + 273; % 1 sec sampling rate
x = zeros(1,n);
% Create a pulse train with varying heights and some noise
x(777:3600:end) = 8 + 2*rand(1,nh);
x = x + rand(1,n);
% Look at lags up to 4 hours
xcorr(x, 4*3600, 'unbiased');

Zero lag is in the middle of the figure.

Note that the consistent heights show reasonably consistent correlation at integer multiples of the base period.
If you want to tolerate timing jitter in the period, you can filter the data with a lowpass FIR filter first, which will smear the peaks in the data and, of course, smear the peaks in the autocorrelation., e.g.,

b = hamming(15);  % Smooth, simple, positive, symmetric, integer delay
x = filter(b,1,x);


Amplitude transfer function Given a LTI-system  $$y[n]=x[n]-2x[n-1]+y[n-1]- \frac{8}{9}y[n-2]$$
The transfer function $H(z)$ is:
$$H(z) = \frac{1-2z^{-1}}{1-z^{-1}+ \frac{8}{9} z^{-2}} $$
How do I calculate the amplitude transfer function so I can plot a Bode diagram like this? (in hand, not in MATLAB)

I thought that I should substitue $z=e^{j\omega}$ in $H(z)$, but it doesn't seem to be right. I hope you guys can help me out.

There are many algebraic paths to get into the required expression some of them being long and tedious, while the others being shorter with a carefully choosen route.

Given :
$$H(z) = \frac{ 1 - 2z^{-1} }{1 - z^{-1} + \frac{8}{9}z^{-2} }$$

To evaluate the corresponding Bode-Plot, by hand, we need to first derive the expression called as the Frequency Spectrum Magnitude $\lvert H(e^{j\omega})\rvert$ of $H(z)$, provided that the system is stable (i.e.all poles of $H(z)$ inside unit circle) And then compute:

$$20\log_{10}\left(\lvert H(e^{j\omega})\rvert\right),$$
to get the Magnitude in dB scale.

Let's proceed into the first step and plug $e^{j\omega}$ into $H(z)$ to get Frequency Response:
$$H(e^{j\omega}) = \frac{ 1 - 2e^{-j\omega} }{1 - e^{-j\omega} + \frac{8}{9}e^{-2j\omega} }$$

As this is a complex valued rational fraction expression, it is equal to:
$$H(e^{j\omega}) = \frac{P(e^{j\omega})}{Q(e^{j\omega})}$$

And from complex algebra, the magnitude of $H(e^{j\omega})$ is equivalent to:
$$ \lvert H(e^{j\omega})\rvert = \frac{\lvert P(e^{j\omega})\rvert}{\lvert Q(e^{j\omega})\rvert}$$

Hence we need to compute the absolute values of the numerator and denumerator expressions, where the absolute value of a complex number is understood as: $$z=a+jb \implies \lvert z\rvert=(a^2+b^2)^{0.5}$$

Lets elaborate the Numerator:
\begin{align}
P(e^{j\omega}) &= 1 - 2e^{-j\omega} = \left(1-2\cos(\omega)\right) + j\left(2\sin(\omega)\right)\\
\lvert P(e^{j\omega})\rvert &= \left[\left(1-2\cos(\omega)\right)^2 + \left(2\sin(\omega)\right)^2\right]^{0.5}\\
\lvert P(e^{j\omega})\rvert &= \left[1+4\cos^2(\omega)-4\cos(\omega) + 4\sin^2(\omega)\right]^{0.5}\\
\lvert P(e^{j\omega})\rvert &= \left[5-4\cos(\omega)\right]^{0.5}\\
\end{align}

Also elaborate the Denumerator (Denominator I mean!):
\begin{align}
Q(e^{j\omega}) &= 1 - e^{-j\omega} + \frac 89 e^{-2j\omega} = e^{-j\omega} \left(e^{j\omega} - 1 + \frac 89 e^{-j\omega}\right)\\
\lvert Q(e^{j\omega})\rvert&=\lvert e^{-j\omega}\rvert\cdot \bigg\lvert \left(e^{j\omega} - 1 + \frac 89 e^{-j\omega}\right)\bigg\rvert\\
\lvert Q(e^{j\omega})\rvert&=\bigg\lvert \left(e^{j\omega} - 1 + \frac 89 e^{-j\omega}\right)\bigg\rvert\\
\lvert Q(e^{j\omega})\rvert&=\left[\left(\cos(\omega) - 1 + \frac 89 \cos(\omega)\right)^2 + \left(\sin(\omega) - \frac 89\sin(\omega)\right)^2\right]^{0.5}\\
\lvert Q(e^{j\omega})\rvert&=\left[\left(\frac{17}{9}\cos(\omega) - 1\right)^2 + \left(\frac 19 \sin(\omega)\right)^2\right]^{0.5}\\
\end{align}

As there seems no more useful simplifications, we can deduce $\lvert H(e^{j\omega})\rvert$ as:
$$\lvert H(e^{j\omega})\lvert = \frac{\left[5-4\cos(\omega)\right]^{0.5}}{\left[\left(\frac{17}{9}\cos(\omega) - 1\right)^2 + \left(\frac 19\sin(\omega)\right)^2\right]^{0.5}}$$

The Bode plot takes the logarithm of this to produce the graph:
$$\begin{align} \lvert H(e^{j\omega})\rvert_\text{dB} & \triangleq 20 \log_{10}\left(\lvert H(e^{j\omega})\rvert\right) \\
& = 20 \log_{10}\left(\frac{\left[5-4\cos(\omega)\right]^{0.5}}{\left[\left(\frac{17}{9}\cos(\omega) - 1\right)^2 + \left(\frac 19\sin(\omega)\right)^2\right]^{0.5}}\right)\\
& = 10 \log_{10}\left(\frac{5-4\cos(\omega)}{\left(\frac{17}{9}\cos(\omega) - 1\right)^2 + \left(\frac 19\sin(\omega)\right)^2}\right) \\
& = 10 \log_{10}\left(5-4\cos(\omega)\right) - 10 \log_{10}\left(\left(\tfrac{17}{9}\cos(\omega) - 1\right)^2 + \left(\tfrac 19\sin(\omega)\right)^2\right)\\
\end{align}$$

Which produces the same graph when compared with freqz([1 -2],[1 -1 8/9]) function of MATLAB to produce discrete time frequency response plot, for a linear frequency axis. In discrete time it is customary to use linear frequency axis rather than a logarithmic one as usual in continuous time Bode-Plots.


When can the impulse response become zero? The article Efficient Use Of Sparse Adaptive Filters (Proc. Asilomar Conference, Khong et al., 2006) introduces adaptive filters for the estimation of channels or systems having a sparse impulse response. Assuming that the channel is modeled by an FIR model, then the output of the channel will be given
by $$y[n] = \sum_{i=0}^{p-1} h[i] x[n-i] \tag{1}$$ where $x$ is the input to the model and $h$ is the MA coefficients. As an example, consider an MA(3) model of order p=3 having coefficients h = [1,0.2,0.5,0.6]

I can expand the equation to write as $y[n] = x[n] + 0.2*x[n-1] + 0.5*x[n-2] + 0.6*x[n-3]$

My problems are:

From a time series or statistics point of view, FIR filter = Moving Average and IIR = Autoregressive model. I cannot understand how to relate the delays, 1,2,3 denoted by the order to the lags terminology often used in statistics. The impulse are known as the coefficients of the model and how to relate the impulse to the coefficients?
Why or rather what is the intuitive meaning of the coefficients becoming zero?
When does that happen/some real life examples of applications?

The notion of sparsity entails that an object, living for instance in an $n$-dimensional space, can be described (in the suitable basis/frame) by a number $k$ of meaningful components (each above a threshold, or whose combination is close enough to the signal)  "much smaller" than $n$.

When talking about filter identification (adaptive or not), the filter is, at the beginning, unknown. To avoid looking inside an infinity of possibilities, one first restricts  the space of solutions to the subspace of length-$n$ or order-$n$ filters. Either for true physical reasons, out of wishful thinking or heuristic motivations, or because  the optimisation problem is not directly tractable (non-smooth, non-convex), one may call for Ockham's razor, or the law of parsimony, one may want to restrict the solution (not the real object) to be $k$-sparse.

Your example could  be sparse per se, if initially the space of solutions could have up to 100 taps. Additionally, it could perhaps be approximated (in practice) by a sparser filter:
$$y[n] = x[n] + 0.6*x[n-3]$$
discarding the two smallest coefficients.

With this introduction, let us go back to your questions:

I would disagree with "FIR filter = Moving Average and IIR = Autoregressive model": an IIR can be AR/MA, some FIR can be rewritten in an AR fashion. The rest of the question is unclear to me: if the model is linear, the impulse response would be some sparse approximation of the actual coefficients.
The actual coefficients are not becoming zero. But, either they form a  filter short enough to be qualified as sparse, or their estimate has very few non-zero taps
It may happen in some echo cancellation problems: either in a room or on delay lines, one expect that the source bounces after some travelling back and forth. If at each bouncing the reflection is scarce, and there is some attenuation, the  "global" echo filter could look like  (two bounces)
$$y[n] = x[n] - 0.25*x[n-99]+ 0.5*x[n-100]+ 0.25*x[n-101]  + 0.125*x[n-199]+ 0.25*x[n-200]+ 0.125*x[n-201]$$
hence a sparse filter, even with order 201. Outside speech processing, we used such priors in seismic processing of remove multiples, approximate models of bouncing waves, assuming that the filters were sparse at boundaries between rock layers, see A primal-dual proximal algorithm for sparse template-based adaptive filtering: Application to seismic multiple removal, 2014, IEEE Transactions on Signal Processing:

  Unveiling meaningful geophysical information from seismic data
  requires to deal with both random and structured “noises”. As their
  amplitude may be greater than signals of interest (primaries),
  additional prior information is especially important in performing
  efficient signal separation. We address here the problem of multiple
  reflections, caused by wave-field bouncing between layers. Since only
  approximate models of these phenomena are available, we propose a
  flexible framework for time-varying adaptive filtering of seismic
  signals, using sparse representations, based on inaccurate templates.
  We recast the joint estimation of adaptive filters and primaries in a
  new convex variational formulation. This approach allows us to
  incorporate plausible knowledge about noise statistics, data sparsity
  and slow filter variation in parsimony-promoting wavelet frames. The
  designed primal-dual algorithm solves a constrained minimization
  problem that alleviates standard regularization issues in finding
  hyperparameters. The approach demonstrates significantly good
  performance in low signal-to-noise ratio conditions, both for
  simulated and real field seismic data.

Note that parsimony, usually resorting to an $\ell_0$ count measure, is not easy to optimize, as non-differentiable, non-convex, etc.


MFCC classification model I have audio samples which MFCCs i want to train, but there is a problem. I can't find a classification model, because the samples have different length and consequently the MFCC matrices will also have different lengths. My question is - which machine learning model can i use?

When utilizing MFCCs, it is common to split the MFCC matrix into pieces. Lets assume your MFCC-matrix $\bf{M}$ has the dimensions $[\bf{M}]$$=N \times 39$, where $N$ denotes your time index and you have 39 MFCC bins. Then you cut your matrix into blocks of length $B$ (I propose $B=39$ so your blocks have a square shape) and save them in a list. For the last block, you can either zero-pad until it has the same length or discard it completely. Each block has to be annotated with the label of the original audio sample.

A standard approach for classification is a Deep Neural Network, that consists of Convolutional and adjacent Pooling layer since MFCCs are like images locally connected/correlated and Convolutional layers take care of that. Of course, you should not forget to add an MLP layer at the end with a softmax activation to get an estimated probability distribution as an output of your net as well as an appropriate regularization technique.

A more advanced approach would be the use of a recurrent network structure like an LSTM-gate or a plain Recurrent layer. Therefor, you have to store the blocks of a sample in chronological order so the network can learn the temporal structure of a whole sample.

EDIT: Actually one often takes only 13 MFCC bins. The $N \times 39$ matrix also contains the Deltas and DeltaDeltas of the MFCC bins.


why transmitting a long pulse its spectrum gets spread? I'm transmitting an RF signal with a CC1101 at 440 MHz.

This signal is a constant of 1, so the spectrum should be a delta.

Using a SDR dongle, I've checked the spectrum and it's wider than a delta but also, after a while, this peak kind of turns in many peaks that move (it looks like sea waves travelling).

Any idea why this happens?

this is John BG

1.- GFSK is not same as single tone

According to the CC1101 specs you refer to, the C1101 uses GFSK.

GFSK varies the varrier frequency, therefore the spectrum cannot be a delta.

Call it a 'shifting delta' but it cannot be same spectrum as just a tone.

2.- you want the carrier frequency to vary

ISM are crowded bands, the carrier has better chances to deliver if varying than if constant,
against jamming, interferrers, and any power level that neighbouring transmitters may radiate.

3.- not even ASK

You can still set modulation to MSK, OOK, or even ASK but even with ASK

Even with PKT_FORMAT[1:0] to Asynchronous mode, CRC disabled, Manchester disabled and PKTCTRL0 setting infinite packet length, PKTLEN is going to cut the packets to ensure syncing takes place.

This is done to avoid persistent propagation of errors along same packet that might generate uncertainty
to possibly overlap next packets.

4.- Temperature drifts, manufacturing tolerances

Such small and really low power devices may be subject to temperature drifts and manufacturing tolerances.

The packet header is needed at a high enough repetition to avoid transmitter and receiver missing frequency lock.

Hope it helps, John BG


Design Filter Banks for every 100 Hz I need to study a wave file for every 100 Hz. However, I dont know how can I design filer for each 100 Hz. Let's say 0-100 Hz, 100-200 Hz, ... , ((fs/2)-100) _  (fs/2) Hz.

Is it any MATLAB code for it ?
How many samples are needed for bandwidth of 100  Hz ? I maen how should calculate it ?
Please guide me.
Thanks.

A windowed FFT with a length that corresponds to 2/100ths of a second will give you overlapping filters each with roughly 100 Hz bandwidth.  Just use every other FFT result bin if you want 100 Hz filter spacing.  The window shape (Von Hann, Nutall, etc.) will allow some tuning of the filter shape, but the main passband lobe of a windowed FFT will be about twice the FFT bin spacing in width.

An unwindowed FFT creates filters with too much stop band ripple to usually be considered useful.

You can also zero-pad and/or use an FFT a higher integer multiple of 0.01S in length to get better control over the windowing (a proportionately wider or narrower Gaussian window, for instance) and the window's frequency response.

The number of samples needed will be the FFT length in samples, or the sample rate multiplied by the FFT's length in time.


How is voice sound on TV or phone distorted through a face mask? After a year of listening to people talk through face masks, it appeared to me that the sound is not merely muffled, but it has a very specific and recognizable tinge and texture to it.
It reminds me less of sound through a muffling pillow, and more of old speech recordings from the 1950s.
This is apparent whether you listen, today, to recorded speech on TV or radio or through the phone. Sometimes you hear it in person too.
I suspect that the sound is not just muffled or low-pass filtered, but that there is more happening. Perhaps it's a band-pass? A notch? Are there non-linear or time-variant effects from turbulence due to the interaction between plosives and the fabric?
Has anyone looked into this?

There are many researches on the acoustic effect of masks since the pandemic. You can find them on the website of The Journal of the Acoustical Society of America and most of them are free to access.
Masks have selectivity on both frequency and spatial.
In frequency domain, masks react as low-pass filters that attenuate high frequencies significantly. Interestingly, different research gives different results. They are saying 1 kHz$^{[1]}$, 2 kHz$^{[2]}$ and 3 kHz$^{[3]}$ maybe due to different understanding about "significantly".
In spatial domain, masks lead to reduced frontal sound radiation compared to the lateral or rear radiation at high frequencies$^{[2]}$.
Researchers also tested the effect of microphone placement. They placed microphones on the lapel, cheek, forehead, and next to the mouth$^{[1]}$. These differences are mostly caused by the spatial radiation pattern with a mask.
And of course, different masks have different acoustic effect according to their materials, weaves and types.
References

[1] Acoustic effects of medical, cloth, and transparent face masks on speech signals
[2] Impact of face masks on voice radiation
[3] Effects of face masks on acoustic analysis and speech perception: Implications for peri-pandemic protocols
[4] Acoustic effects of non-transparent and transparent face coverings

[5] Effect of masks on speech intelligibility in auralized classrooms


FPGA vs DSP for MAC operations I'm new to SP Stack Exchange . I am doing a project on implementing a handwriting recognition machine using Neural Networks in Real Time .

The Image Processing Part involves convolving a 4*5 kernel over a 28*28 image. I want to know the pros and cons involved with using a DSP vs FPGA .

The FPGA I plan to use belongs to Xilinx's Spartan 6 family which has 90 DSP slices. Each DSP Slice has 18*18 bit MAC units which can operate in parallel at 200MHz.

My college has a DSP starter kit TMS320C6713 DSK from Spectrum DIGITAL. I'm not sure how MAC operations are performed in DSP .Do we have parallel MAC units? The datasheet isn't very clear on the MAC units .

Are there any advantages on using the DSP over the FPGA which I might have overlooked?

Any help is appreciated .

At a high level, it's a relatively simple tradeoff: software development for DSP or general-purpose microprocessors is typically much less labor-intensive than FPGA development. At the same time, the amount of computational throughput offered by an FPGA typically far exceeds what a CPU can do, because you can tailor the logic specifically for your application.

Most any processor that is advertised for "DSP" will have multiply-accumulate (MAC) instructions available, and it is common for DSP processors to be able to execute multiple instructions per cycle. If you read the specifications for TI's C6713, you'll see that they advertise that it can complete up to 8 instructions per cycle (not necessarily all MACs, but a mix of several types of instructions). Note that it is pretty difficult to get close to that throughput, but it is possible if you invest enough effort learning how the architecture works and structuring your code accordingly.

One other major difference between an FPGA implementation and the particular DSP that you noted is that the C6713 supports native floating-point capability (instead of just fixed-point). For some applications, this can make algorithm implementation much simpler, as you don't have to worry about the issues that come along with fixed-point computations. While it is possible to implement floating-point computations in an FPGA, such an implementation would be much more demanding on the hardware, so depending on how complex your algorithm is, you might run into issues with resource utilization, achievable clock speed, and so on.

You should start by quantifying what kind of computational demands your algorithm is going to place on your hardware. Then, decide how much engineering time you can allocate to the project; from there, the choice will likely be clear.


Magnitude and phase of frequency response from Z-transform I have been given a signal which I have completed the Z-transform for and have got:
$H(z)=1-z^{-1}$. I need to find the magnitude and phase of this transform but am unsure on how to actually do this. All I have got so far is $H(e^{j\omega}) = 1 - e^{-j\omega}$.

Its frequency response can be decomposed into a real part and an imaginary part as
$$
H(e^{j\omega}) = 1-e^{-j\omega} = 1-(\cos\omega - j\sin\omega) = (1-\cos\omega) + j\sin\omega
$$
The magnitude response is the modulus of the complex frequency response
$$
|H(e^{j\omega})| = \sqrt{\mathcal{Re}\{H(e^{j\omega})\}^2 + \mathcal{Im}\{H(e^{j\omega})\}^2} = \sqrt{(1-\cos\omega)^2 + \sin^2\omega} = \sqrt{2 - 2\cos\omega}
$$
where $\mathcal{Re}$ and $\mathcal{Im}$ represent real part and imaginary part respectively.
The phase response is the phase angle of the complex frequency response
$$
\angle H(e^{j\omega}) = \arctan \frac{\mathcal{Im}\{H(e^{j\omega})\}}{\mathcal{Re}\{H(e^{j\omega})\}} = \arctan \frac{\sin\omega}{1-\cos\omega}
$$


FFT calculation suitable, Wav File, Java If I have a wav file with this characteristics:

16bits (2 Bytes per Sample)
32000 Sampling Rate
60 Seconds
2 Channels

then I have:

Samples Per Channel: 32000*60 = 1920000
Samples Per File: 32000*60*2 = 3840000
Bytes Per File: 32000*60*2*2 = 7680000

It's a huge file!

If I change my file Characteristics to:

8bits (1 Byte per Sample)
8000 Sampling Rate
60 Seconds
1 Channel

then I have:

Samples Per Channel: 8000*60 = 480000
Samples Per File: 8000*60*1 = 480000
Bytes Per File: 8000*60*1*1 = 480000

If my FFT 2^12 = 4096 or 2^16 = 65536

Do I need to calculate the FFT 480000 times?

Does each how many samples I need to calculate the FFT?  (480000 - 65536 = 414464?)

what is the suitable calculus?

How to do to show the spectrum?

To show only one spectrum then FFT of the whole sample could be done.But for better result you could take a window and take FFT. For taking windows of 4096 with 50% overlap then it is equivalent to doing FFT
no of times
$$
\frac{480000}{4096*0.5}
$$
times that is equivalent to
$$
\frac{Total\hspace{1mm} sample\hspace{1mm} length}{(window \hspace{1mm}length)*(amount\hspace{1mm} of\hspace{1mm} overlap)}
$$


Calculating the magnitude spectrum and phase spectrum From a window function $x(t)=u(t+2)-u(t-2)$, we can get the Fourier Transform $X(j\omega)=\frac{2\sin(2\omega)}{\omega}$.

Then, I want to calculate its magnitude spectrum and phase spectrum.

The magnitude spectrum is the magnitude distribution at every $\omega$, so it's simply absolute value of $X(j\omega)$, is this correct?

And how to calculate the phase spectrum?

HINT:

Magnitude $M(\omega)\ge 0$ and phase $\phi(\omega)$ are defined by

$$X(\omega)=M(\omega)e^{j\phi(\omega)}\tag{1}$$

Note that $(1)$ is generally complex-valued. In your example, the Fourier transform $X(\omega)$ is clearly real-valued. This restricts the possible values of the phase $\phi(\omega)$. What are those two values of $\phi(\omega)$ for which $X(\omega)$ in $(1)$ is real-valued? How do they relate to the sign of $X(\omega)$? I'm sure you can take it from here.


How to Plot Summation of Shifted impulses in MATLAB? I am trying to plot the following equation in MATLAB:-

I am not sure how to write the following summation equation in MATLAB.

According to the features of the unit impulse function $\delta(t)$; at t = 0, the function goes toward infinity. Actually, the unit impulse function is meaningful with its area which is 1.
Additionally, MATLAB don’t show the famous arrow notation of the unit impulse function. Instead, it shows an invisible line which can be noticed where the function increases into infinity, i.e. the time value on the time axis where the function increases into infinity is blank. This can be seen from the figure.
On the other hand, the unit sample function $\delta[n]$ value is defined with the area of the unit impulse function which is 1.
We can obtain h[n] function by sampling h(t) function. But, what about those inifinities?
So, let’s do some gamble!
If we code h(t) function in MATLAB, it will actually be a row vector. In this row vector, we can change those infinities with ones in order to sample the function correctly. This process will require a for loop and an if statement.
The MATLAB source code for the complete process is given below:
t = -10: 1: 10; % Time range

Ts = 1; % Sampling period

syms k; % Symbol assingment for summation

function_1 = (t + 1) .* symsum(dirac(t - k), k, 1, 5); % Summation

subplot(2, 1, 1);
plot(t, function_1, "r", "LineWidth", 3);
xlabel("Time (s)");
ylabel("h(t)");
title("h(t) Time Function");

for a = 1: length(function_1) - 1 % For loop-if statement duo for replacing infinity results with unit sample function
    if function_1(a) == Inf
       function_1(a) = 1;
    end
end

function_2 = (t + 1) .* function_1; % Complete function to be sampled

n = t / Ts; % Sample range

subplot(2, 1, 2);
stem(n, function_2, "ro", "LineWidth", 3);
xlabel("Samples [n]");
ylabel("h[n]");
title("h[n] Sample Function");

According to the figure, sampling process has gone well and it matches with analysis that I’ve done on the paper.

Actually, there may be other approaches for the replacements of the infinities with ones but this one does the trick well.


Calculating SNR of a complex-valued signal correctly I have an array of complex-valued data. I would like to modify the array by injecting an artificial signal at a chosen Signal to Noise ratio.
I have seen the signal to noise ratio for complex-valued data calculated as:
$$ SNR = \frac{<SS^*>}{<NN^*>} $$
Right now, I am assuming my data is all noise and calculating the denominator by taking the mean of the real part of the data multiplied by its imaginary part.
One problem I have with this equation is that it would imply that if I added a constant real signal to all the data (for example just adding $ 0.5 + 0j $) the SNR of the added signal would be zero since its imaginary part is zero and hence $ SS^* $ is zero. Obviously this doesn't make sense. So how do I correctly calculate the SNR of an added signal? I'm looking for a way of determining $a$ and $b$ in a complex number $ a+ib $ added to all the data, such that the resulting data has an artificial signal of $SNR = X$ embedded in it.

Just simply compute the variance of the signal and the variance of the noise. (And adjust the noise for sampling factor if necessary, see below). Variance is not a complex result and the  complex conjugate is used in the computation so should be invariant to real or imaginary content; 5+j0 and 0+j5 have the same variance; if the entire signal rotated in phase, how would that have any effect on the power?
Note computing SNR when you are creating a signal is trivial as you can completely separate the signal and noise components, scale appropriately and add. Measuring SNR on a combined signal is a little more challenging and for this purpose I recommend using correlation to measure SNR given the direct relationship between rho (correlation coefficient) and SNR.
Also pay attention to the spectral density of the noise to set the SNR properly such that it is scaled to be just the portion of the noise that is in the bandwidth of the signal.


How to solve this question by using convolution? I had previously posted this question but my method was wrong in it.
Two signals are given:

$$x[n]=2^nu[n]$$
$$h[n]=u[n]$$

I have to find $y[n]$. I have spent hours on this question but I am unable to solve it. Kindly please tell me the way to solve it.

Assuming you mean that $y[n] = h[n] * x[n]$, then
\begin{equation}
\begin{split}
 y[n]  &\stackrel{(a)}{=} h[n] * x[n]\\
 &\stackrel{(b)}{=} \sum_{m=-\infty}^{\infty}
 h[m]x[n-m]\\
 &\stackrel{(c)}{=} \sum_{m=-\infty}^{\infty}
 u[m]2^{n-m}u[n-m]\\
 &\stackrel{(d)}{=}2^n
 \sum_{m=-\infty}^{\infty}
 u[m]2^{-m}u[n-m]\\
 &\stackrel{(e)}{=}2^n
 \sum_{m=0}^{\infty}
 2^{-m}u[n-m]\\
 &\stackrel{(f)}{=}2^n
 \sum_{m=0}^{n}
 2^{-m}\\
 &\stackrel{(g)}{=}2^n
 \frac{1 - (\frac{1}{2})^{n+1}}{1 - \frac{1}{2}}\\
 &\stackrel{(h)}{=}2^{n+1} - 1
\end{split}
\end{equation}

(a) is the convolution symbol.

(b) is using the definition of discrete convolution.

(c) replacing the quantities.

(d) $2^n$ is independent of the summation index.

(e) $u[m]=1$ if $m\geq 0$ and $0$ otherwise.

(f) if $n-m<0$, or if $n<m$ then the summation is zero. So we have (for non-zero quantities), $m = 0 \ldots n$.

(g) Using geometric series

(h) rearrangements.


Data representations for speech data I'm interesting in knowing more about different kinds of data representations that are used for speech data. I know of FFT, but not any others. More precisely, I have the 3 following questions.

What data representations for speech data?
What is considered the most 'raw' representation?
How has the use of different representations evolved over time?

Simply answering 1 or 2 of the questions would also be great. Thanks.

I have seen people use MFCC coefficients for speech related analysis.It is a  done using filter banks and logarithmic scales for representing speech.


How does sampling rate impact Discrete-Time Kalman Filter state space modeling assumptions? Consider a very simple, discrete-time constant position-type model for state updating in a Kalman filter:

$$
x_{k+1} = x_k + w_k
$$

The Kalman filter will be run with update interval $T_s$ such that $x_k$ represents position at time $t=kT_s$.

Typically, the process noise $w_k$ is modeled as a zero-mean, white stationary random process of variance $\sigma^2_w$.

Now, assume that the position time series is strictly band-limited with bandwidth far lower than the Nyquist frequency $1/(2T_s)$.  This corresponds to a high degree of oversampling.

In principle (say, with zero measurement noise), the process noise could be recovered by applying a simple differencing high pass filter with transfer function:

$$
H(z) = 1 - z^{-1}
$$

to the position time series.  However, if the position time series is strictly bandlimited, then the recovered process noise is strictly bandlimited, which violates the white process noise model.

Does this mean that the Kalman filter needs to use a correlated process noise model for highly oversampled systems?

In my case, the state space model (which, in fact, is constant velocity not constant position) describes the time evolution of a certain biological quantity.  The discrete-time state space model has not been chosen from a precise description of the dynamics of the underlying physical (physiological) phenomena,  (This is because such descriptions are complex and require access to additional biological quantities which are rarely knowable in practice).  Rather, the constant velocity model is a common choice in the literature motivated primarily by expedience, simplicity, and good tracking results.

In the literature, all results seem to be reported for near-Nyqvist sampling. Some studies in the literature analyzed real biological data (again, sampled near Nyqvist) and claimed that the process noise for this model was white.  In my case, for historical reasons, I am oversampling by a good factor of 15.  When I analyzed a very clean sample signal, I found that at the oversampled rate, the process noise was highly correlated.  A decimated-to-near-Nyqvist version of the signal gave rise to a significantly lower amount of temporal correlation in the process noise.  This is what motivated my question.

The answer to your question is : maybe. When one is at the discrete time model, it is generally assumed that you solved the state transition equations from the continuous time differential equations where the physics have already been appropriately modeled, which is where your answer lies. Time discretization is based on sampling the continuous time state transition matrix, thus the sampling intervals aren't related to Nyquist, they correspond to the fidelity of the modeling. Shorter intervals tend to be more accurate approximations. The time intervals don't have to be uniform, but uniform intervals provide a less complicated implementation.

While shorter time intervals correspond to a heuristic where shorter has better approximation, you raise a valid concern, sampling too closely, may indeed introduce correlation in the process noise.  This can be solved with an augmented state.

These are all considerations. Kalman Filters need to be tuned.  The measurement noise has been neglected in this discussion.

So, the answer to your question is that it really depends.


Wiener Filter for Image Noise Reduction (Image Denoising) I'm trying to get my head round the operation of the Wiener filter for the purpose of image noise reduction. In my case I'll have used another noise reduction filter first and will then use the result of this as an approximation of the noise characteristics for the Wiener filter.

Regarding information on the Wiener filter, I found the following Matlab code & description useful:

http://www.mathworks.co.uk/help/toolbox/images/f11-12251.html#f11-14272

and a few other good links such as

http://blogs.mathworks.com/steve/2007/11/02/image-deblurring-wiener-filter/

So from a Matlab perspective I can see how to use the inbuilt Matlab function, but I'd like to gain a more fundamental understanding rather than just use the function call, yet at the same time I'd prefer to find something more digestible than the Wikipedia entry on Wiener filtering.

Anyone care to offer a brief explanation on Wiener filtering?

There is another Wikipedia entry on Wiener filtering more applicable to image processing.

To summarize (and convert to 2D), given a system:
$$
y(n,m) = h(n,m) * x(n,m)  + v(n,m)
$$
where

$*$ denotes convolution,
$x$ is the (unknown) true image,
$h$ is the impulse response of a linear, time-invariant filter,
$v$ is additive unknown noise independent of $x$, and
$y$ is the observed image.

We want to find a deconvolution filter $g$ so that we can estimate $x$ as follows:
$$
\hat{x}(n,m) = g(n,m) * y(n,m)
$$
where $\hat{x}$ is an estimate of $x$ that minimizes the mean square error.

In the frequency domain, the transfer function of $g$, $G$ is:
$$
G(\omega_1, \omega_2) = \frac{H^*(\omega_1, \omega_2)S(\omega_1, \omega_2)}{|H(\omega_1, \omega_2)|^2S(\omega_1, \omega_2) + N(\omega_1, \omega_2)}
$$
where

$G$ is the Fourier transform of $g$,
$H$ is the Fourier transform of $h$,
$S$ is the mean power spectral density of the input $x$, and
$N$ is the mean power spectral density of the noise $v$.

The equation for $G$ can be re-written as:
$$
G(\omega_1, \omega_2) = \frac{1}{H(\omega_1, \omega_2)}
\left[ \frac{|H(\omega_1, \omega_2)|^2}{H(\omega_1, \omega_2)|^2 + \frac{N(\omega_1, \omega_2)}{S(\omega_1, \omega_2)}}
\right]
$$
So the Wiener filter has the inverse filter for $H$, but also a frequency-dependent term that attenuates the gain based on the signal to noise ratio.


Laplace transform of a finite duration signal Consider the following signal:
$$ x(t) = e^{-2t}[u(t) - u(t-5)] $$

This signal exists only from 0 to 5 time units. Elsewhere, it is zero.

Now, let's find the laplace transform of this signal using Linearity and Time shift properties.

$$ e^{-2t}u(t) \leftrightarrow \frac{1}{s+2} \ , \ \ Re \{s \} > -2 $$
Also,
$$ e^{-2(t-5)}u(t-5) \leftrightarrow \frac{e^{-5s}}{s+2} \ , \ \ Re \{s \} > -2 $$
$$ \Rightarrow e^{-2t}u(t-5) \leftrightarrow \frac{e^{-5s}e^{-10}}{s+2} \ , \ \ Re \{s \} > -2 $$

Thus, by linearity property,
$$ e^{-2t}[u(t) - u(t-5)] = \frac{1 - e^{-5(s+2)}}{s+2} \ , \ \ Re \{s \} > -2$$

Note: The time shifting property doesn't alter the ROC;

However, the textbooks that i am refering (Oppenheim and Schaum series) both tell that the ROC of a finite duration signal is the entire S-plane, possibly zero and infinity (in some cases).

But the above signal being of finite-duration, possess ROC that is not the entire s-plane. Please help me figure this conceptual error.

Note: The above problem is from Schaum series. Here are the images of the textbook's section relevant to the above question.

Source of the Question and its solution:

Property of finite duration signals:

In Schaum's outline series:

In oppenheim:

The property claimed by Schaum and Oppenheim is also true for the given example. Note that the Laplace transform

$$X(s)=\frac{1-e^{-5(s+2)}}{s+2}\tag{1}$$

has no pole at $s=-2$:

$$\lim_{s\to -2}X(s)=\frac{1-(1-5(s+2))}{s+2}\Big{|}_{s=-2}=5$$

So the ROC is indeed the entire $s$-plane. Even though the ROCs of the two individual signals in your solution have the same right half-plane as ROC, their sum has the entire $s$-plane as ROC because the two signals cancel everywhere except in a finite interval.


retrieving original data from phase and magnitude of Fourier transform I use this snippet of python code to transform data to Fourier phase and magnitude and then retrieving original data.

from random import randint as RI
import numpy.fft as FFT
import math

w = 4
h = 4
random_range = 255

vals = [[] for i in range(h)]

for i in range(h):
    for j in range(w):
        vals[i].append(RI(0,random_range))

fftc = FFT.fft2(vals)

magnitudes = [[] for i in range(h)]
phases = [[] for i in range(h)]

for i in range(h):
    for j in range(w):
        magnitudes[i].append(math.sqrt(fftc[i][j].real**2+fftc[i][j].imag**2))
        phases[i].append(math.atan(fftc[i][j].imag/fftc[i][j].real))

for_ifft = [[] for i in range(h)]

for i in range(h):
    for j in range(w):
        rl = magnitudes[i][j]*math.cos(phases[i][j])
        im = magnitudes[i][j]*math.sin(phases[i][j])
        for_ifft[i].append(complex(rl,im))

ifftc = FFT.ifft2(for_ifft)

print vals
print ifftc

When i compare outputs, results are different. for example two typical outputs are:

//original data
[[115, 4, 33, 91], [228, 123, 46, 111], [83, 227, 81, 55], [229, 188, 138, 41]]

//retrieved data from phase and magnitude informations
[[ 309,182.25,105,57.25],[ 33,118.75,120.50,116.25],[ 148,128.75,79,94.75],[ 165,112.25,-6.50,29.75]]

I know the theory behind it. for example for calculating phase and magnitude we do this :

Magnitude = Square_Root(real*real+imaginary*imaginary)
Phase = Inverse_Tangent(imaginary/real)

and for calculating complex number from magnitude and phase we do this :

Real = Magnitude * cos(Phase)
Imaginary = Magnitude * sin(Phase)

So, what is the cause of problem. Is it a data scaling issue.

You need to change how you're calculating the phase from "atan" to "atan2".  The reason is that atan doesn't work in all quadrants.  More on atan2 here.


Why in formula of power of periodic signal is $2N+1$? Where does $1$ come from, like from $-N$ to $N$ is $2N$ so why is $2N+1$?

Isn't $-N$ to $N$ $2N+1$ index in total (you probably missed the 0 here)
$$ P = \lim_{N\to\infty} {\frac{1}{2N+1} \sum_{n=-N}^N} \lvert x[n] \rvert^2 $$


Correcting Impulse Invariance Method I'm trying to find out if the correction (Jackson, Nelatury, Mecklenbräuker) could improve the (IIM based) filter response near Nyqvist.

Here's my c++ routine (which is used for to calculate various RIAA and non-RIAA filters by just changing the timeconstant values and samplerate):

double a0, a1, a2, b0, b1, b2;
double fs = 44100;
//timeconstants (case RIAA):
// frequency -> time conversion 1/(2*pi*fc) (= R*C)
//poles
double p1 = 3180e-6; // 1/(2*pi*50.05Hz)
double p2 = 75e-6;    // 2212Hz
//zeros
double z1 = 318e-6;  // 500.5Hz
double z2 = 0.0;       // 3.18e-6 for Neumann pole (50kHz)

double pole1= exp(-1.0/(fs*p1));
double pole2 = exp(-1.0/(fs*p2));
double zero1 = exp(-1.0/(fs*z1));
double zero2 = exp(-1.0/(fs*z2));

a0 = 1.0;   // = 1.0
a1 = -pole1 - pole2; // = -0.931176
a2 = pole1 * pole2; // = 0
b0 = 1.0; // = 1.0
b1 = -zero1 - zero2; // = -1.731986
b2 = zero1 * zero2; // = 0.733838

Tried to google "ready to use" solution of this but the only source code I found few papers I could use for correction.

Papers:

Jackson
Mecklenbräuker
Nelatury

Juha

EDIT:

Bypassed the above mentioned papers. Improved the method so that now there's no need for additional correction biquad but by using z2 (which is unused) for the correction. One can decide where the error lies. Result:

You might be surprised to hear that your code implements neither the conventional impulse invariance method (IIM) nor the corrected IIM. Instead it implements the matched Z-transform, which maps the poles and zeros of the analog prototype according to

$$z_k=e^{s_kT}\tag{1}$$

where $T=1/f_s$ is the sampling interval.

The two versions of the IIM can be derived as follows. If we take as an example the RIAA de-emphasis transfer function

$$H(s)=\frac{sT_2+1}{(sT_1+1)(sT_3+1)}\tag{2}$$

with $T_1=3.18\,\text{ms}$, $T_2=75\,\mu\text{s}$, and $T_3=318\,\mu\text{s}$, we first have to do a partial fraction expansion of $(2)$ to get

$$H(s)=\frac{A_1}{s-p_1}+\frac{A_2}{s-p_2}\tag{3}$$

with $p_1=-1/T_1$, $p_2=-1/T_3$, and

$$A_1=\frac{1}{T_1}\frac{T_2-T_1}{T_3-T_1},\qquad A_2=\frac{1}{T_3}\frac{T_3-T_2}{T_3-T_1}$$

With the constants used in $(3)$, the transfer function of the transformed system can directly be written as

$$H(z)=\frac{TA_1}{1-e^{p_1T}z^{-1}}+\frac{TA_2}{1-e^{p_2T}z^{-1}}\tag{4}$$

From $(4)$ it can be seen that the poles are the same as for the matched Z-transform, but the zeros are different. Note that for this example there is only one zero for the matched Z-transform as well as for the conventional IIM.

The corrected IIM according to Mecklenbräuker and Jackson subtracts a constant term from the transfer function $(4)$ of the conventional IIM:

$$H(z)=\frac{TA_1}{1-e^{p_1T}z^{-1}}+\frac{TA_2}{1-e^{p_2T}z^{-1}}-\frac{T}{2}(A_1+A_2)\tag{5}$$

Note that this solution has one more zero than the other two methods. The poles of all three methods are identical.

The figure below shows the approximation of the RIAA de-emphasis filter given by $(2)$ according to the three aforementioned transformations (IIM in red, IIM corrected in green, matched Z in magenta). IIM performs worst, matched Z-transform and the corrected IIM are similar but their errors have different signs. This fact could be used to come up with a filter that combines the corrected IIM and the matched Z-transform. Since their poles are identical, we only need to combine their numerator coefficients. The result of the combined filter is shown in black. Its approximation error is less than 1 dB over the whole frequency range.


How to mix 8bit/8kHz unsigned linear PCM data in AppsScript? I am having the most difficult time mixing simple unsigned 8-bit/8kHz linear PCM audio samples in JavaScript (AppsScript).  I've tried all the basic maths listed here, employing basic number arrays.

Sign, add, clip, unsign.  Subtract 128 from each byte, add them, clip, then add 128.
Produce the Average, each sample just a[x]+b[x]/2
Viktor T. Toth's, maths linked.

Full example is below.  This will ask you to authorize Google Drive to access the audio samples (no it does NOT touch anything of yours):
https://script.google.com/macros/s/AKfycbzMOWV5Z5soG3uinrFc0WcscNXxsDq9elE9rzG0t703vsvDJZMS/exec

As you can see the 3 mixed results are downright horrible.  I've always assumed the unsigned bytes were just offset by 128.  Not the case?  Do I need to fold half the byte to get a true linear representation?  Does a basic number in JavaScript not translate appropriately into signed/unsigned representation?  If that's the case the Java/AppsScript is somehow properly doing it with the two real samples.

Update
Here is some specific code below, and the output of this is in the link above.  Alas, it simply doesn't quite work (byte folding in JavaScript)?

   var vegaArray = DriveApp.getFileById('0B-e9EqGm0pWPQ3RUTXFyUERDVTA').getBlob().getBytes();
   var fdraArray = DriveApp.getFileById('0B-e9EqGm0pWPaUJQUmFRQWctNG8').getBlob().getBytes();

   for(var i=44;i<vegaArray.length;i++)
     vegaArray[i] = Math.round( ( 2 * (vegaArray[i] + fdraArray[i]) )
       - ( vegaArray[i] * fdraArray[i] / 128 )
       - 256
     ); // for
   // vegaArray[] is the mixed output

Assuming that abyte[] and bbyte[] are two arrays with samples valued from 0-255, this is the solution:

combined = [];
for(i = 0; i < abyte.length; i++)
{
    sample =  2*(abyte[i] + bbyte[i]) - ((abyte[i]*bbyte[i])/128) - 256
    combined[i] = Math.round(sample);
}

I have a full example with both wav files extracted into arrays (I did this using matlab).  Also I verified that it plays back fine using matlab.

https://jsfiddle.net/4hsfo96g/


Coherence bandwidth and coherence time in an OFDM system While studying an OFDM system, I have encountered statements like coherence bandwidth is the reciprocal of delay spread and coherence time is the reciprocal of the Doppler spectrum.

Can somebody help me to understand this mathematically?

How is coherence bandwidth related to delay spread, and how is coherence time related to Doppler spread?

This diagram (from Molisch) shows the various Fourier relationships between measurements in wireless systems.


Granulometry by morphological opening and closing I am trying to understand granulometry. I understand that it deals with determining the size distribution of particles in an image and uses the principles of opening and closing.

I have a smoothed image with marbles in different sizes and call it A and a structuring element B, which is a disk with variable radii. I get that when opening the image with the structuring element B, smaller marbles disappear and by increasing the radius of B, bigger marbles start to disappear.

What I do not understand is what will happen if I use closing instead of opening. Do the bigger marbles disappear with a smaller radius B? So, does the exact opposite happen with closing?

When you use the opening, the overall operation is called a granulometry. However, when you use a closing, it's called an anti-granulometry. It's the dual operation. Combined together, they form the Pattern Spectrum.

In your case, if you apply a closing, you will certainly measure the structures between the marbles.

However, it is definitely not recommended to apply a smoothing first, because it will smooth the edges and small particles. Consequently, the particles detection will be less accurate.


Intuitive Meaning of Regularization in Imaging Inverse Problems Hello Every one I have been trying to understand the intuitive meaning of using a regularizer in images.  Specifically what does the Total Variation regularizer do in images and how is it able to prevent the edge smoothing?

One general form of Inverse Problem in Imaging which assumes Linear Operator is given by:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda R \left( x \right) $$

Where $ R \left( x \right) $ is the regularization function.

One way, which is pretty intuitive in my opinion, is to treat the above as a solution to Maximum a Posteriori Estimation (MAP) problem.

The first term, the fidelity term is the result of the assumption additive white Gaussian noise.
The regularization is the term of the Prior Knowledge on the image.

In the case of Total Variation it is usually given by:

$$ R \left( x \right) = \sum_{i, j} \left| {x}_{i + 1, j} - {x}_{i, j} \right| + \left| {x}_{i, j + 1} - {x}_{i, j} \right| $$

If you think on it as a Derivative (Vertical and Horizontal) then the MAP tells you have the assumption of Laplace Distribution for the Derivative.
It is well known that Laplace Distribution is something which promotes sparsity.
Namely the assumption is having few large values for the Derivative but mostly having values of zeros.
This fits perfectly with the Piece Wise Smooth model of Imaging.

Hence, in case of Total Variation Denoising the use of this regularization tries to smooth image while assuming that small variation of the derivative is caused by noise hence should be smoothed out and larger variations should be sparse.


How can I accurately represent an mp3's wave form in respects to time? I'm trying to programmatically plot a sound wave that I have placed in a 2d array, however, I cannot consistently get the sound wave to reach the end of the window. For example, when you load a sound in audacity, you get this:

In contrast, when I load a sound into my program I will often get something like this (red is the play head that moves at speed of one blue tick mark every second. The Green line is the end of the song in seconds):

As you can see, I can plot the complete wave form, but it does not sync up with the end of the song. Oddly enough, there are times during this song that the wave form is an accurate representation of what's being played. There are also a handful of songs that do not seem to have this issue (but more than a handful of others that do have this issue).

This is how I'm drawing the wave form:

    float drawStep = songLengthInSeconds / numberOfSamples;
    for(int i = 0; i < waveForm.length ; i++){
      for(int j =0 ; j < waveForm.length - 1; j++){
         renderer.line(x, waveForm[i][j], x + drawStep, waveForm[i][j+1]);

         x+= drawStep;
}

}

The renderer.line method signature is

renderer.line(beginningX, beginningY, endingX, endingY);

And I get the number of samples by just adding all of the waveForm[x].length together.

Finally, the top and bottom wave forms in the picture below are being drawn differently which is why they don't match. The bottom wave form is the one that is represented by the code above.

I've been fighting with this for a couple weeks now and would be very grateful if anyone has any suggestions that would push me in the right direction.

Storing the samples in a square array is weird. You can't represent all song lengths (in samples) as squares of integer array dimensions but yet the complete array is being plotted.

Your code seems to miss one x+=drawStep when transitioning between rows of the array. Neither does it "connect the dots" across the rows. If the last sample of each row is duplicated at the next row then it's fine as a row represents only waveForm.length - 1 samples.

Additionally, it just might be that x and drawStep are not accurate enough. Try making both doubles instead of floats. If there is an improvement but not yet good enough, calculate x by multiplication before each j loop. You should typecast songLengthInSeconds to double too.


Why does the flat part of the BER curve appear after the Viterbi decoder (MatLab)? I am trying to implement convolutional codes according to this MathWorks example in MatLab 2014a.

clear; close all; clc
rng default
M = 4;                 % Modulation order
k = log2(M);            % Bits per symbol
EbNoVec = (0:20)';       % Eb/No values (dB)
numSymPerFrame = 100000;   % Number of QAM symbols per frame

modul = comm.RectangularQAMModulator(M, 'BitInput', true);
berEstSoft = zeros(size(EbNoVec));

trellis = poly2trellis(7,[171 133]);
tbl = 32;
rate = 1/2;

decoders = comm.ViterbiDecoder(trellis,'TracebackDepth',tbl,...
'TerminationMethod','Continuous','InputFormat','Unquantized');

for n = 1:length(EbNoVec)
    % Convert Eb/No to SNR
    snrdB = EbNoVec(n) + 10*log10(k*rate);
    % Noise variance calculation for unity average signal power.
    noiseVar = 10.^(-snrdB/10);
    % Reset the error and bit counters
    [numErrsSoft, numErrsHard, numBits] = deal(0);

    while numErrsSoft < 100 && numBits < 1e7
        % Generate binary data and convert to symbols
        dataIn = randi([0 1], numSymPerFrame*k, 1);

        % Convolutionally encode the data
        dataEnc = convenc(dataIn, trellis);

        % QAM modulate
        txSig = step(modul, dataEnc);

        % Pass through AWGN channel
        rxSig = awgn(txSig, snrdB, 'measured');

        % Demodulate the noisy signal using hard decision (bit) and
        % soft decision (approximate LLR) approaches.
        demods = comm.RectangularQAMDemodulator(M, 'BitOutput', true, ...
        'DecisionMethod', 'Approximate log-likelihood ratio',...
        'VarianceSource', 'Property', 'Variance', noiseVar);
        rxDataSoft = step(demods, rxSig);

        % Viterbi decode the demodulated data
        dataSoft = step(decoders, rxDataSoft);

        % Calculate the number of bit errors in the frame. Adjust for the
        % decoding delay, which is equal to the traceback depth.
        numErrsInFrameSoft = biterr(dataIn(1:end-tbl), dataSoft(tbl+1:end));

        % Increment the error and bit counters
        numErrsSoft = numErrsSoft + numErrsInFrameSoft;
        numBits = numBits + numSymPerFrame*k;

    end

    % Estimate the BER for both methods
    berEstSoft(n) = numErrsSoft/numBits;
end

semilogy(EbNoVec, berEstSoft,'-o','LineWidth', 1.5)
hold on
legend('Soft','location','best')
grid
xlabel('Eb/No (dB)')
ylabel('Bit Error Rate')

And I am confused about the obtained curve:

Why does not the BER curve decrease exponentially after the 5dB?
Have anybody faced this issue before?

while numErrsSoft < 100 && numBits < 1e7

Nope, you can't stop after 10⁷ bits if you're after error rates < 10⁻⁶. That means you expect to see less than 10 bit errors. That's not a sufficient statistic!

So, this is where I know your curve becomes unreliable, anyway.

You might still be running into some error flow due to the implementation, but I can't look inside Matlab's implementation of neither the QAM decision nor its Viterbi, so these are "meh" at best for scientific considerations about the nature of Viterbi.


Convert fixed frequency bandwidth samples to time domain I'm trying to convert signal (bandwidth 1400 MHz, start freq 9300 MHz, stop freq 10700 MHz, step 3.5 MHz) from frequency domain to time domain. Signal data in frequency domain, as well as signal data in time domain, I got from network analyzer. So after processing data from freq to time I'm able to compare results. The problem is, that my calculation result differs from network analyzer data. You can see it on the picture below. Blue line - data from network analyzer, green line - my calculation result. Picture (a) is a signal in frequency domain, picture (b) - signal in time domain. In time domain step is 5e-10 seconds. Does anyone has any suggestions? Maybe Network analyzer freq data simply has not enough resolution to convert to time domain correctly as Network analyzer do? Thank you.

As I suspect, it was an issue with insufficient number of points in freq domain data. I convert time domain data from the green trace to freq domain and got exactly the same values in freq domain as on the first picture with original freq data from network analyzer.


Detrend data with no clear secular trend prior to Fourier analysis? I am completing Fourier analysis on many different time series of sediment particle flux exiting an experimental flume. Data is collected at a resolution of 1 Hz for durations ranging from ~5,000 to 50,000 seconds. Flux ranges from 0 to over 1,000 particles. The figure below shows an example time series that I am analyzing - the x-axis has units of seconds and the y-axis has units of particles/second:

I am familiar with time series and Fourier analyses, but I am not an expert. A colleague told me that I should detrend all of my time series before performing the Fourier analysis. I have completed the analysis without, and with detrending the time series and I see no obvious difference in the resultant spectra; this does not surprise me because I see no secular trend, for example, in the data shown in the example figure above. This leaves me with several questions. First, should I detrend the example time series shown in the figure above? If I should, why (or why not)? Second, there are at least several ways I can think of to detrend the example time series (e.g. subtracting a polynomial fit, subtracting the series mean, etc.), do people have any suggestions for a time series like the example one shown above? Last, are there clear criteria I can work with in making a decision of whether detrending is necessary prior to completing a Fourier analysis? Thanks.

How you pre-process data is related to the purpose of the analysis and since you haven't mentioned why you are doing spectral analysis, answering your question can't go beyond making some general remarks.

Detrending is a preprocessing operation. You want to compensate or remove some bias from your data. As an example, many kinds of data are collected in the presence of long term phenomena that is extraneous to what your are interested in finding. Many kinds of data have seasonal, diurnal, weekly, monthly, tidal, .... ect cycles associated with them. Unemployment has a seasonal dependence to it as an example.

I once was tasked with looking for jumps in syndrome data that is reported to health departments and there was strong increases on Mondays and after holidays in reports but that had to do with many doctor's offices being closed on weekends.  This was an artifact of how the data is collected and not how epidemics evolve.

Matlab's dtrend function removes a linear fit from data.  If the collection interval is small compared to some much longer trend in the data, it reduces the bias introduced by taking a section that is less than the period of the longer term trend, particularly for data like your which is strictly positive.  One may conceptualize dtrending as reducing low frequency spectral leakage that would otherwise mask features in your spectra.

Again, without knowing what you are looking for in your data, no specific answer can be recomended for your situation.


Fourier Transform Units It is documented that 'one' of the units of the Fourier Transform [of $x(t)$ volt] is volt per Hz. That is $X(\omega)$ components will have units of volt per Hz, where $\omega$ is the angular frequency (along the $x$-axis).

For example, the Fourier Transform of $\cos(\omega_0 t)$ volt has two components - one at $\omega = \omega_0$ and one at $\omega = -\omega_0$, each assigned a 'value' of $\pi$. So it is assumed that the values are '$\pi$ volt per Hz'.

Now, when the plot is adjusted to become a function of cyclic frequency $f$ (in Hz), which means plotting $X(f)$ versus $f$, then the Fourier Transform components of $\cos( 2\pi f_0 t)$ become '1/2 or 0.5' instead of '$\pi$'.

My question is --- why isn't the '$\pi$ volt per hz' preserved when translating between $X(\omega)$ and $X(f)$?

According to tables showing Fourier Transform pairs, the Fourier Transform of  $\cos( 2\pi f_0 t)$ volt has two components -- each with a value of '1/2 or 0.5'. For these values, does anyone know what are their units? I had assumed that the values for each component would have been 'preserved' (to be $\pi$ volt per Hz), regardless of whether we do the plot as a function of '$\omega_0$' rad/s or a function of '$f_0$' Hz. Thanks all.

For a given signal $g(t)$ with Fourier transform $G(\omega)$:

$$g(t) = \frac{1}{2\pi}\int_{-\infty}^{\infty} G(\omega)e^{i\omega t}\textrm{d}\omega$$

If you change variables such that $\omega=2\pi f$

$$g(t) = \frac{1}{2\pi}\int_{-\infty}^{\infty} G(2\pi f)e^{i2\pi ft}\textrm{d}(2\pi f)=\int_{-\infty}^{\infty} G(2\pi f)e^{i2\pi ft}\textrm{d}f$$

So if you present your Fourier transform as a function of $f$ instead of $\omega$, your Fourier pair of $g(t)$ will now be $\hat{G}(f)=\frac{G(2\pi f)}{2\pi}$ and:

$$g(t) =\int_{-\infty}^{\infty} \hat{G}(f)e^{i2\pi ft}\textrm{d}f$$

Note that with this definition, there is a scaling of $2\pi$ in comparison with the first equation I wrote. So now, for example, the Fourier transform of a cosine will consist of two deltas not with amplitude $\pi$, but with amplitude $\frac{\pi}{2\pi}=\frac12$, as correctly stated in your question.


FFT over a fixed and equal numbers I've got a vector of $100 000$ numbers. All numbers are equal ($7000$ for example).

If I perform FFT over this vector, what will I get?
From my understanding, I should receive a fixed DC line. Is this true?
Is there any website (i've searched) which can simulate that?

Because I did tried to run FFT with Python np.fft.rfft and np.fft.rfftfreq and didn't get it...

What am I missing here?

(My question is derived from a broader perspective, so bear with me...)

Yes it is true. Your FFT will compute those DFT values for k=0 to k=99.999 all of which are exactly zero but $X[0] = 7000 \cdot 100000 = 7 \times 10^8 $ for k=0. (the DC)

In order to see why this is the case, we can simply look at the definition of the DFT (as FFT is nothing but a fast algorithm to compute the result of DFT)

The DFT $X[k]$ of a signal $x[n]$ of length N is:

$$ X[k] = \sum_{n=0}^{N-1}{ x[n] e^{-j\frac{2\pi}{N}nk }} ,~~~~~~~~~~~~~~0 \leq n,k \leq N-1 $$

Now assuming the signal $x[n]$ is constant K:
$$x[n] = K ~~~~~~~~~~,~~~~~~ \text{for} ~~~~~~ 0\leq n \leq N-1 $$

The DFT becomes:
$$ X[k] = \sum_{n=0}^{N-1}{ K e^{-j\frac{2\pi}{N}nk }} ,~~~~~~~~~~~~~~0 \leq n,k \leq N-1 $$

There are a number of ways to evaluate or interpret this finite sum of a periodic complex exponential sum all of which will yield the following:
$$X[k]=
\begin{cases}
N \cdot K ~~~~, k = 0 \\
0 ~~~~, 1 \leq k \leq N-1 \\
\end{cases}$$


Why reprojection error is minimized instead of triangulation error for bundle adjustment? For bundle adjustment optmization why is reprojection error minimized instead of traingulation error or ray-to-3d point error? Disadvantage with reprojection error is: for perspective camera, its highly dependent on depth of the camera from the point in 3d right.

Triangulation error is not measurable in most cases - unless you have an alternate way to get at the 3D points being triangulated, you normally have no ground truth about them. Also, in normal situations where global scale is unknown, triangulation error is only defined up to scale.

Reprojection error is the gold standard error measure because it is directly measurable and it is unaffected by global scale,