How to model a generic low frequency signal? I'm trying to apply Fourier analysis to a specific problem I have.
I have essentially an integral like the following

$$
\int_{\Omega} f(t) g(t) dt
$$

And I'm trying to assume that $g$ is a narrow band signal (namely the Fourier transform is defined on a compact set). I want to prove that if $g$ is really narrow band then

$$
\int_{\Omega} f(t) g(t) dt \approx \left( \int_{\Omega} f(t) dt \right) \left( \int_{\Omega} g(t) dt \right)
$$

To prove this, I made this assumption

$$
g(t) = g_a(t) = \frac{1}{a} \int_{-\infty}^{+\infty} \hat{g}(\omega)rect \left( \frac{\omega}{a} \right) e^{-j2\pi \omega t} d \omega
$$

Where $a$ is a positive real parameter. When $a \to 0$ this integral turn out to be

$$
\frac{1}{a} \int_{-\infty}^{+\infty} \hat{g}(\omega)rect \left( \frac{\omega}{a} \right) e^{-j2\pi \omega t} d \omega \to \hat{g}(0)
$$

which means

$$
g(t) = g_a(t) \to \hat{g}(0) = \int_{-\infty}^{+\infty} g(t) dt
$$

and therefore

$$
\int_{\Omega} f(t) g(t) dt \approx \int_{\Omega} f(t) \hat{g}(0) dt = \hat{g}(0) \int_{\Omega} f(t) dt = \int_{-\infty}^{+\infty} g(t) dt \int_{\Omega} f(t) dt
$$

which is not exactly what I want , but it's close enough. The final question is:

is the model

$$
g_a(t) = \frac{1}{a} \int_{-\infty}^{+\infty} \hat{g}(\omega)rect \left( \frac{\omega}{a} \right) e^{-j2\pi \omega t} d \omega
$$

the correct model for a narrow band signal? I'm not entirely sure because of the factor $1/a$ I've introduced. But the idea why it makes sense is the following.

If the signal has narrow band, that means the variation in time is really slow, that means the narrower is the band the more the signal tend to be constant, and in theory if it would be constant then I'd be able to factorize the function $g(t)$ from the integral.

Clarification : probably in this discussion is more correct to say "low frequency".

Note that you can also use Parseval's theorem. Assuming that $f(t)$ and $g(t)$ are both real-valued we have:

$$\int_{-\infty}^{\infty}f(t)g(t)dt=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)G^*(\omega)d\omega\tag{1}$$

If $g(t)$ is a narrow-band low-pass signal, i.e., if we can model $G(\omega)$ as $G(\omega)\approx G(0)\text{rect}(\omega/\Delta\omega)$, and if we can assume that $F(\omega)$ is approximately constant (equal to $F(0)$) in the narrow interval $\Delta\omega$, then we can write

$$\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)G^*(\omega)d\omega\approx\frac{1}{2\pi}F(0)G^*(0)\Delta\omega\tag{2}$$

where $\Delta\omega$ is twice the bandwidth of $g(t)$ (because the bandwidth is defined as the support at positive frequencies). With

$$F(0)=\int_{-\infty}^{\infty}f(t)dt$$

and

$$G^*(0)=G(0)=\int_{-\infty}^{\infty}g(t)dt$$

we obtain from $(1)$ and $(2)$

$$\int_{-\infty}^{\infty}f(t)g(t)dt\approx\frac{\Delta\omega}{2\pi}\int_{-\infty}^{\infty}f(t)dt\int_{-\infty}^{\infty}g(t)dt\tag{3}$$


Circular wrapping of an asymmetric function (in DFT calculations) Convoluting a signal (using discrete FT) for a given interval [a, b] with a  Gaussian can be done by circular wrapping as shown in Numerical Recipes. I found a shortcut for circular wrapping so that the phase remains zero for an even-N data. Consequently, the peaks in the signal will not shift their positions after convolution. The shortcut is to center the Gaussian at the mean (a,b)+dt, where dt is a sampling interval. After that, one can do a circular shift by N/2.
Consider a hypothetical case where we have a skewed Gaussian instead of a Gaussian. Should I have the centroid at the mean (a,b)+dt before circular shift by N/2?
Continuous version

Circularly wrapped version

The key is being even, notice that $f_{N+k} = f_{k}$, and when you do a circular wrapping of an infinite function you should have $f_{k} = \sum_{i=-\infty}^\infty f((i \cdot N + k) T_s)$, then naturally $f_{N-k} = f_{k}$, in practice you want to choose $f$ and $N$ so that $f$ vanishes before $N/2$, then you are safe to use only one term of the summation that define $f_{k}$.
This holds true for asymmetric signals, the important point is that if the convolution result has more than $N$ non-zero points then then you can only recover the samples that are less than $N$ indices from both ends of the signal.
In other words if you have a signal that is non-zero from $L \le i \le H$ then you can only recover the samples at $i$ such that $H - N < i < L + N$. If you have a response of length not greater than $N$ then you can recover all the outputs. If your response has length $N + d$ you can only recover $N - d$ samples from the output. And by circularly shifting by the inputs by different number of samples you can control which samples of the output will be usable. you can also split your operators in multiple chunks to process an infinite duration signal as in overlap-save method.
An example in python
There is a pair of useful function fftshift and ifftshift to perform these operations without unintended code obfuscation.
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import skewnorm;
a = 3

plt.figure(figsize=(6, 6))

# circularly shifted equally spaced grid
t = np.fft.ifftshift(np.linspace(-10, 10, 2048));

# filter
h = skewnorm.pdf(t, a)
# signal without circular shift
x = np.zeros((2048))
x[500:1300] = np.arange(-400, 400);
# apply filter to signal (result is circularly shifted, you can
# recover the signal in the natural order with ifftshift)
y = np.fft.irfft(np.fft.rfft(h) * np.fft.rfft(x)) / len(t)


Uniform sampling of a periodic signal According to the Nyquist sampling theorem, a sampling frequency for a periodic signal $x(t)$ with a period $T$ and bandlimited by $f_m$ Hz is given by

$$ f_s > 2f_m. $$

Is it possible to relax this constraint for periodic signals?

You can, in a sense, if you can also sample your signal at an incommensurate sampling rate for the same period.

Any periodic signal can be decomposed into its constituent tones with a DFT and a proper algorithm.  (I have an algorithm I developed about 10 years ago that I have in the pipeline for my blog that I recently found out that McLeod described about 20 years ago.)

So, reduce your question down to a single pure tone.  The frequency for a pure tone, or an alias tone, can be found exactly (see my blog if you doubt this) for any frequency.

So, reduce your question down to a pure tone with a whole number of cycles in your sample frame.

The limitation with a DFT due to the discrete sampling is the aliasing effect.  If a tone is higher than the Nyquist frequency it will appear in the DFT as a frequency less that the Nyquist limit.  So, if a tone appears in bin $k$, where $k < N/2$, the will be a conjugate value in bin $N-k$, which is over the Nyquist limit.  In actuality, you can't tell with the frequency is $k$ cycles per frame, or $N-k$.  The frequency could also be any multiple of N over those.

The possible frequency can be: $k$, $N-k$, $k+N$, $2N-k$ ...

By sampling at two different truly incommensurate rates, the sets of possible frequencies from each of the DFTs will only have one value in common, and that will be the actual frequency.  If the rates are nearly incommensurate, the next possible common value will be well out of range.

Here is a concrete example:

60 Hz Sampling  N = 60

49 Hz Sampling  N = 49

             49                   60

 5 Hz   5, 44, 54, 93 ...    5, 55, 65, 115 ...

47 Hz   2, 47, 51, 96 ...   13, 47, 63, 107 ...

97 Hz   1, 48, 50, 97 ...   23, 37, 63, 97  ...

Now, suppose you thought you could be clever and use every other point of your 60Hz sampling, and every third point as your two different sampling rates.

One second sample frame

30 Hz Sampling  N = 30

20 Hz Sampling  N = 20

Possible Frequencies:

             20                                30

 5 Hz   5, 15, 25, 35, 45, 55, 65 ...   5, 25, 35, 55, 65 ....

47 Hz   7, 13, 27, 43, 47, 63, 67 ...  13, 17, 43, 47, 73 ....

You can't tell whether the 5 Hz is 5, 25, 35, 55, etc.

You can't tell whether the 47 Hz is 13, 43, 47, etc.

So that doesn't work so well.

Hope this helps.

Ced

Followup:

This technique assumes you have "ideal" instantaneous sampling, so it is more theoretical than practical.


Understanding where the constant $2/N$ comes from in Fourier transformation I'm implementing Fourier transformation in my analysis and I wanted dig a bit deeper on the reasons why the absolute value of Fourier transformation is usually multiplied by the constant $2/N$ to get the peak amplitude value of a sinewave with certain frequency.
In the book Understanding Digital Signal Processing by Lyons, the author states the following relationship between the peak amplitude $A$ of a sinewave and the output magnitude $M_r$ of the discrete Fourier transformation (DFT) for that particular sinewave is:
$$M_r=AN/2,\;\;\;\;\;\; \tag{1}$$
where the $r$ stands for real input values to DFT and $N$ is the number of input values to DFT. From this relationship, I trivially get the amplitude I want to as $A=2M_r/N$, which I see is many times done in many fft-examples in Matlab found throughout the web.
Now my big question was, why is the relationship in $(1)$ true? I started to read more from the book Fourier Analysis and Its Applications by Folland and I found the following in his book (in section about DFT):
$$\widehat{f}\left(\frac{2\pi m}{\Omega}\right)\approx \frac{\Omega}{N}\widehat{a}_m,\;\;\;m=0,1,...,N-1\;\;\;\;\;\;\tag{2}$$
where $\widehat{f}$ is the amplitude function, $\widehat{a}_m$ is the $m^{th}$ output of DFT, $N$ is again the number of inputs and $\Omega$ is the length of the time interval $[0, \Omega]:$
$$\widehat{f}\left(\frac{2\pi m}{\Omega}\right)=\int_0^\Omega e^{-2\pi i m t/\Omega}\;f(t)\;dt,$$
where $f$ is the wave function. Now when I look at $(1)$ and $(2)$, there seems to be a connection between them:
$${\color{red}{\widehat{a}_m}} \approx{\color{blue}{\frac{N}{\Omega}}}{\color{green}{\widehat{f}\left(\frac{2\pi m}{\Omega}\right)}},\;\;\;\;\;\;{\color{red}{M_r}}={\color{blue}{\frac{N}{2}}}{\color{green}{A}}.$$
These two results are almost satisfying but I wondered why it seems to be the case that $\Omega=2$?
My questions: Where does this $2$ come from? Why in Lyons's book there is $2$ instead of $\Omega$?
I thought could it be somehow related to the symmetry of the DFT output? One time unit to left and right: $[-1,1]$ so the length of the interval would be $\Omega=2$? A bit vague this last part but could I be onto something here?
UPDATE:
The definition for DFT in book Understanding Digital Signal Processing is given as:
$$X(m)=\sum_{n=0}^{N-1} x(n) e^{-2\pi i nm/N},$$
where $x(n)$ is some continuous time-domain signal. In the book Fourier Analysis and Its Applications the corresponding definition is:
$$\widehat{a}_m = \sum_{n=0}^{N-1}a_n e^{-2\pi i mn/N}\;\;\;(0\leq m<N),$$
where $a_n = f\left(\frac{n\Omega}{N}\right)$.

The two comes from the fact that a real valued sinusoid is really the average of two complex one.

$$ \cos( \theta ) = \frac{ e^{i \theta} + e^{-i \theta} }{2}  $$

There is the "2" in the denominator.  In the DFT, the other half is located at bin $N-k$, if $k$ is the bin number.

That definition of cosine comes straight from Euler's equation:

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

The $N$ comes straight from the definition of the DFT.  The relationship beteen the amplitude, N, and the magnitude of the DFT bin is why I like to use the $1/N$ normalized form of the the DFT over the more conventional unnormalized form.

I recommend that you read my first blog article The Exponential Nature of the Complex Unit Circle for an understanding of Euler's equation and my blog article DFT Bin Value Formulas for Pure Real Tones for understanding the bin values of a pure tone and what "leakage" really means.  Equation (19) is the answer to your "big question".

The $\Omega$ comes from the continuous FT.  The DFT is the Discrete Fourier Transform.  Too often, they are conflated leaving confusion like yours.  It takes a lot of heavy math to understand the relationship between the two.  The premise of my blog articles is that the DFT can be learned and understood straight from the summation definition without any reference to the continuous case.

Hope this helps.

Ced


Filtering RANSAC Estimated Homographies I am using RANSAC algorithm for homography estimation between pairs of images taken with cameras which do not have any translation between them (pure rotation and change of scale/zoom). It works well in half of the cases. The correct output looks like this:

The red lines are filtered correspondences and the quadrilaterals illustrate how the homography distorts perspective.

Sometimes, however, many bad cases happen, like these:

I already have a simple test in the RANSAC loop. It makes a simple quadrilateral (a unit square) and transforms it with sample transform. Then it looks whether the transform kept its convexity.

Still, however, bunches of concave quadrilaterals come out.

Do you have any idea how to properly test the homography, if it behaves "nicely" and filter out the incorrect solutions?

I found some code where they test that none of the three transformed points are colinear. But this does not seem sufficient as it will not filter out deltoids and other "invalid" quadrilaterals...

There is a problem in checking whether the homography is OK.

The algorithm for checking correct homographies may interest someone, so I will write it down here:

1) Create a quadrilateral $ABDC$ with vertex coordinates (in homogenous coordinates):

$$\begin{eqnarray} A:& (-w/2,-h/2, 1.0) \\ B:& (w/2,-h/2, 1.0) \\ C:& (-w/2,h/2, 1.0) \\ D: &(w/2,h/2, 1.0) \end{eqnarray}$$

where $w,h$ are width and height of the image, respectively. If the whole image frame (a rectangle) is transformed to convex quadrilateral, than any convex quadrilateral within it will also be transformed "neatly".

2) Create transformed quadrilateral $A'B'D'C'$ in which every vertex is transformed using the computed homgraphy (e.g. $C'=HC$ ). From now on, all points will be converted to non-homogenous coordinates.

3) Compute vectors $\vec{u}$, $\vec{v}$ for parametric representation of diagonals:

$$\begin{eqnarray}d_{1} :& A+(D-A)s =A+ \vec{u}s \\ d_{2} :& B + (C-B)t=B+\vec{v}t \end{eqnarray}$$

The intersection of diagonal comes from $d_{1}=d_{2}$:

$$t=\frac{1}{d}\left[(B_{y}-A_{y})\vec{u}_{x} - (B_{x}-A_{x})\vec{u}_{y}\right]$$

$$s=\frac{1}{d}\left[(A_{x}-B_{x})\vec{v}_{y} - (A_{y}-B_{y})\vec{v}_{x}\right]$$

Then the convex quadrilateral satisfies $s,t \in (0,1)$.

In practice, one can introduce a fudge factor to avoid not only non-convex and degenerate quadrilaterals, but also near-degenerate ones, like when three points are near colinear. So the test can be modified such that $s,t \in (\lambda, 1.0-\lambda)$, where lambda is a small number (e.g. $\lambda=0.01$).

Older problem, fixed in the above algorith:

I found the problem here - having a certain homography, the test can pass for a smaller quadrilateral, but not for the larger one. This is why some "ill" homographies passed through.

The green square represents a source image, the orange is a transformed one. As you can see, the left hand one is convex, but starts deforming as the source is larger:

Finally even larger source yield non conver quadrilateral:

I made a mistake with scaling. The points in homogenous coordinates $(x,y,w)$ were scaled in $x$ and $y$ direction, but in $w$. This is why the same transform produced different quadrilaterals.

I have corrected the algorithm accordingly.


Group delay of the integrator I have a set of sensors, some velocity-meters and others accelerometers. To process the two groups together, I integrate the accelerometers signals (in time domain, using MATLAB's cumsum), so I'll have all sensor data now in velocity. Knowing the precise arrival time of a signal at different sensors are very important in  determining the location of the signal source. Wondering if I am introducing some phase shift by integration, I tried to estimate the group delay. I am not sure what I am doing is right though. Given that the integrator can be written as:

$$y(n) = x(n) + y(n-1)$$

the filter coefficients will be

\begin{align}
a&=\begin{bmatrix}1 &-1\end{bmatrix}\\
b&=\begin{bmatrix}1 & 0\end{bmatrix}
\end{align}

Giving these values to the MATLAB's grpdelay function (Group delay), I get a constant value across all frequencies of half a sample.

Does it mean that my signal after integration is shifted by half a sample regardless of frequency ?
And what about zero frequency (the mean value), why is the group delay 0.5 sample at 0?

Your integrator is a perfect accumulator, which is an unstable (marginally stable) system with a pole at $z=1$:

$$H(z)=\frac{1}{1-z^{-1}}\tag{1}$$

This means that the value of the frequency response at $\omega=0$ (corresponding to the pole at $z=1$) is undefined. Consequently, the value of the group delay at $\omega=0$ is also undefined. Note that you can write down the frequency response of an ideal accumulator using a Dirac delta impulse (see this question and its answer), but this doesn't help in finding a value of the frequency response (or the group delay) at $\omega=0$.

For $\omega\neq 0$ (or, more generally, $\omega\neq 2\pi k$), we can easily determine the frequency response and group delay of the system:

$$H(e^{j\omega})=\frac{1}{1-e^{-j\omega}}=e^{j\omega/2}\frac{1}{e^{j\omega/2}-e^{-j\omega/2}}=e^{j\omega/2}\frac{1}{2j\sin(\omega/2)},\quad\omega\neq 2\pi k\tag{2}$$

From $(2)$, we see that the phase $\phi(\omega)$ equals $\omega/2$ plus a constant term which is irrelevant for the group delay. Consequently, the group delay is given by

$$\tau_g(\omega)=-\frac{d\phi(\omega)}{d\omega}=-\frac12,\quad\omega\neq 2\pi k\tag{3}$$

which is what you found using Matlab. However, Matlab doesn't tell you what's going on at $\omega=0$, where the value of the group delay is indeed undefined.


SNR estimation prior to demodulation? What are the standard ways to measure SNR, particularly before demodulation/synchronization?

Perform a 2 dimensional cross correlation against a noise free replica of the signal you expect to receive with delay offset on one axis and frequency offset on the other, and find your maximum correlation. Normalize that to get $\rho$ and from $\rho$ you can get SNR. Do this after channel equalization, otherwise you will get the SNR degradation due to ISI as well.

For further details, see:

Noise detection

and

How can I find SNR, PEAQ, and ODG values by comparing two audios?

If the SNR is high and the noise is well behaved (thermal noise floor, no band selective interference present) then you can make a reasonable estimate directly from a spectrum plot such as the example I show below, where the SNR in this case was 40 dB as is clear from the distance between the signal and the noise. This is signal power to the total noise power that is in the signal bandwidth.


Echo removal from a chunk of PCM signal I have implemented a paper regarding echo detection in real-time using MFCC coefficients. Now, I want to remove the echo from the mixed signal. Let me explain it as:

There are two signals, NearEnd signal (Signal Containing Echo) $N$ and Farend Signal (Original signal which causes echo when played) $F$
Echo is being detected in NearEnd Signal $N$. This signal also contains my voice.
Now if we want to remove the $F$ from $N$, is it simple like the addition of two PCM signals?

For example, If I want to mix two PCM signals, research papers show it can be done in multiple ways but the easiest way is :
$$M=A+B$$
Now, if I want to remove F from N, and we get a new echo-free signal $S_{new}$, then will it be like
$$S_{new} = N -F$$

Your $S_{new}=N-F$ requires that you know both $N$ and $F$.
But since $N=F*R$, with $*$ being convolution, and $R$ being the Room Impulse Response (RIR), you're solving the equivalent problem $S_{new} = F*R-F = F*(R-\delta)$, with $\delta$ being the delta dirac impulse or its discrete equivalent.
The hard part here is estimating $R$, which aligns nicely with your initial statement that you want to do echo detection in itself, not signal reconstruction. So, you're building a RIR estimator for itself, not for usage as part in an equalizer. That makes a lot of sense: Maybe a heuristic based on observing the MFCCs as they are, then "cleaning them up" in MFCC domain (simply nullifying things that shouldn't be there), then estimating a time-domain signal that would have these corrected MFCC signals, then using that as estimate $\hat F(t)$ for $F(t)$ and the received time-domain signal $N$ in a classical IR estimation algorithm, giving you a RIR estimate and thus an echo description.
Since convolution with the RIR is an LTI system, trying to do the estimation itself in MFCC domain feels like a bad idea – what is easy in time domain and linear frequency domain suddenly becomes hard: MFCC's are great if your system is not linear, but make working with linear systems harder!
If you insist on integrating MFCCs into your detection process, you could try to use them  as measure of speech quality, and feed that knowledge back into an adaptive time domain filter.
Generally, if I only had MFCCs, I'd try to reconstruct a time-domain signal from that, and then apply the mature theory of impulse response estimation on that – there's plenty of results that achieve optimum performance under some given conditions, and if you have one of these, whatever else you could do can only be worse.


Downsampling PCM Audio I'm looking everywhere for solutions for audio downsampling. However for one reason or another they are incomplete, sometimes buggy, and sometimes are in a programming language I cannot understand.

So here's my use case:

I have an array containing two channels (left-right, it is a microphone) with Float32 data that gets interleaved and then converted to PCM like this:

function interleave(left, right) {
    var length = left.length + right.length;
    var result = new Float32Array(length);
    var _index = 0;
    for(var index = 0; index < length;) {
        result[index++] = left[_index];
        result[index++] = right[_index];
        _index++;
    }
    return result;
}

function convertToPCM(raw) {
    var output = new Buffer(raw.length*2);
    for (var i = 0, offset = 0; i < raw.length; i++, offset += 2){
        var s = Math.max(-1, Math.min(1, raw[i]));
        output.writeInt16LE(s < 0 ? s*0x8000 : s*0x7FFF, offset);
    }
    return output;
}

This is NodeJS\Javascript, but I think it's pretty clear how it works.

Now here's the issue. The output of these function needs to generate a WAV file downsampled to 22000 or 16000, but the source has a sample rate of 44100.

Could you please highlight what algorithm should I use in order to achieve this?

Thank you very much!

Breaking down the conversion into an upsampling and a downsampling step as explained in mbaz's answer is not necessary - band-limited interpolation achieves the same result in one single process, for any conversion ratio, and its quality is easy to fine-tune (through various approximations of the band-limited pulse).

libsamplerate is an example of open-source implementation that does not rely on upsampling/downsampling. It should be possible to bind it to NodeJS through an add-on, or even to rewrite it in javascript - once stripped of the code handling corner cases (variable rate, compensation of processing delay), the code fits in a 50-line loop.


Downconverting real signal to IQ baseband Despite reading various answers on this site and on other sites I am still confused about how to convert a real-valued signal back down to its original IQ baseband values. Note I am still a beginner at this.

I have a very minimal prototype in Matlab but I'm getting stuck. I begin by generating some IQ data and modulating it onto a carrier signal, ready to be sent as a stream of real valued numbers.

% constants
f = 100; % carrier frequency
Fs=500; % sampling frequency
t=(1/Fs:1000)*T; % time domain

% generate 1000 raw I/Q values to *transmit*
I=randi([0 1],1000,1)';
Q=randi([0 1],1000,1)';

% modulate I/Q values onto a carrier wave
Iwave = I.*cos(2*pi*f*t);
Qwave = Q.*sin(2*pi*f*t);

I assume that adding the two waves together is correct (and not subtracting them like I've seen in some notations)

% add up I/Q. and transmit this stream of real-valued numbers
tx=Iwave+Qwave;

Now the receiver mixes the signals with cos and sin waves

rxI = tx .* (cos(2*pi*f*t));
rxQ = tx .* (sin(2*pi*f*t));

Below is an FFT of the received real-valued signal in black peaking at 100Hz. And the mixed I and Q signals on the receiver, peaking at 200Hz (I know the red peak is hard to see, but it's there at 200Hz).

Now comes the tricky part. Every resource I've read mentions I should low pass filter this I and Q, and I'll be able to get back my baseband data.

cutoff=200/Fs;
order=32;

% a 32 order low-pass FIR filter with cutoff at 200Hz
filter=fir1(order,cutoff);

filtI=filter(filter,1,rxI);
filtQ=filter(filter,1,rxQ);

Here's the FFT of the I and Q signals after low-pass filtering

Now the sent and received I channel as shown below obviously look nothing alike and simply cannot be demodulated.

What am I doing wrong here?
The ideal answer would specify exactly what filter I need to get back my baseband data in Matlab. This I think is reasonable since the example I provided is ideal and simple.

Please also correct any misunderstandings in this post aside from the obviously incorrect filter I am using.

The allowed bandwidth of the data created and used to modulate the IQ signal needs to be below the cutoff frequency of the receiver post-heterodyne low pass filter, which needs to be below the carrier frequency, which needs to be below half the HF sampling rate to avoid aliasing.

So try a lower data rate and/or a much higher frequency carrier, with a suitably corresponding higher frequency sample rate.  You can lower your data rate by using a slower changing rand statement, and filtering that to remove harmonics of the squarish pulses.


Low-pass filter with varying cut-off frequency over time Is there a novel way to create a low-pass filter (Finite Impulse Response Filter) with varying cut off frequency over time?

For example, I have a white noise that is 10 seconds long sampled at 44.1 kHz. I want to create a low-pass filter that will have a varying cut-off frequency over time. At time 0 my cut-off frequency would be 100 Hz and will decrease linearly down to 1 Hz at time 10 seconds.

I wouldn't call this novel but... I'm going to give an example on how to do this for a specific filter with a specific discretization.

A 1st order butterworth low-pass filter is given by:

$H(s) = \frac{Y(s)}{X(s)} =\frac{\omega_c}{s + \omega_c}$

Where $\omega_c$ is the cut-off frequency of the filter.

The backward difference (aka backward euler) approximation is given by:

$s \approx \frac{1-z^{-1}}{T_s}$

Where $T_s$ is the sampling time. By substitution, the discrete transfer function becomes:

$H(z)  = \frac{Y(z)}{X(z)} \frac{\omega_c}{\frac{1-z^{-1}}{T_s} + \omega_c} = \frac{\frac{\omega_c \cdot T_s}{1+\omega_c\cdot T_s}}{1 - \frac{1}{1+\omega_c\cdot T_s}\cdot z^{-1}}$

Let $a = \frac{\omega_c \cdot T_s}{1+\omega_c\cdot T_s}$, then:

$H(z) = \frac{a}{1 - (1-a)\cdot z^{-1}}$

Thus the difference equation is:
$y[n] = a\cdot x[n] + (1-a)\cdot y[n-1]$

Why go through the tedious process of getting an equation for $a$ when I could have just used the above formula and try to guess some values for $a$? Well, now if you want to change the cut-off frequency $\omega_c$ in real time, you simply have to recalculate $a$ each time you compute the signal. This, of course, assumes you have a constant sampling period.

Note that this filter is an IIR filter, not FIR as you requested, but the logic is basically the same: as long as you know how your coefficients relate to the cut-off frequency, then you just simply have to recalculate them every time you want a different cut-off frequency. Nothing stops you from repeating the same steps for different discretization methods, different types of filters, etc...

You may even wish to work purely in the discrete-time domain, but the relevant part is always knowing how your coefficients relate to the cut-off frequency.


How does complex data of an FFT and IQ modulation differ? If I have an arbitrary time domain signal (e.g. a simple cosine) i can sample it with an ADC and get a time discrete representation.
I can perform an FFT on the samples and get a complex output.
The phase of the signal is:

$$\phi=\arctan\frac\Im\Re$$

Using an IQ down modulation of the same signal gives complex data, too, with the phase being:

$$\phi=\arctan\frac QI$$

Is there a difference between these two and when is one preferred over the other?

EDIT

Thanks for the input!

@Dilip Sarwate: Here is what I did: set $\phi=60°$ (see your example)

The results are exactly what I expected. I can see the frequency of 1 Hz in the spectrum (plot 3) and the initial phase of $60°$ in the phase spectrum (plot 4). I get the same result when I use a non-complex function, e.g. $x(t)=cos(2\pi t+\phi)$.

I read everywhere that phase information is lost if I don't use IQ down-modulation, but that does not seem to be the case, as I can see the phase information in the plots.

Or am I missing something here?

If the signal frequency is constant and exactly integer periodic in the FFT width, then the transform basis vector for each FFT result bin is the same as a sinusoidal (complex exponential) IQ down-modulator mixer input that starts with a phase of zero at sample 0.  Thus the atan2() results for a signal within that FFT result bin should be identical.

If you only need one result bin of a DFT, then a complex mixer (complex Goertzel) requires less computation than a full FFT result (assuming equivalent numerical precision, etc.).  If you need O(logN) DFT result bins or more, then the FFT is more efficient.

If the frequency is not known to be exactly integer periodic in the DFT/FFT width, then an FFT might be useful as part of an initial frequency estimator algorithm.


Incorrect Butterworth filter outcome using Rational Transfer Funcion I've been playing with Second-Order Butterworth high-pass filter for a while.
I've calculated the coefficients using Matlab function and I wanted to use them to filter simple sine wave generated by Matlab.
Unfortunately the results are not correct/consistent, verified by Matlab filter(b,a, signal) function.
My Matlab code is as following (ommiting plotting part):

clear all;
fp = 10000; %sampling frequency
ts = 0:1/fp:0.05; %number of samples, length
squareF=100; %frequency of signal
squareWave = sin(2*pi*squareF*ts); %generate the sine
fc = 500;
[b,a] = butter(2 ,fc/(fp/2), 'high'); %numenator, denumenator
for n = 3:501
    filteredWave(n) = (b(1)+b(2)*squareWave(n-1)+b(3)*squareWave(n-2))/(a(1)+a(2)*squareWave(n-1)+a(3)*squareWave(n-2))*squareWave(n);
    %Rational Transfer Funcion
end
filteredCorrectly = filter(b, a, squareWave);

Filter plots are below.

Where is the error?

The output can be calculated using the coefficients as follows:

$$a_0y(n)+a_1y(n-1)+a_2y(n-2)=b_0x(n)+b_1x(n-1)+b_2x(n-2)$$

In your formula you are not doing this. Instead, you are doing a division with terms that only involve the input and not the output, which doesn't make sense as the filter is IIR.

I suspect that you got confused with the transfer function definition of an IIR filter, which is given by

$$H(z)=\frac{Y(z)}{X(z)}=\frac{\sum\limits_{k=0}^2 b_kz^{-k}}{\sum\limits_{k=0}^2 a_kz^{-k}}$$

Note that you are using an equation in the $z$-domain applied to time-domain signals (to put it simply). That's why the loop you wrote is not working.


What are the known use cases for different wavelet families? There are numerous wavelet families that differ by a number of parameters (like the number of vanishing moments, symmetry, etc.).

What are the known use cases where one particular family is known to be more suitable and provide better results?
For example, Daubechies 9/7 are used in the JPEG2000 image format. As far as I understand Daubechies wavelets were designed specially to have a high number of vanishing moments so that makes them more suitable for signal compression tasks. What are the known use cases for other families of wavelets?

And as second part:

How does one decide which one is the most suitable for it's task?
Are there any known heuristics, common knowledge, etc?

Answer 0: ask yourself if you really need wavelets

Say yes. Let us concentrate on 2-band real discrete wavelet first. JPEG 2000 is a special case, where CDF9/7 and 5/3 biorthogonal wavelets are used. For the first one, the analysis wavelet is moderately regular, but its moments nicely concentrate piecewise regular parts of the data. The synthesis wavelet is smoother, and visually compensates quantization artifacts. The 5/3 works similarly, with the addition of exact computations on integers. There are almost symmetric, which is somehow better for edge location. The paper Mathematical Properties of the JPEG2000 Wavelet Filters gathers some of the most notable.

Answer 0.5. apart from that example, there is very few "fit-for-all signals" optimality results.

The quality of wavelets depends a lot on the nature of the data: signals, images, on the task, on the resource constraints. But the most important is the mastery of the wavelets you use. Because most often, you are trying to enhance or restore data, but the target is unknown, and quality metrics are often imperfect,  so you can only get informed guesses on simulated datasets. Theoretical properties are guides, but you need models for signals or noise at least (eg spline wavelets with a fractional Brownian noise).

Orthogonal wavelets, says Daubechies', are widely used because orthogonality is useful to prove stuff, and their "shortest" support interval for given moments is nice. On natural signals, a rule of thumb is that it is rare to observe polynomial parts above degree 3. Hence short DB should suffice. However, longer ones are often used to compensates for the stronger asymmetry of the shortest.

The (overall) symmetry of the wavelet is also observed: anti-symmetric to act like derivatives, symmetric to mimic Laplacians.

There is a no-go theorem of discrete 1D 2-band wavelets: except for Haar, they cannot be real, orthogonal with symmetry and finite support. Haar seems crude, but it is very fast, and since you can almost compute everything (delays, etc.) it can do a marvelous job.

But if you go deeper, you can find that $M$-band wavelets $M=4$ for instance, can have all the required properties, with more degrees of freedom. And they have better aliasing cancelation between channels, which makes them better at textures, for denoising.

Answer 1: a good knowledge of  wavelets properties restricts the candidate list, and a little trial-and-error is often unavoidable

But you can also use a set of different wavelets to enhance the result of a single one: perform the processing with each wavelet, and (wisely) combine the result. One of the most efficient technique is based on a pair of Hilbert-related wavelets: take any wavelet basis you like (for its properties), and build a second basis with an (approximate) Hilbert transform relationship. The resulting primal and dual wavelets share the same moment property,  have an opposite symmetry, and form together a dual-tree wavelet transform, that provides almost shift invariance at little cost. But you can also use different shifted versions (time-invariant wavelets), symmetrized combinations, etc.

Answer 2: a union of wavelet bases is often a powerful combination.


DC value of $x(t)=\frac{1}{t}$ Is it possible to calculate the DC value of signals with undefinable area?

More specificaly, in the case of $x(t)=\frac{1}{t}$. $\int\limits_{-\infty}^\infty\frac1t \, \mathrm dt$ does not converge. Does that mean that its DC can not be determined?

You need to be careful with the definition of the DC value of a signal. The actual time average, which is often called DC value is given by

$$\overline{s(t)}=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{-T/2}^{T/2}s(t)dt\tag{1}$$

whereas the value of the signal's Fourier transform at $\omega=0$ is given by

$$S(0)=\int_{-\infty}^{\infty}s(t)dt\tag{2}$$

If $(2)$ is finite, $(1)$ equals zero, and if $(1)$ is finite but non-zero, $(2)$ doesn't exist, and the Fourier transform $S(\omega)$ has a Dirac impulse at $\omega=0$. See also this related answer.

For the given signal $s(t)=1/t$ both integrals in $(1)$ and $(2)$ do not converge in the conventional sense. However, their Cauchy principal value exists and equals zero, as mentioned in a comment by Robert Bristow-Johnson. Hence, the DC-value of $s(t)=1/t$ equals zero, regardless whether you define it by $(1)$ or by $(2)$.

Note that the Fourier transform of $s(t)=1/t$ equals

$$S(\omega)=-j\pi\,\textrm{sgn}(\omega)\tag{3}$$

which is just the frequency response of a Hilbert transformer (scaled by $\pi$).

From $(3)$, $S(0)$ doesn't exist, but the average of the left-sided and right-sided limits $\frac12 (S(0^-)+S(0^+))$ equals zero, and - using the Cauchy principal value according to the definition of the Hilbert transform - the Hilbert transform of a constant is indeed zero.


Fourier transform possible on non-rectangular part of an image Dear Signal Processing readers,

I want to introduce 'noise' into parts of images.

Until now, I worked with rectangluar images for a similar purpose and did the following, using (inverse) Fourier transform (2D) with MATLAB:

fft = fft2(image);    % run fourier transform on image
amplitude = abs(fft); % extract power spectrum>
newPhase = ( (rand(size(amplitude)) * 2) - 1) * pi;          %generate random phase
newImage = real(ifft2(amplitude. *  exp(1i * newPhase)));    %use both to generate noised image

This randomizes the phase all over the image, preserving the power spectrum and should guarantee the same spatial frequency spectra in input and output.

What I want to do now, is applying this effect to only a part of an image, which is not rectangular (its the "inner" parts of a face and will fit in an ellipse).
So I want to "cut out" this ellipse from an image, turn it to "phase noise" (with power spectrum/spatial frequencies constant) as above, and but it back into the "hole" I cut out of the source image before.

My "prime solution" of course would be that FFT/iFFT could operate on ellipses as well, so I could keep using my method shown above. After some amount of google-ing I am quite pessimistic about that option - is there any hope?

I found the following traces to possibly promising approaches, but as I am not very familiar with the math behind FFT, I'd appreciate any hint/explanation (why) which of these ideas go into a good direction:

Somehow "stretch" the ellipse to a rectangle shape, apply my method as is, strech the resulting rectangle back to the input ellipse shape. Here I fear that the rectangle/ellipse-transformation might destroy the matched power spectrum. Do you know good way of transforming without loosing this information? (This led me to this idea)

Two more hints came from this post:

Using a 3D-fourier transform, in which (as far as I got it) the 3rd dimension gives information about which part of the image shall be used. Any idea / experience how to implement this?
Using "Elliptic Discrete Fourier Transform" - but I do not understand what (and how) it does and how it can be implemented.
In this forum, someone suggested to a related question some sort of window tapering (which, as far as I got it, is made to "filter" the influence that only comes into the image because the ellipse of interest is on a black background) that I did not understand as well, plus i am not sure if it could help me because it did not sound like that question would require an ellipse as output.

Thanks in advance for any suggestion!

Cheers. Jan

Even when the image part is rectangular, the FFT result will be contaminated by boundary effects. The FFT of the part is the Fourier transform of an infinite tiling of the rectangular part. (An elliptical or circular part cannot be tiled. A hexagonal tiling would get pretty close though.) If there is a discontinuity going from one tile to another, this will create high-frequency components that are not descriptive of the insides of the part. Randomizing the phases will spatially translate these components so that they appear also inside the part, which is probably not what you want. Mirroring at the boundaries does not eliminate discontinuities of the derivative, and requires doubling of the FFT dimensions.

This is where a smoothly tapering 2-d window function can help. For a circular part you might like the Tukey window with the distance from the part center as its argument.

For an ellipse, stretch the window function accordingly. Multiply the data by the window function to select the area you want to analyze, and FFT the result. Then do the phase randomization. Mask the noise by multiplying it with an ellipitical disc-like mask that has value 0 outside the disc and value 1 inside the disc. You can also use the window function as the mask to get a fuzzy boundary. I did not find yet a perfect way to normalize the noise, but matching the mean and the standard deviation both weighted by the mask works quite OK:

Multiply the input image by the complementary mask: 1 - mask, and sum the two together. Octave source:

pkg load signal;
r = 120; # Circle radius
n = 512; # Image width and height
tukeyr = 0.5; # Window transition relative size
lena = double(imread("lena.gif"))/255;
tukeyn = round(r/(1-tukeyr));
x = repmat([1:n], n, 1);
y = repmat([1:n]', 1, n);
dist = sqrt((x-n/2).^2 + (y-n/2).^2);
tukey = tukeywin(tukeyn*2, tukeyr)(tukeyn+1:tukeyn*2);
win = reshape(interp1(1:tukeyn, tukey, reshape(dist, n*n, 1), 0), n, n);
randphase = 2*pi*rand(n, n);
noisefft = abs(fft2(lena.*win)).*exp(j*randphase);
noise = real(ifft2(noisefft));

#Circular mask
weight = r-dist;
weight = (weight >= 1)*1 + (weight > 0 & weight < 1).*weight;
#Or use window as mask
#weight = win;

weightsum = sum(sum(weight));
noiseweightedmean = sum(sum(noise.*weight))/weightsum;
noiseweightedstdev = sqrt(sum(sum((noise-noiseweightedmean).^2.*weight))/weightsum);
lenaweightedmean = sum(sum(lena.*weight))/weightsum;
lenaweightedstdev = sqrt(sum(sum((lena-lenaweightedmean).^2.*weight))/weightsum);
normalizednoise = ((noise - noiseweightedmean) / noiseweightedstdev * lenaweightedstdev) + lenaweightedmean;
imshow(lena.*(1-weight).+normalizednoise.*weight);

P.S. I found that highpass filtering the copy of the Lena image used for noise generation, by setting some of its lowest-frequency bins to zero, will reduce the error in the magnitudes of the fft of the final image compared to fft of the source image. I think the reason is that low-frequency amplitude fluctuations will be reduced. Within the circle (ellipse) those appear in practice as zero frequency, and in their absence less normalization is needed. So that's one thing to investigate.


How to choose the most appropriate color space for a computer vision application? A lot of image processing techniques used in computer vision consist (among other things) in switching from RGB to another color space (HSV, YUV, LMS...) : color transfer, visual tracking... It seems that different people use different color spaces for the same application.

In the case of statistical color transfer for example, you can use LMS (which is a device-independent color space) or use YUV (a device-dependant color space) and achieve pretty much the same result. Moreover the transformation from RGB to either of these 2 color spaces is achievable by multiplying the RGB pixel by a 3x3 matrix.

Besides the device-dependence/independence factor, there's the decorrelation between the channels. The RGB channels are inter-correlated but so are LMS.

My questions are:

What are the characteristics to look for in a color space in order to be able to say whether it is appropriate or not to use it?
(If you think that the scope of my question is too broad, take the example of statistical color transfer from an image to another)
And why is the decorrelation between the channels of any color space that important?

Computer vision is a huge topic and your particular example (statistical color transfer) is too narrow to provide guidance on the remaning parts.

That being said, classically most algorithms have worked on black and white "color" space, which was either obtained directly from the sensor or was later converted to it. The reasons behind is: 1-It's much simpler, 2- It provides remarkably good results for most applications unless color is explicitly required.

Color is, however, inevitably required for certain machine vision applications. Then comes the problem of selecting the most suitable color space. RGB is the fundamental color space because those R,G,B signals are physically captured by the sensor. (This also came from human color perception and additive color formation based on it)

Correlation in RGB channels should be understood clearly. In general given an R value you cannot say anything about G or B. But given a brightness value and hue, then RGB will (may) be correlated.

YUV is the most useful color space used in TV broadcast technology, image and video coding applications (codecs, cameras, camcorders) and many technical applications where color is necessary. For some reasons YUV was derived based on special requirements of color tv broadcasting over existing black and white tvs with compatibility in mind. Y channel carries the brightness of the scene and typically used to replace the black and white image. U and V channels carry color value and saturation information.

HSV is a more complicated color space to derive from RGB. It can provide more accurate color discrimination than YUV. The HSV is definining image aspects like Hue, Saturation and Brightness which are fundamentally uncorrelated (or less correlated compared to RGB)

For technical reasons other color spaces also exist but they do not provide any advantages for signal processing and especially for machine vision.


Audio signal processing I am trying to code audio signal processing software in java.

Purpose: It will use audio signal processing library to filter out noise from the raw PCM based input (.wav) of heartbeat recording.

Steps involved

Select input .wav file.
Apply Low Pass filter at 33 Hz.
Apply High Pass Filter at 2 Hz.
Apply Amplification with Max Gain.
Store back it as .wav file

There are many C++ library which will cover most of the requirement. But looking for Java for this particular project. Need input from Audio Signal Processing engineers for their input if these steps will be sufficient to able to hear the heartbeats. Also please let me know which order of filtering and which window will be most effective in this scenario.

Input is taken from

http://www.lehman.edu/faculty/hoffmann/itc/techteach/audio/sound2.html

Your low pass should be higher - probably 40 or so.
Your high pass can be higher - around thirty
Your gain should be variable else you will amplify it into distortion - what you want is called normalisation.  You will need to find the maximum volume of your signal and calculate a scaling factor that will bring this up to full scale.  (Say full scale is 1.0 and your maximum value is 0.01 then multiply all samples by 100.)

When I experiment with heart sounds, I find the beat itself to be around 34Hz - I usually do a band pass from 30 to 40 in order to isolate just the beat.

If you do this, you will have a clean recording of ONLY the heartbeat.  If all you want is to count heartbeats or listen to them you will be fine.  If you need to listen for heart defects or other things, then you will need to change the parameters.  You will need higher frequencies and a wider bandwidth for those types of things.

One other thing you need to consider:  30 to 40 Hz doesn't really come through headphones and ears all that well.  Even when you isolate the hearbeat and amplify it to full scale, it is still hard to hear.  Earphones will attenuate it some, and your ears even more so.

As to actually doing the filtering:

You can use IIR filters at higher orders to isolate the heart sounds.  Avoid using chebyshev filters, though.  The "ring" badly and can make one beat "smear" into the next.
If you use an FIR, you may want to downsample first.  It depends on your sampling rate.  If your sampling rate is high (say the typical 44100) then an FIR than can bandpass something around 30Hz will have to be VERY long.  If you downsample to 1000 Hz, then shorter filters can be used.  You could get away with maybe 100 taps at 1000Hz.  Downsampling consists of doing a low pass filter to remove everything above half of your new sampling rate, then "throwing away" a lot of your samples - downsampling by two would throw away every second sample, downsampling by 4 keeps one out of four samples.

When I do audio in Java, I usually import libpd, and use pd to develop the actually processing stuff.  Then the same pd "patch" can be used with libpd in the java programm.
pd  libpd

The order in which you apply the filters (for your  bandpass) shouldn't matter.


How do I get the volume levels from an audio file? I need a Linux tool to analyze Audio files for their Loudness. Maybe the average for each 100 milliseconds.
output like: (0=>-12,100=>-18,200=>-24,300=>-55)
Any Ideas

I can give you a quick and hacky solution with sox that can be easily installed on any Linux distribution.
sox in.wav -n trim 0 0.1 stats : newfile : restart 2>&1 | grep 'RMS lev dB' | awk '{ print $4 }'
This will return RMS levels for each 0.1 second chunk. Change in.wav to the input file and 0.1 if you require a different frame length.


How to Calculate the number of DFTs taken in a STFT Given variables of window size (W), overlap (O), and length of the Signal (S) (in my case, length of a song) and I am trying to determine an equation that gives me the number of DFTs that would be taken in an STFT as the window slides given those variables above.  I'm not sure if those variables is enough to give me the equation or not, and if it is, I'm not sure how to derive that equation.

I will first produce a lower number based on the covering of all the non-zero sample set.

With a window of size $W$, counting samples from $0$ to $S-1$, the first window containing at least  one sample covers $[-W+1,0]$. If the discrete overlap range is a number of samples $O<W$, the next window will be $$[1-O,W-O]=[-W+1,0]+(W-O)$$ and the others $k$th hops are $$[-W+1,0]+k(W-O)$$
Those intervals exceed the last sample $S-1$ when $0+k(W-O)\ge S-1$ or

$$k\ge \left\lceil\frac{S-1}{W-O}\right\rceil$$

so the covering takes $ \left\lceil\frac{S-1}{W-O}\right\rceil+1$ windows, hence DFTs.

Does this suffice in general, as we did not use any window property? For instance, if the window is non-null at one point only (or vanishes), the $O$ cannot be chosen arbitrarily.

The question in general seems open to me. For randomly chosen complex filterbanks, the smallest oversampling seems to suffice for recovery. For other tasks, who knows?


Concept of energy threshold I am analyzing some program responsible for detecting the correct signal in the channel on the receiver side. It seems to me that such an implementation allows for the detection of the correct packets (because I mean packet transmission, from transmitter to receiver) from the received samples. It is known that under normal conditions the energy of the signal is greater than the energy of the noise, so the simplest algorithm to find the beginning of an incoming packet is to measure the energy of the received signal. When no packet is received, the received signal consists of only noise. When a packet begins to be received, the energy of the received signal is increased by the signal component, so the packet can be detected as a change in the energy level of the received signal.
So, the program I am analyzing looks like this (it's MathScript in LabView, it's almost identical to Matlab):

So, the example only divide the samples into groups with size “Bin Size”, and then take the sum of each
group, so that we can get the average power. The logic function can be interpreted as this:

%Keep the powerin the lower fractions of bins
low_lim = max(1,floor(0.1*bin_num));
%z = sort(power);
%noise_power = mean(z(1:low_lim));
noise_power = min(power);
signal_power = max(power);
%signal_power = mean(z(bin_size-low_lim:end));

%energy_threshold = max(.01, min(0.009, noise_power*1000));
energy_threshold = max([noise_power*100, (noise_power+signal_power)/2, 1e-6]);

noise_estimation_error = 0;
if 20*log(signal_power/noise_power) < 3
   noise_estimation_error = 1;
end

Now the one issue I don't fully understand is "Energy Threshold". It seems to me that the value 1E-6 was chosen based on brute force analysis (but correct me if this is some value that is often used to determine noise or something like that, because I have no idea), as some sort of reference point. Noise_power is just that. well, the noise power (although I don't know why multiplied by 100, it is possible that a higher order of magnitude allows a more fair detection - I don't know, sadly)
But the thing that bothers me the most is the value of (signal_power+noise_power)/2 - as far as I can see, it's a simple arithmetic average of signal power and noise power. Is this how energy threshold is calculated? What is it then responsible for? Is it just information about where the packet begins to be received and where the noise ends? The processed signal is later "chopped" of unnecessary noise in the receiver, so that only the proper signal containing the packet comes to the demodulation blocks (well, usually it is so, but we know how it works in practice). I've heard that such an implementation allows for greater efficiency, because it relieves the algorithms used in synchronization - they don't have to process portions of the signal, which e.g. are not the actual packets but just part of the noise. This speeds up processing of packets and simplifies their reception.
So is this the thing? Does the "Energy Threshold" carry information about the end of the noise and the beginning of the packet due to the calculated arithmetic mean? Or am I not fully understanding something here?

Does the "Energy Threshold" carry information about the end of the noise and the beginning of the packet due to the calculated arithmetic mean?

There's no such thing as "end of noise": Nothing says that noise can't take values that are as large (or larger) than your signal's power; in "benign" noise, these values are just rare. (And that especially applies for large $L$: Exactly when your noise estimator becomes bearable, i.e. when $L$ is large and thus $N/L-1$ is small, you get potential for larger sums.)

Or am I not fully understanding something here?

Yep, you need to actually consider noise a random variable, and consider its power a random variable, too.
You need to define a noise model. For example, the often-used assumption of Gaussian noise leads to a distribution with fat tails – which means that, yes, noise will eventually be larger than any signal power, if you just look often enough.
The threshold was chosen for some reason that's impossible for us to know – maybe putting an upper bound to the probability of detecting noise as packet, or an upper bound to the probability of dismissing signal as noise (notice that in general, you don't know your signal's power, either, and it's a random variable, too).
Maybe someone was just lazy and put any value there. My crystal ball is cloudy – you need to ask whoever designed this. (Or your advisor. Asking your advisor/lab supervisor/professor at this point might be wise. You're missing the whole "Probability in Communications" chapter from the book I linked to in my answer to your question of yesterday. Maybe re-read that chapter? It's not actually sufficient in mathematical depth to derive the result I foreshadow below, but it's a good start.)

Noise_power is just that. well, the noise power (although I don't know why multiplied by 100, it is possible that a higher order of magnitude allows a more fair detection - I don't know, sadly)

It's very much not the noise power!
So, "noise power" is a property of the stochastic process (that means: random variable with a time index); you're only estimating that. This is is hence a noise power estimate, not the noise power.
To make matters worse: this is a biased estimate (in any reasonable noise model I can think of). It's probably not very close to the actual noise power.
Let's quickly derive that. Let $Z_t$ be your noise at (discrete) time $t$, a random variable. I don't know your noise model, and I don't think you've given us that info, so we'll have to just say that $Z_t$ has cumulative density function (cdf)
$$F_{Z_t}(z) := P(Z_t \le z),$$
which is nice enough to also imply a probability density $f_{Z_t}(z):=\frac{\mathrm d\,}{\mathrm d\,z}F_{Z_t}(z)$. We'll also have to assume that your noise has a zero mean, and finite variance $\sigma_{Z_t}^2$ (finite power!) (which applies especially to e.g. Gaussian noise). Because anything else would blow up this answer (and not because that is a good representation of every real-world noise), we'll also have to assume our noise is $L$-blockwise stationary and independent (which applies especially to white noise).
We will now look at the case where you get only noise in all segments.
Let us call the powers in your $N/L$ segments $S_i, i=0, \ldots, N/L-1$ (note: you've got a typo in your figure, if you start at 0, the number of the last segment is $N/L-1$, not $N/L$):
$$S_i =  \sum_{k=i\cdot L}^{i\cdot (L+1)-1}\frac1L \lvert Z_k \rvert^2$$
Out of curiosity, let's consider the cdf of $\lvert Z_k \rvert^2$:
\begin{align}
F_{\lvert Z_k \rvert^2}(v) &= P\left(\lvert Z_k \rvert^2 \le v\right)\\
&=\begin{cases}
0 & v < 0\\
P\left(\lvert Z_k \rvert \le \sqrt v \right) & v\ge 0
\end{cases}\\
&=\begin{cases}
0 & v < 0\\
P\left(-\sqrt z \le Z_k  \le \sqrt v  \right) & v\ge 0
\end{cases}\\
&=\begin{cases}
0 & v < 0\\
P\left( Z_k  \le \sqrt v  \right) -P\left( Z_k  \le -\sqrt v  \right) & v\ge 0
\end{cases}\\
&=\begin{cases}
0 & v < 0\\
F_{Z_k}\left( \sqrt v  \right) - F_{Z_k}\left( -\sqrt v  \right)  & v\ge 0
\end{cases}
\end{align}
with a mean $E\left(\lvert Z_k \rvert^2\right) > 0$ (a absolute square is always larger than zero, unless the original variable is 0, and I'd guess this is not a case you care about, because it means "noiseless"). With a little more derivation it's possible to show that the expectation of the random variable $\lvert Z_k\rvert^2$ is the variance of the random variable $Z_k$, i.e. the actual noise power $\sigma_Z^2$; the variance of $\lvert Z_k\rvert^2$ is twice that, at least for Gaussian $Z_k$.
Now, the cdf of that $S_i$ is simply:
\begin{align}
F_{S_i}(p) &:= P(S_i \le p)\\
&=P\left(\sum_{k=i\cdot L}^{i\cdot (L+1)-1}\frac 1L \lvert Z_k \rvert^2 \le p \right)\\
S_i &\overset{\text{CLT}}\sim \mathcal N\left(E\left( \lvert Z_{iL} \rvert^2 \right);
\frac 1L E\left(\left[\lvert Z_{iL} \rvert^2 - E\left(\lvert Z_{iL} \rvert^2\right)\right]^2\right) \right)\\
&= \mathcal N\left(\sigma_Z^2; \frac{2\sigma_Z^2}L\right)
\end{align}
assuming that we can apply the central limit theorem (CLT), which says that for sufficiently large $L$, the sum over a random variable (such as $\lvert Z_k \rvert^2$) has a distribution that approaches the Normal distribution $\mathcal N(\mu_S;\sigma_S^2)$. The factor of $1/L$ in $\sigma_S^2$ stem from the linearity of the expectation (and the resultant quadratic influence of constant factors on variance).
Now, this just shows that your $S_i$ are just a bunch of approximately normal variables.
So, let us consider the minimum of these $S_i$,
$$M := \min_{i=1,\ldots, N/L-1} S_i,$$
which itself is a random variable.
Now, the cdf of that is pretty simple:
\begin{align}
F_M(m) &:= P(M\le m)\\
&= P\left(\min_{i=1,\ldots, N/L-1} S_i \le m\right)\\
&= P\left(S_1 \le m \wedge S_2 \le m \wedge \ldots \wedge S_{L/N-1} \le m\right)
\end{align}
Because we defined the noise to be independent, and thus especially between blocks, $P_1$, $P_2$ and so on are independent variables. Thus,
\begin{align}
F_M(m) &= P\left(S_0 \le m \wedge S_1 \le m \wedge S_2 \le m \wedge \ldots \wedge S_{L/N-1} \le m\right)\\
&= P\left(S_0 \le m \right)P\left(S_1 \le m \right) P\left( S_2 \le m \right)  \ldots P\left( S_{L/N-1} \le m\right)\\
&= P\left(S_0 \le m \right)^{N/L}\\
&= \left(F_{S_0}(m)\right)^{N/L}.
\end{align}
We want to know what the expecation of $M$ is, because that's what we expect our noise power estimator estimates as noise power.
For that, we first find the pdf of $M$, and then just insert that into the definition of expectation.
\begin{align}
f_M(m) &= \frac{\mathrm d\,}{\mathrm d\,m} F_M(m)\\
&= \frac NL \left( F_{S_0}(m)\right)^{N/L-1} \frac{\mathrm d\,}{\mathrm d\,m} F_{S_0}(m)\\
&= \frac NL \left( F_{S_0}(m)\right)^{N/L-1} f_{S_0}(m)\\
\end{align}
Normalizing $m$ allows us to use the cdf of the standardnormal distribution:
\begin{align}
f_M(m)
&= \frac NL \left( F_{S_0}(m)\right)^{N/L-1} f_{S_0}(m)\\
&= \frac NL \left( \Phi\left(\sqrt{\frac{L}{2}}\frac{m-\sigma_Z^2}{\sigma_Z}\right)\right)^{N/L-1} f_{S_0}(m)\\[2em]
E(M) &= \int_{-\infty}^\infty m\cdot f_M(m) \,\mathrm dm \\
&= \int_{-\infty}^\infty m \left(
\frac NL \left( F_{S_0}(m)\right)^{N/L-1} f_{S_0}(m)
\right)\,\mathrm dm
&|\text{ partial integr.}
\end{align}
In the partial integration step we see that we indeed get a component that is the actual $\sigma_Z^2$, i.e. the actual noise power, but we subtract a strictly positive term. This bias means we're *systematically underestimating the noise power with what you call noise_power. This isn't a good estimator for noise power.
(There's a solution to the above integral; if you really want to know more, the $N/L$th power of a $\chi^2$-distributed variable involves a ratio of $\Gamma$ functions. Don't really think that result helps you much, but it proves (aside from the general consideration yielding the same) that you're underestimating noise power.)

But the thing that bothers me the most is the value of (signal_power+noise_power)/2 - as far as I can see, it's a simple arithmetic average of signal power and noise power. Is this how energy threshold is calculated?

No. Not in general, and I doubt it's  a good choice here.
That threshold $\gamma$ needs to be chosen based on how much missed packets you can live with vs how much noise misidentified as packet you can accept. A common choice for a threshold would be one that leads to maximum likelihood estimation of whether there's a packet $D$ or only noise considering received segment power $S_i$.


How to classify filters? In class, we have defined four types of filters: low-pass, high-pass, band-pass and band-stop. I have understood that in order to classify a filter into one of the fourth, one can use the signals:

DC: $1,1,1,\ldots$
Nyquist: $-1,1,-1,1,\ldots$
half-Nyquist: $1, 0, -1, 0, 1,\ldots$

Can you help me understand how to classify $y_n = x_n + x_{n-2}$ as a band stop filter?

Can you help me understand how to classify

Yes!

Just draw the (discrete) spectra of your three test signals; then multiply them  piece-wise with the frequency response of your candidates, and compare with the response of your filter.


ROC of inverse system function If the region of convergence (ROC) for system function $H(z)$ is $R_h$, what is the ROC of the inverse function $G(z)=\frac{1}{H(z)}$?

$$$$

I'm going to go out on a limb and answer this based on my gut instinct.

Let's suppose $H(z)$ is casual and has a single zero (a) and a single pole (b):

$$
H(z) = \frac{(1-az^{-1})}{(1-bz^{-1})}
$$

if the system is casual, then we know that the sequence is right-sided with ROC that is exterior of the circle. Thus, the ROC of H(z) is:

$$b < |z|$$

Now for the part I was uncertain about: the ROC of the inverse system G(z).

$$
G(z) = \frac{1}{H(z)} = \frac{(1-bz^{-1})}{(1-az^{-1})}
$$

When we take the inverse of H(z), it is necessary for a portion of the ROC of G(z) to overlap with a portion of the ROC of H(z) in order for an inverse function G(z) to exist. The overlapping doesn't necessary need to include the entire ROC just a portion of it.  This is the constraint that determines if G(z) has ROC that is casual / exterior-of-circle,   or anti-casual / interior-of-circle.  Its a matter of noting the poles of G(z) and checking the following cases to determine direction of the inequality for ROC:

(1) does the exterior of the largest pole of G(z) overlap with a portion of ROC of H(z)? if so then this represents a valid causal ROC (that is an exterior circle). ROC: a < |z|

(2) does the interior of the smallest pole of G(z) overlap with a portion of the ROC of H(z)? if so then this represents a valid anti-casual ROC (that is an interior circle). ROC: |z| < a

it is entirely possible that both a causal inverse and an anti-casual inverse both exist and are valid.

Another consideration for choosing a valid inverse function of H(z) is to make sure that the inverse function is stable, ie. the ROC includes the unit circle, which enables the system to be realizable in real hardware, and is more than a mere mathematical abstraction.


Need intuition on the frequency response of a signal I'm currently taking a course on signals and systems, and even though we are much ahead of frequency response, I still feel like I am missing an intuitive grasp of what the frequency response of a signal is. Specifically, when I look at the various CTFT pairs, it's difficult for me to intuitively understand why the pairs make sense.

If someone could provide an intuitive explanation of what frequency response is, that would be much appreciated!

Thanks

The frequency response tells you how a system reacts to a sinusoidal excitation of a given frequency, in terms of amplitude and phase.

As any signal, periodic or not, can be described as a sum of sinusoids, this allows you to deduce the response to any signal, by the superposition principle.

Interesting features of the frequency response are the continuous gain, the presence of resonances and the bandwidth. A flat curve reproduces the original excitation without deformations.

You get a feeling of a frequency response when looking at how a slowing spin-drier vibrates.


Why can't we just make all wireless networks use integer multiples of base frequency? I always wondered why transmission capacity depends on bandwidth. For example, let us say that there is an isolated island. In this island, people decide that all wireless networks use frequencies that are positive integer multiples of one base frequency. Then when we do Fourier transform, we can just use least common multiples of all used frequencies and data in each frequency will exactly be determined except background noises. Or am I wrong here? In this sense, bandwidth seems just to be useless....

It's not that easy to distinguish different frequencies when they are close to each other so you'll need to make sure that there is sufficient spacing between the frequencies. Also, you can't just keep transmitting higher and higher frequencies. At some point there will be a limit of the highest frequency that can go through a wire or whatever transmission mechanism you are using and also the fastest computer that you have to process it. So once you cap the highest allowed frequency your method comes down to dividing up the remaining space among a finite set of frequencies. In other words, you are proposing to allot a uniform bandwidth to each of your channels, which is exactly what you said you were trying to avoid.


Analysis of vibration signals I have a vibration signal coming from a motor measured from accelerometer (Irregular Time series) and I need to use these signals for analysis purpose.

I am in need to transform this signal to some form for analysis purpose. I read few blogs and websites containing vibration signal processing technique. The basic thing is to transform the signal from time domain to frequency domain.

So I used the fft function to transform using R

X.k <- fft(signal)

plot.frequency.spectrum() plot a frequency spectrum of a given fft [Link: here]

plot.frequency.spectrum <- function(X.k, xlimits=c(0,length(X.k))) {
  plot.data  <- cbind(0:(length(X.k)-1), Mod(X.k))

  plot.data[2:length(X.k),2] <- 2*plot.data[2:length(X.k),2]

  plot(plot.data, t="h", lwd=2, main="",
       xlab="Frequency (Hz)", ylab="Strength",
       xlim=xlimits, ylim=c(0,max(Mod(plot.data[,2]))))
}

The frequency spectrum of this fft looks like this

The first peak is extremely high. I don't know why this happens as am new to signal processing techniques. I cannot use this for analysis purpose and do prediction as it will end in bad results. Any other ways to do this method? Or do I need to use other techniques?

The first peak is the DC component of your spectrum which is large compared to the AC components (your signal doesn't have negative value and its always above zero which lead to large DC bias)

Find average of your signal over time and subtract it from your signal to remove the DC component.

Good luck.


What is the type of noise in this image? what reason could make the noise/effect in this image?

Below is the original

Thanks!

Looking at it, my thought is that the top image has been printed out (onto poor quality paper), then later scanned.  The artefacts would then be a mixture of paper texture and poor printing.  It could also include ink transfer from the facing page (if part of a book or newspaper, and/or bleed-through from something printed on the back of the paper.


Why is there only one integration in the solution if there is two integral in the formula? In this problem the random variable is theta and according to the formula there should be two integrations but in the solution there is only one . Nor am i able to understand the meaning of x1 and x2 in the formula. this problem is driving me nuts..

You have a single random variable $\Theta$ in this example, so taking the expectation with respect to that random variable results in a single integral. The formula for $R_X(\tau)$ in your question is rather messy and full of mistakes; it should be as follows:

$$\begin{align}R_X(\tau)&=E[X(t+\tau)X(t)]\\&=E[a\sin(\omega_0t+\omega_0\tau+\Theta)\cdot a\sin(\omega_0t+\Theta)]\\&=\frac{a^2}{2}E[\cos(\omega_0\tau)-\cos(2\omega_0t+\omega_0\tau+2\Theta)]\\&=\frac{a^2}{2}\frac{1}{2\pi}\int_{0}^{2\pi}[\cos(\omega_0\tau)-\cos(2\omega_0t+\omega_0\tau+2\Theta)]d\Theta\\&=\frac{a^2}{2}\cos(\omega_0\tau)\end{align}$$


Reason of circular conv in OFDM In OFDM system, after IFFT a CP is added at the beginning of the OFDM symbol by copying it from the end of the symbol.My question is that with CP why we use cyclic convolution instead of linear convolution. I would appreciate if someone please could explain the reason of cyclic convolution here.

Many thanks!

OFDM systems are usually implemented using the FFT and its inverse. At the transmitter, information symbol values are assigned to the various OFDM subcarriers, resulting in a single OFDM symbol in the frequency domain. To obtain a time domain that can actually be transmitted over the channel, an inverse FFT is applied. The receiver on the other end of the link performs the reverse (it computes a forward FFT) on its observed signal, resulting (hopefully) in the original frequency-domain symbols, thus accomplishing OFDM demodulation.

Why add a prefix?

The short answer: multipath through the communication channel. You can think of multipath as a linear filter applied to the transmitted signal via propagation through the communication channel. Thus, the signal observed at the receiver is the original one from the transmitter, convolved with the channel response (and corrupted by noise). This results in intersymbol interference (ISI), which degrades system performance.

Techniques to undo the effects of this channel response are called equalizers, and are widely used in practical communications systems. However, their implementation and analysis can be somewhat complicated. It would be nice if, for your OFDM system, you could come up with a simpler equalization method.

Another way to mitigate multipath would be to insert some sort of guard time between symbols (a prefix) in order to allow the delay spread of the channel from one symbol to die out before the next. This strategy is used in many OFDM systems; the added prefix allows time between symbols for the channel's impulse response to die out, reducing or eliminating ISI. For most other systems, adding this overhead on each symbol would be prohibitively inefficient. However, since OFDM's effective symbol rate is much lower (since many modulated subcarriers are transmitted in each OFDM symbol time in parallel), the overhead of this approach is proportionally lower as well.

Why use cyclic prefixes?

Recall the fact that the OFDM symbols are modulated onto subcarriers in the frequency domain. Also recall the relationship between convolution and the discrete Fourier transform (DFT): convolution in the time domain corresponds to multiplication in the frequency domain. Therefore, the effects of the convolution with the channel response show up as pointwise multiplication of the frequency-domain subcarriers. This makes for a very simple equalizer structure: we just have to multiply each subcarrier by a (per-carrier) constant value at the receiver in order to undo the effects of multipath. This is typically called a "one-tap" equalizer, and is equivalent to assuming that any fading is flat across each subcarrier's bandwidth.

However, recall that multiplication in the DFT domain corresponds to cyclic convolution in the time domain. So, in order for this strategy to work, we would need some way of turning the channel response into a cyclic convolution, not a linear convolution as typically happens in the real world. That's where the cyclic prefix comes in.

The implementation: for each time-domain symbol generated by the transmitter, a cyclic prefix is inserted before it. That cyclic prefix takes the last $L$ samples of the $N$-sample time-domain symbol and concatenates it to the front of the symbol to be transmitted. Therefore, for each $N$-sample OFDM symbol, the transmitter actually emits $N+L$ samples. The length of the prefix $L$ is chosen to be greater than the expected multipath spread.

Why does it work?

With the cyclic prefix inserted, you can see conceptually how the linear convolution induced by the channel turns effectively into circular convolution; the cyclic prefix, which includes the end of the symbol, is smeared onto the beginning of the symbol by the channel's impulse response. This effect is quite literally defined as circular convolution, so the goal of achieving simple equalization is accomplished.


$\tt awgn$ and $\tt wgn$ functions in MATLAB always output the same signal, no matter what I input as SNR I'm very new to all of this.

Right now I'm just trying to add a bunch of different noises to an input signal and see how they affect my signal.

This is the current output of my code:

 signal,Fs = audioread('~/Documents/Methlab/sp03-0.wav');
 %generate white noise from input signal
 [m,n] = size(signal);
 snr = input('SNR: ');
 power = input('Power[DbW] ');
 whitenoise = wgn(m,n,power)

noiseSignal = addnoise(signal,noise,snr);
 %FILTER THE NOISE
%filteredSignal =
 VSLMS(filter_size,ss_uper,ss_lower,input_file,iterations)

subplot(2,2,1);
plot(signal)
title('Input signal')

subplot(2,2,2);
plot(noiseSignal)
title('White noise and signal sum')

If it helps anyone, I'm using this addnoise function to add the noise.

From what I understand if I change power and SNR the white noise should look different based on what I input right? What am I missing ?

I've tried a bunch of powers and SNR combinations. SNR being 0.0004 or 0.9 or power being 30 or 0.0003 had no effect on the signal output.

So basically you have discrete samples of an audio signal previously recorded and you want to add noise to it. More specifically, you want to add white Gaussian noise (WGN) with a given power to obtain a noisy signal with a specific SNR. That's what the MATLAB function awgn do. You can use it this way:

noiseSignal  = awgn(signal,snr,'measured');

The input snr should be un dBW and the 'measured' config allows MATLAB to automatically measured the power of the input signal (signal should be a vector). The function generated then noise samples (complex or real depending on the signal input), with a power adapted to the SNR value specified, and add them to the signal samples. If you need the noise samples for a reason, you can simply do noise = noiseSignal - signal;.

For more details on the awgn function, you can look at the official documentation.

Furthermore, the noise is generated with pseudo random MATLAB functions, so if you want to always have the same noise samples you will need to have the "same" randomness each time you run the script. For example, you can do that by adding at the beginning of your script,

s = RandStream.create('mt19937ar','seed',529558);
prevStream = RandStream.setGlobalStream(s);

and restoring the default randomness by adding RandStream.setGlobalStream(prevStream); at the end of your script.

Finally, typical values for SNR in dB are integers between [-10;30]. A negative SNR just means that the power of the noise is higher than the power of the signal. The values you took for the SNR (0.004 or 0.9 dB) are really closed, and you will not be able to see the difference on the graphic (but if you look at the noise samples you should). Then, based on your code, changing the power variable is useless since the addnoise function seems to modify the noise power to match the SNR value.

Hope it helped.


How does subpixel image shifting using DFT really work? I am trying to assess the quality of several image interpolation methods for an application that involves generating subpixel-shifted images. I thought I could compare the results of a subpixel shift using all these interpolation variants with some perfectly shifted image, but it's probably not possible to obtain it (what would be the need for interpolation then?).

I was thinking about using DFT + shifting in the frequency domain, and I am unsure as to how it actually works in comparison to explicitly interpolating the image (using bilinear, bicubic, etc...). I'm sure it can't possibly generate a perfectly shifted image, but I can't put my finger on it. Is subpixel shifting with DFT equivalent to applying interpolation and if so, which one? What is the bias of the pixel values in images obtained using this method? Thanks!

EDIT:
After thinking the matter over, I figured since the FFT is an approximation (even more so the DFT) of the original function in terms of harmonics (sine functions), that it would amount to some sort of trigonometric interpolation. I recall a "Fourier series interpolation" formula for discrete data which was a trigonometric interpolation, but not sure if it's connected.

There are several key insights you need in order to understand how DFT allows you to shift an image.

First, Fourier's theorum: It's probably easier to look at the continuous (i.e., analog) case first.  Imagine you have some function, call it g(t).  For simplicity, let's say that g(t) is an analog audio recording, so it's a one-dimensional function, which is continuous, and represents the instantaneous pressure as a function of time.

Now, g(t) is one way we can represent our audio recording.  Another is G(f).  G(f) is the Fourier transform of g(t).  So, G(f) == FT(g(t)).   G(f) has all the same information as g(t), but it represents that information in the frequency domain instead of the time domain.  There are some nit-picky details about Fourier Transforms, which I won't mention.

You can sort of think of G(f) as the "distribution of frequencies" contained in g(t).  So, if g(t) is a sine-wave (i.e., a pure tone), then G(f) will be zero everywhere, except at the frequency of that tone.  This is probably a good point to mention that G(f) is in general a complex function--that is to say that it returns complex numbers, which can be thought of having a real and imaginary component or a magnitude and phase.

One small digression here:  Since g(t) is continuous (both in domain and range), G(f) is also continuous. So, how can G(f) be zero everywhere except for the tone frequency?  Well, FT(sin(wt)) = $\delta(w)$.  Where $\delta$ is the Dirac delta function.

Ok, so now we've got continuous FT's under our belt.

Here's the second insight:  A Discrete Fourier Transform is to a Fourier Transform as a sampled signal is to an analog signal.  In this case, the
"discrete" refers to quantization of the function's domain (time or frequency), not it's range.  (The sampled digital signal you get from your sound card is quantized in both domain and range.)

The digital byte-stream you get from your sound card contains "samples" of the original continuous (analog) signal from the microphone.  If we take the DFT of our sampled g(t), we still get a G(f).  G(f), remember, is just a different way of representing the information contained in g(t). If we obeyed Nyquist's theorum, the sampled signal g(t) contains all the "intelligence" of the original continuous signal, so our discrete G(f) must contain all the information from our original continuous signal. Parenthetically, G(f) is still a complex function.

This is where the magic of sub-pixel shifting comes in, but in this case I'm going to write about shifting the audio signal in time by less than a sample, since it's the same thing.

Remember how G(f) is a complex function?  It needs to be complex in order to represent frequencies that aren't zero at t=0.  Remember, sin(0)=0, so sin(2*0)=0, etc.  But what if we started recording a quarter of the way through a cycle of the tone?  That's where the phase part of G(f) comes from.  In this case, the phase would be 90 degrees or pi/2 radians, depending on your preference for representing a quarter of a cycle.  So G(tone_frequency) = 0 + i or $e^{{i \pi}\over{2}}$.

That means we can shift our audio recording in time (by any amount we choose, including a fraction of a sample time) simply by modifying the phase of G(t). Actually, that statement is perhaps a bit too casual.  For an un-quantized, sampled signal the phase can be adjusted arbitrarily (this is part of the reason I made the distinction between quantization of domain and range earlier).  However, for a quantized sampled signal (our byte-stream of audio, for example) the quantization step size (i.e., number of bits) determines the resolution with which we can adjust the phase.  When we Inverse Fourier Transform G(f) (or DIFT it, for this sampled signal), the new set of samples g'(t) = DIFT(G(F)) will all be shifted in time by the amount we pick.

Applying this to your pixels simply means using a 2-dimensional FT instead of the 1-dimensional FT discussed here.


How to compute beam pattern given a delay-and-sum filter in spectral domain? Given a two-mic system (the mic spacing $d$ is known), how can I compute the beam pattern (see the Figure below) over frequency and angle $g(f,\theta)$ for a pair of spectral filters ($H_0$ for channel 0, and $H_1$ for channel 1)?

Both $H_0$ and $H_1$ are complex vectors.

One of the method for calculating beam pattern could be, you simulate/capture audio from different directions on array and filter these captures with filters($H$'s) generated. the audio coming from direction of beam pattern will be passed as it is, but the audio coming from other directions will be attenuated. the amount of attenuation with respect to spatial location is nothing  but the beam pattern .
for implementation,
Take audio in time domain, window it, take n-point DFT of this windowed time domain audio after this do element wise multiplication with filter coefficients to get beamformed output.


Kalman Filter for estimating position with nonconstant velocity & acceleration I am trying to estimate the position & head direction of a rodent going through a 2D environment (a circular surface of 1m radius). Above his head is an overhead camera which records 4 LEDs mounted on his head. The position is calculated as the centroid of the LEDs positions and the head direction as the angle of the LEDs (see images).

The issue is that wires attached to his head occasionally occlude 1-4 of the LEDs and thus prevent/complicate measurement. Furthermore, the head direction estimation is particularly subject to lots of noise due to nonhorizontal head movements which add a 3rd dimension to the LED movements which are assumed to be parrelel to the surface (thus introducing a sort of parallax).

I thought a kalman filter would be a great way too smooth this data. However, in Matlab it seems that to implement this I would need to assume either constant acceleration or velocity which is not the case since the rodent is freely moving.

Is a standard Kalman filter still possible? Maybe a nonlinear (extended) KF? Are there any other better solutions?

Thanks
Roman

Rodent in a square environment. I included this image since here the overhead wires and LEDs are more easily visible

Rodent in the circular environment. Notice that 1 of the LEDs (green) is occluded by the wire:

So this is just the start of an answer. I'll have to keep updating it as I go.

The first attempt is to say that the quantities you are interested in are the location of the center of the four LEDs, and the roll, pitch, and yaw (rotation angles) of the LEDs.

That means your Kalman FIlter state will be:
$$
\mathbf{x}_k = \left[x_k\ y_k\ \alpha_k\ \beta_k\ \gamma_k  \right]^T
$$
where $\mathbf{x}$ is the state vector, $x_k$ is the $x$ location, $y_k$ is the $y$ location, $\alpha_k$ is the roll, $\beta_k$ is the pitch, and $\gamma_k$ is the yaw.

For now, I'm assuming that we are not interested in estimating velocity, $v_k$., or acceleration, $a_k$.

Let's assume that the state update equation is then:
$$
\mathbf{x}_{k+1} = f(\mathbf{x}_{k}) + v_k = \mathbf{x}_{k} + v_k
$$
where $f$ is the state update function, here $f(\mathbf{x}_{k}) = \mathbf{x}_{k}$ so $f$ is just multiplication by the identity, $\mathbf{I}$, $v_k$ is a $5 \times 1$ Gaussian noise vector with covariance matrix
$$
R = \left[
\begin{array}{ccccc}
\sigma^2_x & 0 & 0  & 0 & 0\\
0 & \sigma^2_y & 0 & 0  & 0\\
0 & 0 & \sigma^2_\alpha & 0 & 0\\
0  & 0 & 0 & \sigma^2_\beta & 0 \\
0 & 0  & 0 & 0 & \sigma^2_\gamma\\
\end{array}
\right]
$$
i.e. the changes in all variables are just Gaussian (bandlimited) white noise. Note that this may not give the best results, but let's start there.

Then, the noise measurements will be the $(x,y)$ locations of the four LEDs:
$$
\mathbf{z}_k = \left[x^1_k\ y^1_k\ x^2_k\ y^2_k\ x^3_k\ y^3_k\ x^4_k\ y^4_k  \right]^T
$$
where $x^i_k$ for $i=1,2,3,4$ is the $x$ location of the $i^\mbox{th}$ LED at time $k$ and $y^i_k$ is the $y$ location of the $i^\mbox{th}$ LED at time $k$.

Then the output equation would be:
$$
\mathbf{z}_{k} = h(\mathbf{x}_k, r) + w_k
$$
where $h$ is a function of the state and the radius of the circle that the LEDs sit on $r$.  I have to go offline for a bit, so will fill in $h$ based on these equations later. Here, $w_k$ is an $8 \times 1$ noise vector with diagonal covariance matrix $Q_k$.

The idea where would be to make $Q_k$ time varying so that, when an LED disappears, the entries in the $Q_k$ matrix associated with those measurements have extremely large variances (meaning we are very unsure where the LED is). That should make the update just (effectively) keep the previous value with a minor update based on the new state estimate from the locations of the LEDs we CAN see.

In any event, it looks like your filter is going to be an extended Kalman filter (EKF) rather than a straight Kalman filter if we pursue this route.

The other option would be to use $(x^i_k,y^i_k)$ as the state variables so there's a more direct relationship between them and the measurements. However, then I'm not sure how to apply the fact that the four LEDs are in a set relationship to each other.

This is taking way longer than I expected, so I thought I'd just let you know I am trying to work on it...this week (and last weekend) were a dead loss for non-critical stuff.

The code below mainly just attempts to model the original measurements. No real attempt is made to use the EKF to do the estimation. Just using the mean location of the LEDs is currently all that it does.

Picture of an example track is below.

I'll try and get the EKF version going, but we'll see how this weekend goes.

R Code Below
    #27668

# Plot the data generated
plot_track <- function (xx,zz)
{
  plot(xx[,1],xx[,2], type="l")

  for (t in 1:T)
  {
    points(zz[t,1], zz[t,2], col="red", pch=18)
    points(zz[t,3], zz[t,4], col="green", pch=18)
    points(zz[t,5], zz[t,6], col="yellow", pch=18)
    points(zz[t,7], zz[t,8], col="blue", pch=18)
  }

  points(xx[1,1],xx[1,2],col="magenta",lwd=20, pch=5)
  points(xx[T,1],xx[T,2],col="orange",lwd=20)
}

# From a single x-y track, create all the locations for all the LED location measurements.
# Note that this just takes account of YAW, not pitch and roll.
create_measurements <- function(x,y,radius,yaw)
{
  measurements <- c(x,y+radius,x,y-radius,x+radius,y,x-radius,y)
  dim(measurements) <- c(2,4)

  yaw <- c(cos(yaw), -sin(yaw), sin(yaw), cos(yaw))
  dim(yaw) <- c(2,2)

  measurements <- yaw %*% measurements

  # Reshape for later use in rbind
  dim(measurements) <- c(1,8)

  return(measurements)
}

# GIven the true (x,y) tracks in x, create the locations of the four LEDs
create_all_tracks <- function(tracks, radius)
{
   led_tracks <- data.frame(x1 = numeric(0), y1 = numeric(0), x2 = numeric(0), y2 = numeric(0), x3 = numeric(0), y3 = numeric(0), x4 = numeric(0), y4 = numeric(0))
  for (t in 1:T)
  {
    led_tracks[t,] <- create_measurements(tracks[t,1], tracks[t,2], radius, tracks[t,3]);
  }

  return(led_tracks)
}

# Generation of the original data, Don't do it if z exists.
#if (FALSE == exists("led_tracks_true"))
{
  radius <- 0.1 # Radius of each LED ring

  Nstate <- 3 # States are X,Y, RYaw Angle
  x0 <- rep(0,Nstate) #The state vector starts at all zeros

  sigmas <- c(0.1,0.1,0.001)
  Q <- diag(sigmas)

  T <- 1000
  x_true <- rep(0,Nstate*T) # The true states
  dim(x_true) <- c(T,Nstate)

  #A = [ 1 dt dt^2/2 0
  #      0  1     dt 0
  #      0  0      1 0
  #      0  0      0 1];
  dt <- 0.1
  A <- matrix(c(1,0,0,0,dt,1,0,0,dt^2/2,dt,1,0,0,0,0,1),4,4)
  B <- matrix(c(0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1),4,4)

  true_acceleration <- rnorm(T,0,1)

  for (t in 1:(T-1))
  {
    for (n in 1:Nstate)
    {
      x_true[t+1,n] <- x_true[t,n] + rnorm(1,0, sigmas[n])
    }
  }

  led_tracks_true <- create_all_tracks(x_true,radius)

  plot_track(x_true,led_tracks_true)
}

# Generate the measurements by adding some noise to the true led tracks
n_measurements <- 8
measurement_sigmas <- rep(0.1,n_measurements)
R <- diag(measurement_sigmas)

led_tracks_measured <-  rep(0,n_measurements*T)
dim(led_tracks_measured) <- c(T,n_measurements)

for (t in 1:(T-1))
{
  for (n in 1:n_measurements)
  {
    led_tracks_measured[t+1,n] <- led_tracks_true[t,n] + rnorm(1,0, measurement_sigmas[n])
  }
}

# Pick some times where we have an occluded LED
n_occlusions <- 50
occlusion_times <- sample(1:T, n_occlusions)
occlusion_leds <-  sample(1:4,T, replace = TRUE)

led_tracks_measured[occlusion_times, occlusion_leds] <- NA

# Attempt #1 : Just average the locations of the four LED tracks.
# When an LED is occluded, just ignore that measurement and average
# the remaining three LED locations.
est_avg <- rep(0,2*T)
dim(est_avg) <- c(T,2)

for (t in 1:T)
{
  est_avg[t,1] <- mean(led_tracks_measured[t,c(1,3,5,7)], na.rm = TRUE)
  est_avg[t,2] <- mean(led_tracks_measured[t,c(2,4,6,8)], na.rm = TRUE)
}

plot(est_avg - x_true[,1:2])

# Attempt #2 : An EKF approach?


Are all least square filters adaptive? Are least square filters, or filters that minimize error energy, the same as least mean square adaptive filters?

TL;DR: No, they are not necessarily the same.

Gory Details

Least squares is just an optimization technique. It is used in a variety of ways.

For filter design it is used to select that realizable filter $H_r(e^{j\omega})$ that most closely matches, in the least squares sense, the ideal required filter response $H_i(e^{j\omega})$:
$$
H_r(e^{j\omega}) = \arg \min \parallel H_r - H_i \parallel_2
$$
where $\parallel \cdot \parallel_2$ is the 2-norm or least-squares norm.

This sort of filter $H_r$ is not adaptive. That is, it doesn't change once it has been designed.

Adaptive filters may also use the least squares criterion, but in a different way: as part of the adaptation step.

Adaptive filters start off with initial filter coefficients $\vec{w}_o[0]$ and then use an update:
$$
\vec{w}_o[n] =\vec{w}_o[n-1] + \mu g[n-1]
$$
where $\mu$ is the step-size and $g$ is the gradient of the least squared error  surface in the direction of the minimum (from our current "location" of $\vec{w}_o[n-1]$).

Here, $g$ is determined by our error criterion: least squares. This means:
$$
\parallel \vec{w}_{\tt opt} - \vec{w}_o \parallel_2
$$
where $ \vec{w}_{\tt opt}$ is the unknown optimal (minimizing) solution.


Fundamental period of $x[kn]$ Let $x[n]$ be a periodic with fundamental period $N$ and $y[n] = x[kn]$ where $k \in \mathbb{N}$ and $k\ge2$. Is $y[n]$ periodic? What's the fundamental period of $y[n]$?

Here is my answer:
If $N = mk$ where $m \in \mathbb{N}$ then $y[n+m] = x[k(n + \frac{N}{k})] = x[kn+N] = x[kn] = y[n]$. So in this case the fundamental period is $\frac{N}{k}$. Otherwise simplify $\frac{N}{k}$ such that $\frac{N}{k} = \frac{m}{n}$ where $\gcd(m , n) = 1$. The fundamental period is $m$ since $y[n+m] = x[k(n+\frac{nN}{k}) ] = x[kn + nN] = x[kn] = y[n]$.
I don't know whether my solution is correct and rigorous enough. For instance how can we prove there is no smaller period?

You are correct. The following is what you have shown.
Let $x[n]$ be a periodic sequence with period N :
$$ x[n] ~=~ x[n + rN] ~~~~~,~~~\forall r,N \in Z^+ $$
Define $y[n] = x[M n]$, where $M$ is a positive integer; then
$$ y[n] = x[Mn] = x[Mn + rN] = x[M(n + \frac{rN}{M})] = y[n+K] $$
$y[n]$ will be periodic with $K$, if $\frac{r ~N}{M}$ is an integer. Hence, set $r$ to any value that makes $\frac{rN}{M}$ an integer.
If $N,M$ are already coprime, then setting $r=M$ makes $K=N$ as the period of $y[n]$. If they are not already coprime, then make them coprime $M_c,N_c$, and then set $r=M_c$ to make $K = N_c$ as the period of $y[n]$.
Make $N,M$ coprime by dividing them by their greatest common divisor $g = \text{gcd}(N,M)$ :
$$ N_c = N / g ~~~,~~~~ M_c = M / g $$
Then, for the general case, the period of $y[n]$ is:
$$ \boxed{ K = \frac{N}{ \text{gcd}(N,M) } } $$


What equation predicts the amplitudes of harmonics from a square/triangle/sawtooth/pulse oscillator? I have seen pictures like this which depict the shapes of amplitudes from the various common types of audio oscillators:

Similar pictures of spectra are shown here.
I am attempting recreating these with additive synthesis using sine waves, and I thus need to know what levels to set the component sine waves to.
I'm wondering if there are simple equations that can predict the level of a given partial/harmonic relative to a level of 1 for the fundamental for each type of wave.
I think that in theory these should all follow some sort of $y=\frac{1}{x^c}$ distribution, where $x$ is partial number. If so, what would $c$ be for each of these?
Assume that it is easy enough to just set the even harmonics to null when needed. I think the pulse will be hard to predict too. But I figure the active saw/triangle/square partials should have pretty simple governing relationships.
Any help is appreciated. Thanks.

Update:
In response to Matt's suggestion to look at the Fourier series, I did so, and I think I have the solutions for all but the pulse.
1) Triangle
Fourier series: https://mathworld.wolfram.com/FourierSeriesTriangleWave.html
Amplitude: $y=\frac{1}{n^2}$ for odd partials, even partials are zero
2) Square
Fourier series: https://mathworld.wolfram.com/FourierSeriesSquareWave.html
Amplitude: $y=\frac{1}{n}$ for odd partials, even partials are zero
3) Sawtooth
Fourier series: https://mathworld.wolfram.com/FourierSeriesSawtoothWave.html
Amplitude: $y=\frac{1}{n}$ for all partials
4) Pulse
Fourier series: https://lpsa.swarthmore.edu/Fourier/Series/ExFS.html#EvenPulse
This becomes complex once it's no longer a 50% duty cycle square and I'm not sure what that amplitude expression would be though I'd love to know.
Did I get that right and if so, any thoughts on the pulse?

This is what Fourier series are all about. Under relatively mild conditions, a $T$-periodic function $f(t)$ can be represented as an infinite sum of complex exponentials:
$$f(t)=\sum_{k=-\infty}^{\infty}c_ke^{j2\pi kt/T}\tag{1}$$
Note that if $f(t)$ is real-valued, $(1)$ can equivalently be written as a sum of real-valued sinusoids:
$$f(t)=c_0+2\sum_{k=1}^{\infty}|c_k|\cos\left(\frac{2\pi kt}{T}+\phi_k\right)\tag{2}$$
with $c_k=|c_k|e^{j\phi_k}$. From $(2)$ it should be clear that not only the amplitudes $2|c_k|$ of the individual sinusoids but also the phases $\phi_k$ are important.
The complex Fourier coefficients $c_k$ are given by
$$c_k=\frac{1}{T}\int_{0}^{T}f(t)e^{-j2\pi kt/T}dt\tag{3}$$


Strange step response In the figure there is represented a step response of a mechanical system. I created a step variation in the excitation of a DC generator. This lead to a step variation in the torque developed by that generator. I measured this torque with a arm connected to the stator and to a load cell.
Is this thing normal?

I tried to identify the transfer function of this system have as output the plot from the figure and a step change as input. I used IDENT from MATLAB and got a transfer function with 3 zeros and 3 poles with a fit of 85%. I used this transfer function with a sin wave input and computed the phase shift between the input and the output. This shift is very small, 2-3°. Is this thing normal?
I see that my system has a big delay from the step response but with sin wave response the delay is very small.
Is it possible that the procedure for obtaining the transfer function is not correct?

The step response you have is non-minimum phase. That means various things:

There are zeros (or poles, if the system is unstable) outside the unit circle (in discrete time).
The peak of the impulse response is delayed from the origin.

Below is an example showing an "ideal" delayed step, a non-minimum phase filter impulse response and the output.

As you can see, the output exhibits similar characteristics to what you have plotted.