How to prepare different input size in CNN I'm doing a research about personality identification based on their signature using CNN method, however the learning feature for the personality traits have a different input size. I understand that in CNN we have to make the same size for the input, but can I make the input size in different size?

Here is the journal that I use for my references for the learning feature.

You can research "Spatial Pyramid Pooling (SPP)"
Btw, as far as I know, you don't need make same size for input for convolutional layer but FCN layer


Fourier Transform of image convoluting with kernel edit: clarifying question. ...

The function
$$
\frac{\sin(x)}{x}
$$
has a well-known Fourier transform:
$$
\int_{-\infty}^{\infty}  \: \frac{\sin{x}}{x} e^{-j 2 \pi f x} dx = \begin{cases} \pi & |f| < \frac{1}{2 \pi} \\ 0 & |f| > \frac{1}{2 \pi} \\ \end{cases}
$$
which is derived here.

Your function:

$$
\frac{\sin(x)\sin(y)}{xy}
$$

is a simple separable two dimensional version of it. The 2D Fourier transform will just be:

$$
\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}  \: \frac{\sin{x}}{x} \frac{\sin{y}}{y} e^{-j 2 \pi f_1 x} e^{-j 2 \pi f_2 y} dx dy =
\int_{-\infty}^{\infty} \: \frac{\sin{x}}{x} e^{-j 2 \pi f_1 x}  dx \int_{-\infty}^{\infty}  \:\frac{\sin{y}}{y} e^{-j 2 \pi f_2 y}  dy =\begin{cases} \pi & |f_1| \mbox{ and } |f_2|  < \frac{1}{2 \pi} \\ 0 & |f_1| \mbox{ or }  |f_2| > \frac{1}{2 \pi} \\ \end{cases}
$$

Below is a plot that takes a random image with a white square in the middle (top left), FFTs it, multiplies the FFT by a brick filter (top right), inverse FFTs it and plots the result (bottom right).

R Code Below

# 26383

img <- runif(128*128, 0, 1)
dim(img) <- c(128, 128)

img[32:96,32:96] <- 1

N <- 50

imgFFT <- fft(img)
imgFFT[N:(128-N+2),] <- 0.0000000000000
imgFFT[,N:(128-N+2)] <- 0.0000000000000
img2 <- fft(imgFFT, inverse = TRUE) / 128 / 128

img2[Mod(img2) > 1] <- 1

mx <- max(log(Mod(imgFFT + 0.00001)))
mn <- min(log(Mod(imgFFT + 0.00001)))

scaled_log_imgFFT <- (log(Mod(imgFFT + 0.00001)) - mn)/(mx - mn)

par(pty="s") # make the plit square
plot(0:256, type='n')
rasterImage(as.raster(img),0,128,128,256)
rasterImage(as.raster(scaled_log_imgFFT),128,128,256,256)
rasterImage(as.raster(Mod(img2)),128,0,256,128)


How does google map put together satellite pictures? Just out of curiousity, does anyone know how does google map put together all those massive quantity of satellite pictures with with high accuracy? What advanced computer science technology has it employed?

This certainly is a signal processing question albeit a little vague. I can't say for certain that Google use it but the SIFT (shift invariant feature transform) algorithm http://www.cs.ubc.ca/~lowe/keypoints/ has been used in variant form in technologies such as Microsoft's Photosynth. It allows key features of multiple images containing overlapping content (supject to arbitrary X,Y,Z shifts) to be registered and matched. In simple terms...very clever image stitching.  The link provides all the paper references to the original algorithm. For another example of such technology check out this lecture about photosynth > http://www.ted.com/talks/blaise_aguera_y_arcas_demos_photosynth.html


how to evaluate derivative of convolution integral? If the signal $x(t)$ has ordinary first derivative $\dot x(t)$, then $\dfrac{d}{dt}\big(x(t)\star y(t)\big)$ is:

  (a) $\dot x(t)y(t)$
  (b) $x(t)\dot y(t)$
  (c) $\dot x(t)\star y(t)$
  (d) $\dot x(t)\star \dot y(t)$

I solved it as follows: we know,

$$I(t)=x(t)\star y(t)=\displaystyle\int_{-\infty}^{\infty}x(\tau).y(t-\tau) d{\tau}\implies$$
$$\begin{align}\dfrac{dI(t)}{dt}&=\dfrac{d}{dt}\displaystyle\int_{-\infty}^{\infty}x(\tau)y(t-\tau) d{\tau}\\&=\displaystyle\int_{-\infty}^{\infty}\dfrac{\partial}{\partial t}x(\tau).y(t-\tau) d{\tau} \\&=\displaystyle\int_{-\infty}^{\infty}x(\tau).\dfrac{\partial}{\partial t}y(t-\tau) d{\tau}=x(t)\star \dot y(t)\end{align}$$

So option $(b)$ should be correct but answer in the book is given as $(c)$ but I think both options $(b)$ and $(c)$ should be correct because in above
method I fixed $x(t)$ and shifted $y(t)$ but we can also do that other way
round fixing $y(t)$ and shifting $x(t)$,and we know by doing so,convolution
integral remain unaffected,so tell where i'm wrong ,any help would be appreciated

Note that option (b) is not correct, and that it is also not equal to what you came up with. Option (b) is just the multiplication of $x(t)$ and $y'(t)$, not the convolution. Your solution and option (c) are both correct, assuming that all derivatives exist and that the convolution integrals converge, because with that assumption the following holds:

$$\frac{d}{dt}\left(x(t)\star y(t)\right)=x'(t)\star y(t)=x(t)\star y'(t)\tag{1}$$

where $x'(t)$ and $y'(t)$ denote the derivatives of $x(t)$ and $y(t)$, respectively.

An ideal differentiator is just an LTI system, so it can be written as a convolution, and since convolution is associative, $(1)$ must hold.

This can also most easily be seen in the frequency domain. Differentiation corresponds to multiplication with $j\omega$, so we have

$$j\omega\left(X(\omega)Y(\omega)\right)=\left(j\omega X(\omega)\right)Y(\omega)=X(\omega)\left(j\omega Y(\omega)\right)\tag{2}$$

which is equivalent to $(1)$.


Is the inverse of a causal system also causal? If I have a causal system H(z) and I find the inverse of this system:

$$
G(z) = \frac{1}{H(z)}
$$

Is G(z) also causal?

It's not sufficient to only consider causality, you also need to check whether the inverse system is stable, otherwise it can't be implemented. If $G(z)$ has zeros on the unit circle, it cannot be inverted. If $G(z)$ has no zeros on the unit circle, but if there are zeros outside the unit circle, then there is no causal and stable inverse, because the zeros of $G(z)$ become the poles of the inverse system. There is, however, a stable inverse but it is not causal. If $G(z)$ has all its zeros (and poles) inside the unit circle, then there exists a causal and stable inverse (with all its poles and zeros inside the unit circle). Such a system is called strictly minimum-phase, and its inverse is also a strictly minimum-phase system.


Deciding on the FFT length for water waves in flume (wave facility) Soon I'll be conducting experiments in a wave facility (flume) for my MSc thesis. And I would like to analyse wave spectra (amplitude and energy density).

The time step (or 'accuracy') of the wave maker, is by default dt = 0.001 s.

The sampling frequency of the measurement instruments (wave gauges) is 100 Hz (previously 25 Hz).

I'll start with bichromatic waves, which can be more or less be treated as being deterministic. I would like to prevent spectral leakage as much as possible.

The fundamental frequency is different for different experiments: short waves with a wave period of around 2.5 s. The wave maker includes bound waves (subharmonics), with the wave period Tbound one over the difference frequency of the short waves in the bichromatic wave train.

So far, previous students have set the FFT length to a power of 2. From what I read this is not necessary with modern implementations of the FFT.

One example of a previous student:

Fs = 25 Hz,

N = D*Fs = 2^13 (power of 2)

This results in a duration of the experiment of
D = N/Fs = 327.68 s

with a spectral resolution of
df = 1/D = 0.00305175781 Hz.

To avoid spectral leakage, frequencies of all my wave components should be a multiple of this number if I am correct.

Since I've set my wave gauges at a sampling frequency 100 Hz, I'm wondering if it's fine to increase my duration to D = 1000 s, so N = 100*1000 and is not a power of 2.

The resulting spectral resolution is:
df = 1/D = 0.001 Hz
and I'll round the frequencies of all my wave components to a multiple of this number.

So, the best news first:

Nyquist is your friend! As long as your sampling rate is sufficiently high, all the information contained in a band-limited analog signal can be represented by the digitized signal.

In other words: you can't sample too fast, and even if the ratio of sampling rate and signal frequency is irrational, you don't lose information.

So, you can follow the sample first, ask questions later philosophy¹.

Then, yes, if you want your periodic signal to occupy singular bins in your DFT, the DFT length needs to be an integer multiple of your period. And yes, FFTs can nowadays be implemented with a lot less restricted lengths². Your $10^N = 2^N\cdot5^N$ certainly is covered by the FFTW.

I'd like to encourage you to not only consider the DFT as frequency estimator.

MUSIC is able to calculate the spectrum at arbitrary frequencies (ie. it's an algorithm that you can ask "hey, what's the amplitude at f=0.123141414288?") instead of giving you an frequency-uniform amplitude vector, and if you know that you're looking for a fixed number of "peaks" in the spectrum, ESPRIT might be your thing. The accuracy of these algorithms is not bound by an orthogonal basis of frequencies, but on your ability to estimate the signals autocorrelation – which you can increase by adding more observation.

So, maybe you'd want to do the following:

Apply the DFT to your signal to get a "rough idea" of where your energy is in frequency
use MUSIC to generate a (pseudo-)spectrum around these frequencies of interest with a very fine resolution
As soon as you know how many tones to expect in the spectrum, throw ESPRIT at the problem to find exactly that number of dominant frequencies in your signal

¹: I'll need to make "SDR: Sample first, ask questions later" t-shirts.

²: The FFTW documentation (as far as I know) doesn't really specify, but in a half-sentence claims that things work best for lengths that can be factorized into prime factors $\le7$; one of the papers quoted on the homepage somewhere mentioned you can also have up to two factors of 11 and one factor of 13 in there to get an optimized FFT. That info is from the back of my head, might be wrong or outdated, and even if totally wrong: an absolutely naive DFT done by a matrix-vector multiplication will not put an unbearable load on your PC.
So simply use the DFT lengths you need, unless the vectors you need to transforms are in the millions of values.
And, directly from the fftw.org home page: Arbitrary-size transforms. (Sizes with small prime factors are best, but FFTW uses O(N log N) algorithms even for prime sizes.)
So, assuming you're doing a $10^6$ transform, and FFTW couldn't do better than $\mathcal O(N\log N)$, that'd be in the order of millions of floating point operations per transform. Assuming that your PC is about 20 years old and can do let's say 10 Million of these per second – do you have to worry at all?


Peak Signal to Noise Ratio (PSNR) in Python for an Image How do I calculate the maximum signal to noise ratio (PSNR) in Python? Is there any library that can calculate PSNR for Image?

The steps for calculation of PSNR value of two images:

For details click here

import numpy
import math
import cv2
original = cv2.imread("original.png")
contrast = cv2.imread("photoshopped.png",1)
def psnr(img1, img2):
    mse = numpy.mean( (img1 - img2) ** 2 )
    if mse == 0:
    return 100
    PIXEL_MAX = 255.0
    return 20 * math.log10(PIXEL_MAX / math.sqrt(mse))

d=psnr(original,contrast)
print(d)


Why is not Fourier Transform Good for Non-linear Processes I was reading through slides about Hilbert Huang Transform. In slide 14, which talks about the motivations of a new method instead of Fourier Transform (FT), the author provides those two reasons in addition to other reasons

Physical processes are mostly nonstationary.
Physical processes are mostly nonlinear.

I understand the first reason, where the FT cannot locate the frequencies in time. However, I do not understand why the FT does not work for the nonlinear processes. Based on my understanding, the FT is a method to decompose the signal into different and simple components. Why does condition of non-linearity play a role here?

My first thought for reasoning that assumption, was that we cannot represent a non-linearity in the signal by a linear combinations of sinusoids. However, I think sinusoids themselves will handle the non-linearity.

Because complex exponentials $e^{\jmath \omega t}$, which are results of Fourier transform, are the eigenfunctions for linear, time invariant (LTI) systems. See eigenfunction of LTI. Also see this answer on SP.SE.

Thus, Fourier transform is useful for analyzing linear (not suitable for non-linear one) time invariant (can be intepreted as stationnary) system.


How to determine relative distance between 3 in-phase transmitters by analysis of received signal I would like to determine with high accuracy the relative distances between three RF transmitters a few meters apart by analysis of phase interference in the received wave, the receiver being at an undefined distance a few hundred meters away. Is there a simple method for this? The three RF transmitters can be arranged to transmit the same phase and frequency, or differing frequencies, whatever is required. Is there some kind of radar or reflection interferometry technique that exists and provides the information required, from reflectors?

Check this out, it does calculate the similarities of 2 or more signals
http://www.mathworks.com/help/signal/examples/measuring-signal-similarities.html


Understanding DFT-ODD operation in Gardner Efficient Convolution paper As a newbie to DSP I am trying to understand and implement an efficient convolution engine as per JAES V43.3 1995/03  - Gardner - Efficient Convolution without I/O Delay.

I'm pretty much there with a working proof of concept, but I am struggling to get my head around appendix two - which talks about DFT-ODD operations and computing spectra using the spectra of two existing blocks of results.

Can someone point me to some additional (readable for newbie) documentations that perhaps covers this topic. This is the last piece in my jigsaw puzzle.

The task is to compute the DFT of a sequence from the DFTs of its two halves:

$$X[k]=\text{DFT}_N\{x[n]\}=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}$$
$$\begin{align}x_1[n]&=x[n],\quad n=0,\ldots,N/2-1\\x_2[n]&=x[n],\quad n=N/2,\ldots,N-1\end{align}$$
$$X_1[k]=\text{DFT}_{N/2}\{x_1[n]\},\quad X_2[k]=\text{DFT}_{N/2}\{x_2[n]\}$$

i.e., compute $X[k]$ from $X_1[k]$ and $X_2[k]$. This exact question is treated in this answer, where it is shown (cf. Eq.$(2)$) that the even indices of $X[k]$ are easily computed as

$$X[2k]=X_1[k]+X_2[k]\tag{1}$$

but where it is also shown that there is no efficient way of computing the odd indices $X[2k+1]$ from $X_1[k]$ and $X_2[k]$, even though it is of course possible, albeit inefficiently.

For this reason, in his paper, Gardner suggests another method for computing $X[2k+1]$, called "DFT-ODD", which evaluates the DFT at odd indices. First, we need two equivalent representations of $X[2k+1]$:

$$\begin{align}1.\;X[2k+1]&=\text{DFT-ODD}_N\{x[n]\}=\text{DFT}_{N/2}\{(x_1[n]-x_2[n])W_N^n\}\\
2.\;X[2k+1]&=\text{DFT-ODD}_{N/2}\{x[2n]\}+W_N^{2k+1}\text{DFT-ODD}_{N/2}\{x[2n+1]\}\end{align}\tag{2}$$

where $\text{DFT-ODD}_N\{x[n]\}$ denotes the DFT of the length $N$ signal $x[n]$ evaluated at odd indices (i.e., its length is $N/2$), and $W_N=e^{-j2\pi/N}$. The first representation is identical to the one given in the above mentioned answer, and the second one follows immediately from the definition of the DFT, after splitting the signal in its even and odd samples.

We need one more expression relating the DFT-ODD of a complex-valued sequence of length $N/2$ to the DFT-ODDs of its real and imaginary parts, respectively:

$$\begin{align}g[n]&=g_R[n]+jg_I[n]\\
\text{DFT-ODD}_{N/2}\{g[n]\}&=G[2k+1]=G'[k]\\
\text{DFT-ODD}_{N/2}\{g_R[n]\}&=\frac12\left(G'[k]+G'^*[N/4-k-1]\right)\\
\text{DFT-ODD}_{N/2}\{g_I[n]\}&=\frac{1}{2j}\left(G'[k]-G'^*[N/4-k-1]\right)\end{align}\tag{3}$$

where $*$ denotes complex conjugation.

The DFT-ODD of the original sequence $x[n]$ is now computed as follows:

generate a complex sequence from the even and odd indices of $x[n]$: $g[n]=x[2n]+jx[2n+1]$
compute DFT-ODD of $g[n]$ using representation $1$ of Eq. $(2)$: $G'[k]=\text{DFT}_{N/4}\{(g[n]-g[n+N/4])W_{N/2}^n\}$
compute the DFT-ODDs of $x[2n]$ and $x[2n+1]$ from $G'[k]$ according to Eq. $(3)$: $\text{DFT_ODD}_{N/2}\{x[2n]\}=\frac12\left(G'[k]+G'^*[N/4-k-1]\right)\\\text{DFT_ODD}_{N/2}\{x[2n+1]\}=\frac{1}{2j}\left(G'[k]-G'^*[N/4-k-1]\right)$
compute $X[2k+1]$ from these two DFT-ODDs according to the second representation in Eq. $(2)$.

This little Matlab/Octave script shows how it works:

N = 256;    % some multiple of 4
x = rand(N,1);
N2 = N/2; N4 = N/4;
x1 = x(1:N2); x2 = x(N2+1:N);
X1 = fft(x1); X2 = fft(x2);
Xeven = X1 + X2;

xe = x(1:2:N); xo = x(2:2:N); g = xe + 1i*xo;
G = fft( ( g(1:N4) - g(N4+1:N2) ) .* exp( -1i*2*pi/N2*(0:N4-1)' ) );  % DFT_odd of g
Xe = .5 * ( G + conj( G(N4:-1:1) ) );                          % DFT_odd of x[2n]
Xo = -1i*.5 * ( G - conj( G(N4:-1:1) ) );                      % DFT_odd of x[2n+1]

Xodd = [Xe;Xe] + exp(-1i*2*pi/N*(1:2:N-1)') .* [Xo;Xo];

X = [Xeven.'; Xodd.']; X = X(:);

% check result
X2 = fft(x);
max(abs(X2-X))    % 2.8422e-14


Comparing the FFT to numerical integration in Matlab The result of calculating the fourier transform using numerical integration is:

the result of using Matlab's FFT is:

So where did I go wrong here? I know the FT of a Gaussian should be another Gaussian.  Shouldn't they be the same? Furthermore I would have thought the
Edit:
Following a suggestion using the fftshift function gets me to the following plot, but this looks very different to the original fft despite documentation stating that it just centres the result around zero frequency.

which is more resembling what I would expect.
omega = linspace(-30,30,1000);

t = linspace(-10,10,1000);
y = exp(-t.^2./2)/sqrt(2*pi);
Yft = fft(y);

% Pre-allocate results of integration.
Y = zeros(size(omega));

%Integrate.
for k=1:length(omega)
    Y(k) = trapz(t, exp(-t.^2./2)/sqrt(2*pi).*exp(1j*omega(k).*t));
end

figure()
plot(omega, real(Y), omega, imag(Y));

figure()
plot(omega, real(Yft), omega, imag(Yft));

There are two reasons why you get different results. The first is that you define a vector t between $-10$ and $10$, but the fft command doesn't know that. So the data are interpreted as starting at $0$. This means the first difference is that the FFT computes the transform of a delayed Gaussian impulse, not of an impulse centered at time zero. This is why you get a non-zero imaginary part in the result of the FFT.

The second problem is - as pointed out by nivag in a comment - that the first index of the result of the FFT corresponds to frequency zero, and the negative frequencies appear in the second half of the resulting vector.


Continuous frequency grid to discrete grid I am designing a controller for a system using several methods, one of which is a discrete method. In the continuous time I am looking at a logarithmically spaced frequency grid in the range $$\omega \in [0.1,1]$$Now I want to design a controller in the discrete time, but I want to evaluate the h-infinity norm (which I am minimizing) over the same frequency grid as in the continuous time.

How can I minimize the norm over the same frequency grid in discrete time?

If I understood that right, you are trying to minimize the infinity norm over a discrete set of points. If that is the case you set over which you are trying to minimize is not convex and this is very much a NP hard problem, unless the set is very small, so that we can rigorously loop over it. Those are my thoughts on this.


How to do continuous signal processing (i.e without windowing)? I'm working in speech recognition research, and I wonder if there is a way to analyse sound just like how ears do (i.e. without windowing), so for example to take some feature continuously with an update for each new sample.

It is similar to talk about windowing with overlapping ratio of 99% or something like this,  you can look at the question like: How does ear working without windowing concept?

Yes, it's possible to analyse sound the way ears do.

For example, you could compute the DFT of a signal continuously using several Goertzel filters.

$$
y_k[n] = e^{j2\pi k/N} y_k[n-1] + x[n]
$$

where $k= 0,1,\ldots, N-1$, so that $y_0$ is the DC or zero frequency term.

Of course, this is an unstable filter, so some resetting or forgetting factor is needed to keep it stable. There are ways to make it better behaved numerically (see pp84-93 of Stoer & Bulirsch).

The reason why we don't generally do things this way is that we are usually trying to do them in the most computationally efficient way possible.  And nature doesn't necessarily do this, so we find ways to do similar things but that work better on a computer.


Carrier frequency offset (CFO) estimation using cyclic prefix I understand the cyclic prefix, but in (1) it's said that: "If we consider that the channel effect is minimal and can be neglected. then the phase difference of the CP and the OFDM symbol which the victim of CFO can be found by $2\pi\beta N/N= 2\pi\beta$." and then, the author derives the equation for estimating the CFO $\beta$. My questions are:

How is this phase difference between the cyclic prefix (CP) and the corresponding rear part of an OFDM symbol obtained?
I'm confused about how the cyclic prefix is used to estimate the carrier frequency offset. Can you explain how the formula for $\beta$ is derieved?

Reference

(1)
W.aziz, E.Ahmed, G.Abbas, S.Saleem, Q.Islam, "Performance analysis of carrier frequency offset (CFO) in OFDM usign MATLAB", Journal of Engineering, Vol. 1,No.1, 2012

To understand the equation in question consider how frequency offset $f_\mathrm{off}$ affects the received signal $r(t)$. For simplicity the authors assume the noise-free case and a frequency flat channel, i.e. the received signal before frequency offset is equal to the transmitted signal $s(t)$. Frequency offset is a frequency shift and can be expressed in time domain by
$$
r(t) = s(t) \mathrm{e}^{j2\pi f_\mathrm{off}t}.
$$
The received signal is now sampled at $t = n/F_\mathrm{s},\, n\in\mathbb{Z}$, where $F_\mathrm{s}$ is the sampling frequency. This yields
$$
r(n/F_\mathrm{s}) = s(n/F_\mathrm{s}) \mathrm{e}^{j2\pi n f_\mathrm{off}/F_\mathrm{s}}.
$$
In the following we write $r_n$ and $s_n$ for $r(n/F_\mathrm{s})$ and $s(n/F_\mathrm{s})$, respectively. The authors introduce the normalized frequency offset $\beta = f_\mathrm{off}/\Delta f$ but they don't define what $\Delta f$ is. From the context I understand that it is the subcarrier spacing $\Delta f = F_\mathrm{s}/N$, where $N$ is the number of subcarriers. Using this relationship we can rewrite the above equation as a function of $\beta$:
$$
r_n = s_n \exp\left(j2\pi \frac{\beta}{N}n \right).
$$
Next we consider a received cyclic prefix sample at $n$ (see above) and the corresponding received payload data sample $N$ samples apart at $n+N$:
$$
r_{n+N} = s_{n+N} \exp\left(j2\pi \frac{\beta}{N}(n+N) \right)
= s_{n} \exp\left(j2\pi\beta + j2\pi \frac{\beta}{N}n \right)
$$
In the above equation we have exploited that the cyclic extension is just a copy of the last $N_\mathrm{g}$ samples of an OFDM symbol and therefore $s_{n+N} = s_{n}$ (for all $s_{n}$ that are in the cyclic prefix). We can now see that the phase difference between the cyclic prefix and its corresponding data samples is $\arg(r_n^* r_{n+N})=2\pi\beta$ which is the answer to your first question. It is now easy to obtain an estimate $\hat\beta$ of $\beta$
$$
\hat\beta = \frac{1}{2\pi}\arg(r_n^* r_{n+N})
$$
It's just an estimation because in reality the received signal is corrupted by noise and channel distortion. The estimation can be improved by not only comparing one cyclic prefix and data sample but all $N_\mathrm{g}$ cyclic prefix samples with their respective copies. This is done by adding all $r_n^* r_{n+N}$ in the cyclic prefix and answers your second question about how the formula for $\beta$ is derived.

I have to say that the paper you're referring to is poorly written and that I probably  wouldn't have understood it neither if I hadn't had some prior knowledge about this topic. Thus I suggest you also read the original work about this topic van de Beek et al.:

J. J. van de Beek, M. Sandell, and P. O. Borjesson, “ML estimation of time and frequency offset in OFDM systems,” IEEE Trans. Signal Process., vol. 45, no. 7, pp. 1800–1805, Jul. 1997.


Filter odd or even harmonics with notch or inverse notch filter Hi i had the following question.
I have a signal containing a 200Hz sine wave and it's odd and even harmonics (no other frequencys or disturbing signals are contained).
What i'm looking for is a kind of filter which is able to seperate the odd harmonics from the even ones.
The end goal is to be able to determine the amplitude and phase of each single of the odd harmonics as accurate as possible.
I know that you normally would use an FFT for this purpose but i'm trying to figure out if this wouldn't also be possible with a filter instead.
So i thought about two options:

Use a parallel circuit of notch filters (one for each harmonic to be filtered) to cut out the even harmonics i don't want to use and seperate the odd harmonics from eachother.

Do exactly the opposite and use "inverse notch filters" to filter out the odd harmonics i want to use and seperate them from eachother by using a parallel circuit of "inverse notch filters" (one for each harmonic i want to filter out and take amplitude and phase from).

I've added a little sketch as png to make things clearer.

The problem is that it's very difficult to find understandable information about notch filter or inverse notch filter. (Inverse notch filter is not the same as a peak filter)
For notch filters i only find spurious and short information like these:
https://zone.ni.com/reference/en-XX/help/371325F-01/lvdfdtconcepts/notch_peak_filters/
https://www.onesdr.com/2019/11/21/fm-notch-filters-why-you-need-one-with-most-sdrs/
https://www.youtube.com/watch?v=tpAA5eUb6eo
And the only source of information i found about inverse notch-filter was this paper i mentioned (on top of the second page of the paper):
https://www.scirp.org/pdf/_2013060615443504.pdf
The problem is i would need way more information about both filter types and how derive them from simpler filter types and calculate their transfer function and Z-transform.
So my question is (if i'm not completely on the wood path) which of the both approaches described you would prefer and can you give me some good sources of clear and understandable information about the two mentioned filter types and thier implementation in hardware ?
Our would you recommend me a completely different filter type for this task ?
Greetings

What you are looking for are what we, in the audio space, call comb filters.  Comb filters may or may not have a feedback path, just like FIR and IIR filters. In fact, there is a generalized theory about designing comb filters based simply on designing any digital filter.
In fact, you can think of every LTI digital filter as being a comb filter where the teeth of the comb are spaced apart by the sampling frequency.  To space the teeth apart by a frequency spacing much less than the sampling frequency, you need a delay element that is longer than a single-sample delay ($z^{-1}$).  If a delay with precision of an integer number of samples is good enough to tune your comb filter, good for you.  But if it isn't, you need to deal with interpolation and precision delay.


Power Spectral Density computation and units I want to make some calculs of power spectral densité of signal. For example a real voltage signal (physical unit : $V$) in time $g(t)$, its fourier transform $G(f)$ and $S_g(f)$.
   As far as I know, the power spectral density units for g(t) is $V^2/Hz$. However, I find in various sources the following equation to calculate the power spectral density :
  $S_g(f)=|G(f)|^2=G(f)G(-f)$

I do not understand this "definition" (of course the most fundamental definition of PSD is given by the Wienner-Kintchine theorem) , as it leads to a spectral density in $V^2/Hz^2$.

Regards,
      Mike

I don't know if I get your question, but here is something you can look into.

Fourier transform is basically converting a time domain function to frequency domain. That being said, Fourier transform basic formula gives you the amplitudes of various components of a given signal in frequency domain. Note the word amplitude. That's what Fourier transform provides you, amplitudes. These amplitudes are functions of, of course, frequency. By exponential Fourier transform formula, we get these amplitudes as complex quantities. That might be your $g(t) \iff G(f)$, where $g(t)$ gives the amplitudes of the signal vs. time, and $G(f)$ gives the amplitudes of the different components of the signal vs. their individual frequencies.

Now you are interested in PSD. That's power. So now you don't want amplitudes vs. frequency, you want power vs. frequency, right? Power is proportional to the square of amplitude. So how do you convert amplitude to power? You square the amplitudes.

But, amplitudes occur as complex quantities in exponential formula, remember? So if you multiply one complex quantity with its conjugate, you get the square of its amplitude. That's why you might use

$$S_g(f)=\lvert G(f)\rvert^2=G(f)G(-f)$$

Here, $G(-f)$ is the conjugate of the complex quantity $G(f)$. At $f = 0$, you of course get a real quantity. In that case, you have,

$$G(f)G(-f)=G(f)G(f)=\lvert G(f)\rvert^2$$

Hope that helps.


My preamble bits not enough for locking the timing recovery for a gaussian channel I have chosen the preamble to be 48 bits (16 bits for carrier recovery+16  bits for timing recovery+ 16 bits for AGC of 30dB + 16 bits  for phase ambiguity) As channel is gaussian and there is no much frequency offset i was thinking number of bits would be ok? but i am seeing sometimes they is timing error? I choose the preamble length based on some sample application? Is there any criteria to choose the preamble properly? Does it depend on snr, data rate ...

Timing recovery is implemented as a control loop (measure the timing error and adjust the sample offset) and therefore would have an acquisition time depending on the control loop parameters (loop BW). A reasonable estimate for acquisition time assuming the start of acquisition is in the linear capture range of the loop is the first order 10% to 90% rise time estimate of $t_r = 0.35/BW$ where $t_r$ is the estimated rise time from 10% to 90% and $BW$ is the loop BW in Hz. This is accurate for a first order loop and a reasonable estimate whenever there is a dominant pole as that would set the settling time.

If the Initial state is not within the linear operating range of the loop or if there is lower SNR, this can take much longer. The best approach is to model this in a time domain simulation and observe the timing error versus time with step response testing to evaluate convergence time/ behavior and validate that the loop BW is as predicted from a more detailed control loop model from which these parameters can be computed and help identify issues such as capture range and non-linear behavior. A simulation can also help identify issues related the SNR as well as noise enhancement in the loop. For this purpose the preamble would be treated as infinitely long until it is certain the loop can work in the desired length in all conditions (will this loop operate with carrier offsets or will that be addressed before timing recovery, etc)


What is the effect of the natural logarithm in the frequency domain? I have read that taking the natural logarithm of a signal $x$ in the time domain will introduce a spectral broadening and hence "up-sampling" process before the logarithm is required. I tried to test that by plotting the spectrum (for example: from $-B$ to $B$) for an arbitrary limited band signal from $-B/2$ to $B/2$ and taking the natural logarithm in the time domain but I didn't see any effect of the broadening in the frequency domain, I tried the same method to show the effect of up-sampling and down-sampling (decimation in time) and it was working.

First note that when you use a logarithmic function you shall avoid negative arguments if it's output should be real valued. Then consider the following relation:
$$ y[n] = \ln( 1 + x[n] )$$
where the input $x[n]$ is conditioned such that $ 0 < 1 + x[n] \leq2$ is maintained. Then a Maclaurin series expansion of this expression will be

$$ y[n] = \ln(1 + x[n]) = \sum_{k=1}^{\infty} \frac{(-1)^{k-1}}{k} x^k[n] .$$

For very small $x[n]$ the following approximation can also be taken:
$$ y[n] \approx x[n] - x^2[n]/2 $$

Applying DTFT on this last line yields:

$$ Y(e^{j\omega}) \approx X(e^{j\omega}) - \frac{1/2}{2\pi} \left( X(e^{j\omega}) \star X(e^{j\omega}) \right) $$

The spectrum, $Y(e^{j\omega})$, of the logarithm of $1+x[n]$ will include multiple convolutions of $X(e^{j\omega})$ with itself. The relation becomes exact (with added terms which were ignored) as $k$ approaches infinity, showing that the spectrum will be broadened.

Indeed for any nonzero bandwidth input $x[n]$, the exact output $y[n]$ would have infinite bandwidth (or equivalently unavoidable aliasing) and no finite upsampling is enough to properly represent it. However in practice as $|x[n]|<1$ is maintained, its high powers, $x^k$ will tend to go zero quickly and can be ignored compared to some suitable threshold...


Kalman Filter Estimate vs ACF Least Squares Estimate I am currently reading Chapter 5, Applications to the Gas Markets, in Stochastic Modelling of Electricity and Related Markets by Benth, Benth and Koekebakker, World Scientific, 2008.

In the subchapter 5.3.1 Kalman filtering they assume that the detrended and deseasonalized logarithmic spot price dynamics for gas follow the process $X_1(t) + X_2(t)$ with
$$dX_i(t) = -\alpha_i X_i(t) + \sigma_i dB_i(t)$$
for $i=1,2$ with two independent Brownian motions $B_1, B_2$.

They cite Barlow, Gusev and Lai, 2004, which have shown that the confidence bands on the mean reversion estimates are quite wide. For this reason, they are proposing to use the empirical ACF to estimate the mean reversion parameters, assuming stationarity.

The stationary ACF of $X_1(t) + X_2(t)$ is given by
$$ \rho(\tau) = \hat{\omega}_1 e^{-\alpha_1 \tau} + \hat{\omega}_2 e^{-\alpha_2 \tau}$$
with the lag $\tau$ and the coefficients
$$\hat{\omega}_i = \frac{\frac{\sigma_i^2}{2\alpha_i}}{\frac{\sigma_1^2}{2\alpha_1}+\frac{\sigma_2^2}{2\alpha_2}}$$
for $i=1,2$. Using least squares they arrive at the values $\hat{\omega}_1 = 0.73$, $\hat{\omega}_2 = 0.27$, $\alpha_1 = 0.02$, $\alpha_2 = 0.28$. Using
$$\sigma_i^2 = 2\alpha_i \times \text{Var}(X_1(t)+X_2(t))$$
for $i = 1,2$ and estimating $\text{Var}(X_1(t)+X_2(t))$ to be $0.11$ they get as a starting point for the Kalman filter $\sigma_1 = 0.057$ and $\sigma_2 = 0.129$, which then yields $\hat{\sigma}_1 = 0.065$ and $\hat{\sigma}_2 = 0.573$.

How can the difference between the Kalman filter estimates $\hat{\sigma}_1$ and $\hat{\sigma}_2$  and the implicit least squares ACF values for $\sigma_1$ and $\sigma_2$ be explained? Where does it come from and why would the Kalman filter estimate be the right one in contrast to the ACF least squares estimate? In general, is the difference large?

As noted on p. 11 of this site: http://www.pages.drexel.edu/~pyo22/mem640/lecture04-Estimation/mem640Lecture-Estimation.pdf  least squares is a "special" case of Kalman Filtering. Additional details can also be found on other web sites, such as: http://www.insidegnss.com/node/3445

Also, as noted in both links above - although least squares is a special case of Kalman Filtering - the two algorithms differ from each other.  Note that the 'weighting matrix' W reduces to the identity matrix in the least squares case.  So you might expect that the results would also differ when you process the data.

As to which one will be more accurate - it all depends on how closely the 'signal processing model' of your algorithm matches the process that generates your input data.  I've programmed Kalman Filters whose results were completely wrong because: 1) it started with a bad first guess, or 2) the state model being used was wrong, or 3) the input data was far too noisy.  I've also programmed 'hypothesis testing' arrays of Kalman Filters, where all filters operated on the same data, but each filter was set up to match a particular target motion (ie: zig right, zig left, go straight, etc.) - in that case, only the filter that closely predicted the actual input data was the correct one (ie: residuals were very small), and the residuals of all the others were way out of range.


Minimum Samples Per Chip? Is there a 'Nyquist Criteria' for the minimum samples per chip that are needed for a correlator? I've heard 2, but will 1 work? What are the implications of having less than the minimum (in the context of radar and communications)?

Actually much less than 1 sample/chip will work for DSSS. Doing so erodes the spreading gain and makes chip level detection impractical, but it works. The undersampled spreading sequence simply acts as a vector (subspace) for the undersampled signal.  Matched filtering is still matched filtering -- as long as you are still sampled well above the symbol (not chip) rate.


How to distinguish between FFT power leak or close frequencies This week I was having an interview with a Professor for Grad School. He asked a bunch of questions related to signal processing that lead up to some questions about power leak on the FFT spectrum. So far so good, but then he caught me when he asked something around those lines:
"given just an FFT spectrum, with no prior information about the signal or the sampling frequency, how one can differentiate if two bins close together are a case of power leakage or two distinct frequencies?"
At the time of the interview nothing came to mind so I just answered that it was not possible but since then I was wondering if there is a way of making this distinction. I tried looking online but couldn't find any material that answered this so any help would be greatly appreciated.
Thanks!

Assuming the spectrum is sufficiently interpolated to see the detail within the resolution bandwidth (which is a reasonable assumption given the question), and the waveform isn't further windowed beyond the rectangular window (as I would argue that any such windowing that is applied is equivalent to the presence of multiple tones closely spaced given windowing is an applied amplitude modulation) then one obvious way to quickly distinguish these two cases is by observing the width of the main lobe (as well as all other features such as null locations, or non-existence of the nulls, but the main lobe width would be robust under noise conditions)-- given the kernel of the rectangular window (the DFT of the window ) convolves in frequency with each tone. If there is one tone we would see the kernel alone. If there are two tones spaced closer to together, then both tones would convolve with the kernel, and therefore widen the result.  (The kernel for a rectangular window in time is the Dirichlet Kernel in frequency, basically an aliased Sinc function, so if the spectral peaks are the same width as the Dirichlet Kernel, which is known, not computed, then we can easily declare that peak to be a single tone).
If sufficiently interpolated and noise free then we can also quickly detect this by the lack of frequency nulls (although the nulls are soon buried by any noise so this would not be a robust test). If the spectrum wasn't interpolated and the frequencies are so close that a bin width increase is not obvious, then we could compare the phase and magnitudes of the leakage values in all other bins to what would be expected with the underlying Dirichlet Kernel, under noise conditions we could use all bins and a correlation to what would be expected for a single tone to provide a confidence of one tone or not.
Further Details on using Main Lobe Width Measurement
This process of main lobe width measurement would be precise in confirming a single tone down to very small frequency offsets in the case of a single exponential tones of the form $e^{j\omega t}$, but will have a lower offset limit for the case of real tones due to the interaction of the two exponential tones that each real tone contains ($cos(\omega t) = 0.5(e^{j\omega t} + e^{-j\omega t}))$. In which case we would see even for a single tone a small variation in main lobe width versus frequency that becomes more pronounced as we approach the Nyquist boundaries. The effects of this, and interaction from any other tones further away for the same reason, become significantly reduced if the waveform was windowed prior to computing the DFT (which can't be assumed based on what the OP has stated, and as I already mentioned I would argue the result of windowing a single tone results in the presence of multiple closely spaced tones!). To show the practical extent of using the main lobe width approach, I plot the variability of the main lobe width versus frequency for the rectangular window case at the bottom of the post, which provides insight into the practical limitation of how close we can distinguish multiple tones (pretty close! down to 0.02 bins for most of the frequency range, and 0.2 bins for all frequencies except the end points).
To demonstrate this effect of simply measuring the width of the main lobe for the case of an interpolated spectrum, and its relative insensitivity to noise to address hotpaws comment, see this example below with a rectangular window and two equally leveled tones that are within one bin width (the resolution bandwidth of the rectangular window is 1 bin, so more than that would not be sufficient for "close"):

Here is the comparison of the composite signal with and without noise, where the noise is only 6 dB lower:

(And in this case, even if the spectrum was not interpolated, meaning the only samples available are at the integer bin locations; we can see how the increased width of the main lobe would still be obvious through inspection)
For further details on the limitations of the main-lobe width measurement for a real tone, the following plot below shows the variation of the main-lobe bin width at -3 dB for a rectangular windowed 128 point DFT (interpolated by zero-padding out to 12800 points), where the "truth" given by a single exponential tone approaches 0.886 bins for large $N$ (where $N$ is the number of samples, derived from $2\omega/\pi$ for $Sinc(\omega) = \sin(\omega)/\omega = 0.707$, and can be exactly determined for any size N from the Dirichlet Kernel, $D(\omega)$, from $2\omega/\pi$ for $D(\omega) = 0.707$, where:
$$D(\omega) = \frac{\sin(N\omega/2)}{N\sin(\omega/2)}$$
and $\omega$ is the normalized radian frequency in units of radians/sample : $\omega \in [0, 2\pi)$ for $N \in [0, N)$
With the results shown in the plot immediately below:

As given by the convolution, this bin width will be increased by the frequency separation of two tones that are less than 1 bin apart, and as long as that separation is greater than twice the variability of a single tone, it would be detectable (for example the variability when we are more than 10 bins away from the Nyquist boundaries is less than 0.01 bins, which holds as N is increased). This plot below was simulated but could be analytically predicted by the interference pattern of the underlying Sinc function that becomes the Dirichlet Kernel through aliasing. This demonstrates the practicality of using this for distinguishing a single tone from multiple tones closely spaced.


Scanning from Book - Gradient Removal This sample image:

is warped on the left side.

I don't want to dewarp, but just to remove shadow gradient (caused by book's spine). I guess gradient can be calculated from upper part of the image and that information can be used for reconstruction.

First I thought to apply some layer transformation in Photoshop. I took upper part with clear gradient and scaled it vertically to cover whole image in new layer. Then I did layer difference:

but while upper part is acceptably reconstructed, colored part is not.

Does anyone have an idea how to approach this in Photoshop, or Python/SciPy/OpenCV, or Matlab?

The simplest approach would be to divide the gradient rather than subtract it. Here's what the result looks like:


Augmentation for EEG signal classification using Deep Learning Augmentation is a technique that we use in deep learning for expanding the training dataset. It includes different ways of modifying an image and adding it to the training dataset. My question is if we are using augmentation for EEG classification and say we generate spectrograms and feed them as images to a CNN can we use Augmentation to increase the training dataset?

...if we are using augmentation for EEG classification and say we generate spectrograms and feed them as images to a CNN

If you are approaching this as an image based problem, it could possibly work. There are not enough details about the problem to provide a more accurate answer on this.

In any case, EEG records brain activity typically as a set of signals. Later, each one of those signals would lead to a time-frequency spectrogram representation. This means, one "image" per EEG channel which does not seem very practical from an augmentation point of view.

  ...can we use Augmentation to increase the training dataset?

A better approach to "augmentation" would be to generate surrogate data via suitable models.

This is closer to "augmentation" but still does not solve all of the problems you might be facing.

Jansen's model (for example) will generate realistic EEG waveforms with similar dynamics but only on a per channel basis (out of the box). If you are interested in replicating an EEG recording as a whole (that is, not only the dynamics of each channel but also their correlations), you would need to couple instances of Jansen's model and at that point you are still left with the question of "which coefficients are supposed to be associated with this task"?

So, if you are interested in creating a classifier that discriminates between "This is most likely an EEG signal rather than an electrode just dangling in the air", augmenting your dataset will work as expected. But if you are interested in creating a classifier that discriminates between "The subject is now more likely day dreaming rather than sleeping", then you need either very complex models or a large amount of data (under a similar condition) to infer some kind of a rule to exploit it later.

Hope this helps.


Generate projection image by OpenGL or Matlab I want to generate the "projection" image described as bellow. The left picture is the original picture and the right image is what I want. But, I have no idea about how to generate the right image. I know that OpenCV can only generate affine effect, so I know that only OpenGL or Matlab can generate this effect. However, I have no experience about OpenGL combined with Visual C++, so I don't know how to generate the right picture by OpenGL. Besides, according to other webside, I know that function "affine3d" in Matlab can achieve this effect, but my Matlab doesn't provide this function. Is anyone who can help me? Please. Thank you.

You would not need OpenGL for such perspective transforms.

First create a 3D cylinder with Coca-Cola texture on it. (This is found even in 3D warehouse).

Then, start by creating a pseudo-camera matrix $\mathbf{K}$ and then generate arbitrary rotation angles with desired range for each axis. Then compose the final rotation by $\mathbf{R}=\mathbf{R}_x\mathbf{R}_y\mathbf{R}_z$. Also synthesize a translation $\mathbf{t}$. Your projection matrix then reads $\mathbf{P}=\mathbf{K}[\mathbf{R} | \mathbf{t}]$. Project your 3D points with this matrix and you will get the warping effect.


Pulse Doppler Radar: how to filter its Range-Doppler matrix? I want to simulate a Pulse Doppler Radar with medium pulse repetition frequency (PRF about $100\textrm{ kHz}$, about $500$ pulses in pulse train).

It's rather powerful, and should discriminate targets with $1\textrm{ m}^2$ RCS in ranges up to $200\textrm{ km}$. Okay, let's build its range-speed image (about $\pm 750\textrm{ m/s}$ in speed dimension, about $1500\textrm{ m}$ in range).

But there is a problem: when there are some targets nearby (about $10\textrm{ km}$) they give not just a peak that can exceed the discrimination threshold, but a cross-like figure (sidelobes in range dimension, end, especially, sidelobes in speed dimension). And if the target is close enough, freaky sidelobes still can exceed the threshold and be treated as target candidates.

Even if I use window function with heaviest sidelobe-cancelling effect (Blackman-Harris), the sidelobes became smaller bus weren't wided out completely (and still give me ghost targets).

I think that I missed some important step in Range-Doppler matrix processing. Maybe some additional filtering. Cannot you tell me, what shoud I do to get rid of the sidelobes?

Here's some illustration: on the top row there are Range-Doppler matrices for different PRF ratios, and on the bottom row there is the range-speed matrix slice by speed (left) and range(right) in the target peak. The speed dimension annoyes me most = )

P.S. The number of DFT bins: 540 by speed, about 1200 by range (1200, 1280 and 1360, because PRFs varies).

About the size of the sidelobes -- I don't know how to add the colorbar to such a multi-plot figure, but I've added two additional plots to the bottom. In the left-bottom corner there is the signal power with different speeds (with the fixed range of the max peak range) and in the right-bottom corner there is the signal power with different ranges (with the fixed speed of the max peak speed). The units are power dB.

Two suggestions that may independently, or combined, help to solve your problem.

Consider using a different window/weighting function for your FFT processing.  There are window functions with low uniform sidelobes such as the Chebyshev window. This window function allows for selection of a constant sideline level response.
It sounds like you are using a fixed detection threshold, consider using an adaptive detection threshold so that the threshold can vary with the combination of background noise and sidelobe artifacts.  Mathworks has a pretty good description for a simple case at Constant False Alarm Rate (CFAR) Detection.


Correct magnitude spectra of a cosine DFT? I've just started my course on DSP and haven't laid my hands on MATLAB yet. I was wondering if the plot of the magnitude spectra was correct for the below shown $x(n)$:

you're missing the frequency component at $k=N-k_0$.


How does the intuitive notion of causality fit in with control systems? Edit: By causality, in this question, I do not mean the traditional mathematical definition in the theory o signals and systems; I mean causality as in an intuitive 'what's moving/pushing what notion'...
I'm having trouble understanding how causality fits into our intuitive ideas of inputting our output back into the system as feedback. What does it mean physically that y, the output, is dependent on the previous y values?
I understand the control/feedback engineering intuition, but once I get into the math, what I see are neat differential equations that govern how every single component should vary under such constraints. But, I don't see how we get to talk about how a certain component's behavior is caused by others.
Take exponential growth, for example, the output's behavior is caused by the behaviors of its rate of change and of the input force. Is this causality statement really the case, or is it just nature doing its thing, or am I just mixing ideas and being way too pedantic?
Either way, I want to understand this concept to be able to fully grasp what it really means when, for example, a servomotor gives feedback information about its rotation, and how this relates to the mathematical/physical differential equation framework.

As near as we can tell by experiment, causality is nature's way of doing its thing.
Causality says that if you have a system $y(t) = h\left(x(t), t\right)$, and it is causal, then $y(t_0)$ is dependent only on values of $x(t)$ for $t < t_0$*.
Causality doesn't say you can't combine systems into larger systems (as a rather pertinent example, you can certainly combine systems in a way that makes a control loop).
Causality doesn't imply stability, or instability.  There's a rich set of behaviors that causality allows -- it just doesn't allow the output to follow future inputs.
So in your exponential growth example, if you have a system described by $\dot x = a x + u$, and $a > 0$, the output will tend to grow exponentially (unless $u$ magically or through feedback damps the tendency).  We can say that the growth is caused by $x$ feeding back on itself -- but that causation still happens without violating causality.  Any change in $x$ or $u$ at time $t_0$ will affect the derivative of $x$ at time $t_0$, but (in the absence of infinities) it won't change the immediate value of $x$ -- so, the system is still causal.
* I may be overly strict here -- I'm ruling out $x = u$ as a causal system.  I'm not sure what the popular definition is in control system terms, but I'm lazy so I'm going to let it stand.


Even and odd signals We've studied in signal and systems about odd and even signals. But what is the actual implementation or application of this property?
like for energy and power signals we know If you're using a signal to communicate, then the more energy it has, the farther it will reach (I'm oversimplifying a bit, but in general this is true). The more energy a signal has when it is received, the easier it is for the receiver to recover the information in the signal without (or with few) errors.

  In similar fashion whats with odd and even signal?

I am not sure that you would call this a practical application as much as a practical implication, but knowing (or forcing) a signal to be odd or even enables you to process it in a certain way.

For example, if you know (or force) a signal to be even, then the imaginary part of its Discrete Fourier Transform (DFT) goes to zero and you can represent this signal just as a sum of $\cos(\cdot)$ functions. This relates the DFT to the Discrete Cosine Transform (DCT).

Specifically:

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

This is a compact way of saying:

$$X[k] = \sum_{n=0}^{N-1} x[n] \cdot \left ( cos\left ( \frac{2 \pi k n}{N} \right ) - i \sin\left ( \frac{2 \pi k n}{N} \right ) \right ) $$

This makes it clear that $X \in \mathbb{C}$ (omplex) and is composed of a real part (where $x[n]$ is multiplied by the $\cos(\cdot)$) and an imaginary part (where $x[n]$ is multiplied by the $\sin(\cdot)$).

But, $\cos(\cdot), \sin(\cdot)$ are even and odd functions respectively and the way their symmetry works when combined with integration means that if a signal is an even periodic function, then it can be represented by a series of cosines.

OK, but I only have a sound recording that could be anything...how do I know (or force) it to be an even function? (or an odd function).

One of the key assumptions of the Discrete Fourier Transform is that $x[n]$ is a periodic function. That is, even if you pass some random recording of a person counting "one two three four", all that the DFT sees at its input is "one two three four one two three four one two three four one two three four...". In other words, once you reach the length of $x[n]$ (which is $N$), you wrap back to its start.

So, how do we make this an even signal?

Notice that in even signals $f(x) = f(-x)$. We could therefore generate a new signal that is double the length of $x[n]$ and once we reach the end of $x[n]$ we repeat it by playing it backwards. In this way, $x[N+n] = x[N-n]$ or "one two three fourouf eerht owt eno one two three fourouf eerht owt eno one two three fourouf eerht owt eno..." and so on.

In terms of a "demo", the DFT of [1,2,3,4] has an imaginary component while the DFT of [1,2,3,4,3,2] doesn't.

I do not think that you will find a direct application as such, it will probably be something that "exploits" the implication that a function is even or odd.

Hope this helps.


How to generate such sounds as in example? I found these files on Internet:

1.wav
2.wav
3.wav
4.wav

How to generate such sounds?

Geno Chen wrote "Acoustic Grand Piano (Instrument #0)"

What to do next?

https://ccrma.stanford.edu/

There is a Center focused on Music and Acoustics at Stanford University and their website has many resources. There are several papers on computer generation of musical instrument sounds.


Performing DFT of streaming audio problem. Is there a limit? I am trying to write software that will perform the discrete fourier transform of real time data coming from the microphone into the sound card on a computer. I am using Java with the javax.sound APIs.

I am capturing 250 ms of data each time. That means, as the sampling rate is 44100, I am buffering 11025 of 16-bit shorts and performing the DFT on this data in each iteration.

My DFT code correlates this window of data with sine and cosine waves from 19 Hz up to 198 Hz in steps of 1, so as this window is 250 ms of data , this means it is actually correlating the incoming audio with frequencies from 76 Hz up to 792Hz in steps of 4.

If I now wish to take a smaller window of data each time of 125ms instead of 250ms in order to speed things up. If i correlate this data with sine and cosine waves of incrementing frequencies by 1 like before, this means that it is actually correlating the incoming audio with waves in steps of 8Hz. This is not good for my problem as I want to be able to look for frequency differences in the signal of less than 8Hz.

My question is, is this the normal trade-off, speed versus frequency resolution , or is there some way around this ?

Many thanks in advance. I am new to audio processing so forgive me if i'm missing out on some common sense :)

The window size of your FFT limits your resolution and there is no way round this (i.e. you can't zero pad or zero stuff to get round it). The easiest way to think about it is: lower frequencies have longer periods so need a longer window (more samples) to differentiate.

I believe it would be impossible to increase your resolution even if you used an alternaitve method like a series of very steep band passes narrower than the frequency resolution permitted by FFT for your window size. Yes, you would get some measures out from your band passes but if the window size was not long enough, the data would be highly inaccurate in differentiating very close lower frequencies as there would not be enough data to disentangle them.

Consider a window of 1 sample or 2 samples, in these cases it is easy to see that you just can't recover resolution from small windows without making them bigger by gathering more independent samples.


Why Successive Approximation, Analog to Digital converter is most commonly used? There are various types of ADC in used, like,

ramp-compare ADC
Wilkinson ADC
pipeline ADC
Flash ADCs
Sub-Ranging ADCs
Folding ADCs
Pipelined ADCs
Successive Approximation (Algorithmic) ADCs
Integrating (serial) ADCs

but  Successive Approximation (Algorithmic) ADCs (ADC0808) is most commonly used.
I do not understand the reason as why it is better than any other ADC ?

According to this comparison,

There is absolutely nothing so special about it as compared to other, then why it is most commonly used for education purpose and small electronic projects etc ?

For many applications, the technology often selected is the lowest cost one that meets the most common requirements (audio sample rates, etc.)


Question on the update step of Bayesian filter The update step in Bayesian filtering is given as follows, in a textbook

$p(x_k|y_{1:k}) = \frac{1}{Z_k}p(y_k | x_k)p(x_k|y_{1:k-1})$,

$Z_k$ is defined as $\int p(y_k|x_{k})p(x_k|y_{1:k-1})dx_k $

I was able to derive the second term on the right .The steps are given below
$p(x_k|y_{1:k}) = \frac{p(x_k , y_{1:k})}{p( y_{1:k})}=\frac{p(y_{1:k}|x_k)p(x_k)}{p( y_{1:k})} $

$=\frac{p(y_k,y_{1:k-1}|x_k)p(x_k)}{p( y_{1:k})}$

$=\frac{p(y_k|y_{1:k-1},x_k)p(y_{1:k-1}|x_k)p(x_k)}{p( y_{1:k})}$

By Markov property and Bayes rule

$=\frac{p(y_k|x_k)p(x_k|y_{1:k-1})p(y_{1:k-1})}{p( y_{1:k})}$
$=\frac{p(y_k|x_k)p(x_k|y_{1:k-1})p(y_{1:k-1})}{p(y_k, y_{1:k-1})}$

Applying bayes rule to the denominator

$=\frac{p(y_k|x_k)p(x_k|y_{1:k-1})p(y_{1:k-1})}{p(y_k| y_{1:k-1})p(y_{1:k-1})}$

$=\frac{p(y_k|x_k)p(x_k|y_{1:k-1})}{p(y_k| y_{1:k-1})}$

Can someone explain how the denominator in the above equation can be written as

$Z_k$ = $\int p(y_k|x_{k})p(x_k|y_{1:k-1})dx_k $

I would prefer to see the definitions, but based on your prior work (which I haven't checked thoroughly), we can already reach the final stage. You left off at $$ p(x_k | y_{1:k}) = \frac{p(y_k|x_k) p(x_k|y_{1:k-1})}{p(y_k|y_{1:k-1})} $$ which is a probability distribution of $x_k$ based on the left hand side. The term in the denominator is independent of $x_k$ (if I have guessed your model correctly). Therefore, it is an unknown positive constant. That is, there exists some constant $Z_k$ such that $$ p(x_k | y_{1:k}) = \frac{1}{Z_k} p(y_k|x_k) p(x_k|y_{1:k-1})$$ where $Z_k = p(y_k|y_{1:k-1})$.

Now let's integrate both sides.
$$ \int p(x_k | y_{1:k}) d x_k = \int \frac{1}{Z_k} p(y_k|x_k) p(x_k|y_{1:k-1}) d x_k $$
The left hand side is a probability distribution with respect to $x_k$, therefore it sums to $1$. Because $Z_k$ is independent of $x_k$, we can pull it out of the integral. Multiply both sides by $Z_k$ to obtain your answer,
$$Z_k = \int p(y_k|x_k) p(x_k | y_{1:k-1}) d x_k.$$

Is $Z_k$ a constant with respect to $x_k$?

We know that $Z_k = p(y_k|y_{1:k-1})$ in other words $Z_k$ is the probability of $y_k$ given all previous $y_{1:k-1}$ after integrating out $x_k$. Note that this is the same thing as $$Z_k = \int p(y_k|y_{1:k-1},x_k) p(x_k) d x_k$$
which (after integrating) removes all dependency on $x_k$. Therefore, as noted in the original response, $Z_k$ is constant with respect to $x_k$ and can be factored out of the integral.

This practice of marginalizing out a variable based on its prior is a common theme in Bayesian systems.


Delay in Savitzky-Golay filtering I am applying a Savitzky-Golay filter to a signal, using the scipy function.

I need to calculate the lag of the filtered signal, and how much is it behind the original signal.

Could someone shed some light on this matter? How could I calculate it with scipy? How should I interpret the result correctly?

I would be very grateful!

Standard Savitzky-Golay filters are linear phase (type I) FIR filters. So they have an odd number of filter coefficients $2N+1$, and the delay equals $N$.

For a good overview of Savitzky-Golay filters see this article by Ronald Schafer.

For the definition of the four types of linear phase FIR filters see this answer.


Which domain used for denoising additive and multiplicative noises I want to know what is the difference between additive noise, multiplicative noise.. In what domain these noises are handled to remove?

I am interested particularly in image dataset

Does different domains are used in case of satellite imaging where noise due to reflection of frequencies is statistically independent?

I want to know what is the difference between additive noise,
  multiplicative noise.. In what domain these noises are handled to
  remove?

Additive noise and multiplicative noise are just models of how noise corrupts our data.

One very common model is the additive noise model, where we have our 'true' data vector $s[n]$, (which we are trying to ascertain), being corrupted by a noise vector, $v[n]$. What we are given, is $x[n]$, where:

$$
x[n] = s[n] + v[n]
$$

This is called the 'additive' noise model, because as you can see, noise is added to our true signal, giving us $x[n]$. There are many ways in which we can remove this noise in an additive model, such as filtering, (which is a form of averaging if you like). This is a very common type of noise model. If I am talking, and you are talking over me, your voice can be modeled as additive to my voice. You voice would be additive 'noise', that is added to my voice, which in this vain example is the 'true' signal, (although this would be disputed in a heated argument between two people). In a more objective example, thermal noise from a microphone's electronics can also be modeled as an additive noise, which is added to a voice signal that it has received. Many things can be modeled as additive types of noises.

Multiplicative noise on the other hand is still a model, but in this model our true data samples are being multiplied by noise samples, like so:

$$
x[n] = s[n]v[n]
$$

One common way to remove multiplicative noise is to actually transform it into an additive model, and then apply everything we know from the additive noise reduction field. We can do this easily via taking the logarithm of the signal, filtering, then inverse log transform. Thus we can do:

$$
x[n] = s[n]v[n] \implies \log(x[n]) = \log(s[n]v[n]) = \log(s[n]) + \log(v[n])
$$

At this point we have an additive model once more. Now, we can filter as we normally might, to remove or reduce $\log(v[n])$, and then simply take $\log^{-1}$ of the result, yielding an estimate of $s[n]$, our true signal.

  I am interested particularly in image dataset

One example of multiplicative noise is in illumination differences between images, which is solved in the above way. The non-uniform illumination across an image can be modeled as a pixel-by-pixel multiplication of the image by an illumination mask. This is also known as homomorphic filtering by the way. Anytime you can model a corrupting phenomenon as multiplying your clean data signal, you can use this multiplicative model.


GCC-PHAT for Multiple Sound Sources Can GCC-PHAT be used to detect multiple time-delays (that is, having multiple acoustic sources simultaneously)? The documentation for the gccphat(sig,refsig) function says that it assumes both sig and refsig came from a single source.

If  it is the case that gccphat(sig,refsig) can indeed only take into account a single source, can you recommend other approaches I could try to detect multiple sources? I am fairly new to this field so I could appreciate all the help. Thanks.

The theory behind the Generalised Cross-Correlation - GCC - (in general, not just the PHAT version of it) uses the cross-correlation of the recorded signals to estimate the delay between them. I am sure you are already aware of that, but if you consider the fact that cross-correlation cannot always differentiate between the different (temporal) components of the signal stream (i.e. different sources), it means that the method works best for single sources.

The method Jonas Schwarz proposes is one possible solution to that, but it may or may not work all the time. I created two different signals, each comprised of three (different) harmonics and delayed one of them in respect to the other by 450 samples. Then I summed them and duplicated the sum with an addition of a delay to the duplicate of 50 samples. Finally, I added a little bit of noise (doesn't really change much but still a nice thing to try). The sampling rate I used is 1 KHz, so it is easy to convert from samples to seconds. You can see the plots below

As you can see the main pick is very well localised and it is right on spot. Now, the next higher is annotated. The value given is about 560, which is definitely off. This, of course, is an oversimplified example here, but as you can see, you can't always be sure that the second peak will necessarily correspond to the second source.

Another approach that has been used in other techniques (such as the broadband MUSIC algorithm) is to divide the spectrum into bands (such as with the use Fourier analysis, or filterbanks) and perform the GCC algorithm in each band. This assumes that in each band only one source will dominate (in a sense, this is the same as considering incoherent sources). From there, you could either use histograms (see 3D Localization of Multiple Audio Sources Utilizing 2D DOA Histograms by Delikaris-Manias et al., DOA Estimation with Histogram Analysis of Spatially Constrained Active Intensity Vector by Delikaris-Manias et al. and references therein) or any other way to estimate the direction of the sources (or even the number of them).

However, you have to keep in mind that, as mentioned above, this will not work under all circumstances. Various approaches as to how to use the acquired cross-correlation functions can be found in Considering the Second Peak
in the GCC Function for Multi-Source TDOA estimation with Microphone Array by D. Bechler and K. Kroschel and in Time-Delay Estimation via Linear Interpolation and Cross Correlation by J. Benesty et al.

Finally, a well structured (in my opinion of course) paper with quite some useful information is Real-Time Multiple Sound Source Localization and Counting using a Circular Microphone Array by Pavlidi et al. The authors in this paper have proposed various methods on how one could approach the problem of localisation (angle estimation only, not distance) and various improvements that have found to work quite well. Keep in mind though that this refers to a circular microphone array.


E of a signal using Rayleigh I have to find The energy of a signal using Rayleigh th. the signal is $$ x(t) = A e^{-At } u(t) $$ assuming A>0

Using the classic definition of E , I found that it should be $$ \frac{A}{2} $$

Using Rayleigh I should do
$$  \int_{-inf}^{inf} | \frac{A}{A+i2\pi f}|^2 df $$
This because I previously found X(f)
From this I obtained
$$\frac{A^2}{8 \pi^2  f}  log (\sqrt (A^2 +4\pi^2f^2)) $$ with log from -infinite to infinite.

$$ log \sqrt x = log \frac{x}{2} $$ but now I have no idea how to continue , and if this is correct.. thank you

If I understand correctly, you want to verify the energy calculation in the frequency domain by computing the energy as

$$E_x=\int_{-\infty}^{\infty}|X(f)|^2df\tag{1}$$

with

$$X(f)=\mathcal{F}\big\{x(t)\big\}=\frac{A}{A+i2\pi f}\tag{2}$$

From $(2)$ we get

$$|X(f)|^2=\frac{A^2}{A^2+(2\pi f)^2}=\frac{1}{1+\left(\frac{2\pi f}{A}\right)^2}\tag{3}$$

With $(3)$ and with the substitution $x=2\pi f/A$, the integral $(1)$ becomes

$$E_x=\frac{A}{2\pi}\int_{-\infty}^{\infty}\frac{1}{1+x^2}dx=\frac{A}{2\pi}\arctan(x){\huge|}_{-\infty}^{\infty}=\frac{A}{2\pi}\cdot\pi=\frac{A}{2}\tag{4}$$


Positive or negative sign on Fourier transform formula I have seen both the formula of Fourier transform with positive and negative sign on exponential as  $$ X(\omega)=\int_{-\infty}^{\infty} x(t)e^{-j\omega t}dt$$ and $$ X(\omega)=\int_{-\infty}^{\infty} x(t)e^{j\omega t}dt$$ I am confused which one is the correct formula. I also solved for Fourier transform by taking the following example $$x(t)=\begin{cases} 1, \hspace{5mm} \text{for} \hspace{2mm} |t|<1 \\0, \hspace{5mm} \text{for} \hspace{2mm} |t|>1  \end{cases}$$ and got the same result as $$ X(\omega)=\begin{cases} 2\frac{\text{sin}\omega}{\omega}, \hspace{5mm} \text{when} \hspace{2mm} \omega \neq 0 \\2, \hspace{13mm} \text{when} \hspace{2mm} \omega = 0\end{cases}$$ Can anyone explain whether both the formula for Fourier transform are correct or not?

The definition with the negative in the exponent is the accepted definition of the Fourier transform... however, this is an arbitrary choice.  It could just as easily be defined with $e^{jw}$ and the inverse transform with $e^{-jw}$.


3D Position estimation using 2D camera I have a camera (iPhone), I have a 3D control object in the image that I know its properties very well. (My control object). There is also a secondary object in motion. The ultimate goal is to establish the 3D trajectory of the moving object for a given period of time. (Tracking)

I like to ask, could I find out?

Phone's distance to the control object (for sake of discussion, let's assume camera is at certain height and certain distance neither of them is known but camera is perpendicular to the surface which is known)
The secondary object where I can locate the object in each subsequent frame.  My goal is to estimate its 3D trajectory as I have indicated above.

Bonus question, we can make the system such that the phones distance to the control object can be set (although not preferred), would this help me with the second point?

If you object has 6 known points (known 3D coordinates, $X, Y$ and $Z$) you can compute the location of the camera related to the objects coordinate system.

First some basics.

Homogenous coordinate is vector presentation of euclidean coordinate $(X,Y,Z)$ in which we have appended so called scale factor $\omega$ such that the homogenous coordinate is $\textbf{X} = \omega \begin{bmatrix}X & Y & Z & 1\end{bmatrix}^T$. In your own calculations try to keep $\omega=1$ as often as possible (meaning that you "normalize" the homogenous coordinate by dividing it with its last element: $\textbf{X} \leftarrow \frac{\textbf{X}}{\omega}$). We can also use homogenous presentation for 2D points such that $\textbf{x} = \omega\begin{bmatrix}X & Y & 1\end{bmatrix}$ (remeber that these $\omega, X,Y$ and $Z$ are different for each point, be it 2D or 3D point). Homogenous coordinate presentation makes the math easier.

Camera matrix is $3\times4$ projection matrix from the 3D world to the image sensor:

$$
\textbf{x} = P\textbf{X}
$$

Where $\textbf{x}$ is the point on image sensor (with pixels units) and $\textbf{X}$ is the projected 3D point (lets say that it has millimeters as its units).

We remember that cross product between two 3-vectors can be defined as matrix-vector-multiplication such that:

$$
\textbf{v} \times \textbf{u} = \\
 ( \textbf{v} )_x \textbf{u} = \\
\begin{bmatrix}
        0 & -v_3& v_2 \\
        v_3 & 0 & -v_1 \\
        -v_2 & v_1 & 0
\end{bmatrix}
\textbf{u}
$$

It is also useful to note that cross production $\textbf{v} \times \textbf{v} = \textbf{0}$.

Now lets try to solve the projection matrix $P$ from the previous equations. Lets multiply the projection equation from the left side with $\textbf{x}$s cross product matrix:

$$
(\textbf{x})_x\textbf{x} = (\textbf{x})_xP\textbf{X} = \textbf{0}
$$

Aha! The result must be zero vector. If we now open the equation we get:

$$
    \begin{bmatrix}
        0 & -w& y \\
        w & 0 & -x \\
        -y & x & 0
    \end{bmatrix}
    \begin{bmatrix}
            P_{1,1} & P_{1,2} & P_{1,3} & P_{1,4} \\
            P_{2,1} & P_{2,2} & P_{2,3} & P_{2,4} \\
            P_{3,1} & P_{3,2} & P_{3,3} & P_{3,4}
    \end{bmatrix}
    \textbf{X}
\\
=

\begin{bmatrix}
 P_{3,4} W y - P_{2,1} X w - P_{2,2} Y w - P_{2,4} W w + P_{3,1} X y - P_{2,3} Z w + P_{3,2} Y y + P_{3,3} Z y \\
 P_{1,4} W w + P_{1,1} X w - P_{3,4} W x + P_{1,2} Y w - P_{3,1} X x + P_{1,3} Z w - P_{3,2} Y x - P_{3,3} Z x \\
 P_{2,4} W x + P_{2,1} X x - P_{1,4} W y - P_{1,1} X y + P_{2,2} Y x - P_{1,2} Y y + P_{2,3} Z x - P_{1,3} Z y
\end{bmatrix} = \textbf{0}
$$

With little bit of refactoring we can get the projection matrix $P$ outside of the matrix:

$$
\tiny
\begin{bmatrix} 0 & 0 & 0 & 0 & - X\, w & - Y\, w & - Z\, w & - W\, w & X\, y & Y\, y & Z\, y & W\, y\\ X\, w & Y\, w & Z\, w & W\, w & 0 & 0 & 0 & 0 & - X\, x & - Y\, x & - Z\, x & - W\, x\\ - X\, y & - Y\, y & - Z\, y & - W\, y & X\, x & Y\, x & Z\, x & W\, x & 0 & 0 & 0 & 0 \end{bmatrix}
\begin{bmatrix}
    \textbf{P}_1 \\
    \textbf{P}_2 \\
    \textbf{P}_3 \\
\end{bmatrix}
=
\textbf{0}
$$

Where $\textbf{P}_n$ is the transpose of $n$:th row of the camera matrix $P$. Last row of the previous (big) matrix equation is linear combination of the first two rows, so it does not bring any additional information and it can be left out.

Small pause so we can gather our toughs. Note that the previous matrix equation has to be formed for each known 3D->2D correspondence (there must be at least 6 of them).

Now, for each point correspondence, calculate the first two rows of the matrix above, stack the $2\times12$ matrices on top of each other and you get new matrix $A$ for which

$$
A\begin{bmatrix}
    \textbf{P}_1 \\
    \textbf{P}_2 \\
    \textbf{P}_3 \\
\end{bmatrix}
=
\textbf{0}
$$

As we have 12 unknowns and (at least) 12 equations this can be solved. Only problem is that we don't want to have the trivial answer where
$$
\begin{bmatrix}
    \textbf{P}_1 \\
    \textbf{P}_2 \\
    \textbf{P}_3 \\
\end{bmatrix} = \textbf{0}
$$

Fortunately we can use singular value decomposition (SVD) to force

$$
\|
\begin{bmatrix}
    \textbf{P}_1 \\
    \textbf{P}_2 \\
    \textbf{P}_3 \\
\end{bmatrix}
\|=1
$$

So to solve the the equations, calculate SVD of matrix $A$ and pick the singular vector corresponding to the smallest eigen value. This vector is the null vector of matrix A and also the solution for the camera matrix $P$. Just unstack the $\begin{bmatrix} \textbf{P}_1 & \textbf{P}_2 & \textbf{P}_3 \end{bmatrix}^T$ and form $P$.

Now you wanted to know the distance to the object. $P$ is defined as:

$$
P = K\begin{bmatrix}R & -R\textbf{C}\end{bmatrix}
$$

where $\textbf{C}$ is the camera location relative to the objects origin. It can be solved from the $P$ by calculating $P$s null vector.

(Hartley, Zisserman - Multiple View Geometry in Computer Vision)

Finally, when you calculate the cameras location for two frames, you can calculate the unknown objects locations (or locations of some of the points of the object) by solving two equations for $X$:

$$
\textbf{x}_1 = P_1 \textbf{X} \\
\textbf{x}_2 = P_2 \textbf{X} \\
$$

Which goes pretty much the same way as how we solved the camera matrices:
$$
(\textbf{x}_1)_xP_1\textbf{X} = \textbf{0} \\
(\textbf{x}_2)_xP_2\textbf{X} = \textbf{0} \\
$$

And so on.


Low Level Audio Engine Audio Data Reformatting - Sample Rate and Bit Depth Conversion First off, there seems to be a sparse number of low level audio programming tutorials out there so if anyone knows of some good resources I would be very appreciative. Now on to the topic at hand!
I am using WASAPI for the audio engine that I am trying to develop and I am having difficulties reformatting read audio file data into a format that the audio device is prepared to receive. Specifically I am wondering how to on the fly convert the sampling rate and bit depth of a read audio file into the requested format of the device. I saw an example where linear interpolation was used but read elsewhere that such a filter introduces too many artifacts such as aliasing, so again I'm wondering what practical algorithms would be for reformatting audio data in a real time application.

Sinc interpolation introduces no aliasing artifacts to data that was originally properly bandlimited before sampling.  (But needs to be slightly modified, as a Sinc or Dirichlet kernel has an infinite duration.)
One canonical reference on why reconstruction works and how to interpolate is here: https://ccrma.stanford.edu/~jos/resample/
I have a simple implementation in Basic here: http://www.nicholson.com/rhn/dsp.html#3
The other common method is low-pass filtered integer ratio upsampling followed by low-pass filtered integer ratio downsampling by the appropriate two ratios.  Another common method is to generate a polyphase FIR filter bank (of multiple phases of windowed Sinc (or otherwise optimized FIR filters) interpolation kernels.  Another common method is to use polynomial interpolaters of each of the Sinc kernel's humps, resulting in a Farrow filter bank.
You need to carry enough bits of arithmetic intermediate values in order to keep your quantization noise floor below whatever your spec might be.


Help with filter stability, ROC and output of a system I am trying to solve these two questions, but I am stuck with the filter stability and ROC. Hope somebody could help me with it.

According to what I understand, the third option is correct. Fourth option is also incorrect as it only contains one pole and on zero. I think the first and second options are incorrect; the system is unstable and ROC doesn't include unit circle because the pole lies at z = -2 and the system is casual. Not sure though. I'd appreciate some help in this regard or if any of the other options are incorrect.

For this question, the pole is outside the unit circle, so I thought the system must be unstable, but then I realize the output depends on future input (x[n+1]), which makes it non-casual, also number of zeros>poles. Since for a non-casual, the ROC moves inward, this system includes unit circle and is stable? But doubtful about the second part of the option regarding being dependent on input.
I wrote the following code to find y[n],
n = 100
w = np.linspace(-np.pi,np.pi,n)
z = np.exp(1j * w)
H = z* (z - 0.5)/(z + 2)
t = np.arange(-n/2, n/2)
x  = sp.unit_impulse(n) + 0.5 * sp.unit_impulse(n,1)
plt.plot(t,ifftshift(ifft( H * fft(x))))

Got the following result,

I don't think, 1st, 2nd and 5th options are correct, looking at figure. As I increased the value of n, the peak became narrower, so for the 4th option, I don't know if it is showing diverging oscillation. The only peak in this question is at zero.
Any help is appreciated.
Edit: Question 1 seems to work with the above logic. Stuck with no. 5 now.

For any filter with AR(1) component the coefficient must be smaller than in order to be stable. Otherwise any input is infinitely multiplied by a non decreasing factor.


Matlab: What should be the BER performance for BPSK using Constant Modulus Algorithm equalizer I need help in plotting the Bit error curve or the symbol error curve for BPSK modulation scheme for varying Signal to Noise ratios or Eb/N0. The plot should show the simulated versus the theoretical curve, but I cannot figure out how to mitigate the problems when using the Constant Modulus Algorithm as an Equalizer which are

(1) I want to use the filter function instead of conv in order to model a moving average channel model, chanOut = filter(ht,1,s). But, when I use filter, I am getting an error. How can I use filter function here?

(2) Bit error rate calculation

I am still unable to use filter and unsure if BER curve is proper or not.
Below is the code I wrote:

UPDATE after the answer posted, I have included all the suggestions. The error has been mitigated but now the curve for CMA is coming out to be flat instead of the one shown in the answer. Would like to know what I should do to get similar curve.

Well, I took a look at your code, and spent a lot of time on it, and I have discovered some mistakes, some practicals and some theoretics. Here are my answers:

(1) You can use the function filter. Just remember that it returns a vector of the same length of the input, while conv returns the $ N_x + N_h - 1 $, where $ N_x $ is the length of your input and $ N_h $ is the length of your channel. This difference needs to be take in consideration when you create your noise vector. You won't be able to add the noise vector to your signal if you don't pay attention to this. Read the documentation to get more information.

(2) As I said, there is a lot of mistakes at your code:

First I changed the vector noise to only add a real number to the signal, and not a complex number.
Second, I changed the equation

w = w - mu * e(i) * yy * x;

to

w = w - mu * x * e(i) * yy.'

that I took from here.

There is some mistakes in your CMA algorithm. Here are some lines and the coments:

yy = w(:,i).'*x; Remeber to conjugate w in this. In BPSK we only have real numbers so this is not a problem, but if you had complex numbers, you must conjugate w.

e(i) = abs(yy)^2 - 1; In your original code, you didn't use abs(y) to calculate the error. Again, it is not a problem if you are only working with real numbers, but it will become a problem with complex numbers.

The way that you are getting the vector x inside the CMA loop is not getting all the samples since the beginning. So I iniciated the vector x before the loop as:

x = zeros(Le, 1);

and inside the loop I used:

x = [r(i); x(1:end-1)];

It will start CMA from the beginning of the receiving symbols.

This line doesn't make sense:

sb=w'*x;

This only equalizes ONE sample. The correct way to perform the equalization once you have the filter coefficients is:

sb = filter(w(:,N+1), 1, y);

Remember what I said about the function filter and the function conv. Read the documentation.

The next problem is in your number of errors computing. This is a very trick, but you need to remember (or learn) that every equalizer generates a delay in your signal. So you cant do ip- ipHat directly. You need to compensate this delay. I used:

nErr_CMA(ii) = size(find([ip(1:end-Le+1) - ipHat(Le:end)]),2);

If you don't know about this delay, you need to study a litle bit more.

In the plot part, I used the following:

figure
semilogy(Eb_N0_dB,theoryBer,'b'),grid;
hold on;
semilogy(Eb_N0_dB,simBer_CMA,'r-','Linewidth',2);
%axis([0 14 10^-5 0.5])
grid on
legend('Theoretical BER','sim-CMA');
xlabel('Eb/No, dB');
ylabel('Bit Error Rate');
title('Bit error probability curve for BPSK in ISI with CMA equalizer');

I also increased N to N = 1e5, so as to perform a Monte Carlo simulation properly.
I changed Le to Le = 5;
I changed the w initializations for w = ones(Le, 1); to avoid that the update part of the CMA equation results in zero.
Also, change your channel. I didn't have time to verify that, but I strongly suspect that the channel [0.2 0.9 0.3] has some properties that makes it difficult for the CMA to equalize. I used ht = [1 0.5 0.2]. Try using this channel.

Ok, I gave you a lot of answers. I don't know if there is any more mistakes, but here I was able to get a BER curve that looks right. If it is right or not, it is up to you. I suggest you to look for articles about communications systems using CMA algorithm and compare your results.

Let me know if I need to let anything I said more clear. Good luck on your research.


Maximising each element in a matlab array Is it possible to construct an array A, without looping, where A(i,j) = max(B(i),C(j)) and B and C are column vectors?

If the matrices are not too big, repmat could work:

B = [2 3];
C = [0 1 2 3];
% Vectorize the vectors for a simplication
B = B(:);
C = C(:);

A = max(repmat(B,1,length(C)),repmat(C',length(B),1));

For those interested, Comparing BSXFUN and REPMAT and Matlab - bsxfun no longer faster than repmat? address its relative efficiency with respect to repmat. Apparently,  bsxfun used to be much faster, sometime ago. The difference seems shallower with the recent Matlab versions.


Set of all images with a certain MSE between themselves and a given image Mean squared error between two $M\times N$ images, $I[m, n]$ and $D[m, n]$ is defined as:
$$
MSE(I, D) = \frac{1}{M\times N}\displaystyle \sum_{m=1}^{M}\sum_{n=1}^{N}(I[m, n]-D[m, n])^2
$$
For a given image, $I[m, n]$, and a particular value, $\epsilon $, how to calculate all images, $D[m, n]$, that:
$$
MSE(I, D) = \epsilon
$$

Well, you know, there are infinitely many. Start with $D=I$, change one particular pixel to $I[m,n]+\sqrt{\epsilon}$. You have your first candidate that satisfies ${\rm MSE}(I,D)=\epsilon$. Now you can generate a new image $D$ by adding a perturbation to an arbitrary subset of pixels (say, add $\Delta$ to a single pixel) and correct the resulting difference in another pixel (say, subtract $\Delta$ from another pixel). As this works for any subset and any perturbation $\Delta$, the number of images you can generate is infinite.

Of course, if your images are INT8, this makes the set finite, yet still impractibly large. Only the amount of perturbations you can apply to single pairs of pixels is in the order of $2^{7MN}$, considering larger groups of pixels makes this even more larger.

Another way of looking at it: an $M \times N$ INT8 image has $256^{MN}$ degrees of freedom. ${\rm MSE}(I,D)=\epsilon$ is a single contraint. It's a quadratic one yes, but still it does not significantly reduce the degrees of freedom you have.


Oversampled polyphase filter banks i would a confirm of this: with polyphase structure is possible only design a filter bank with a INTEGER oversampling ratio? For non-integer i've seen the weighted overlap-add metod, is right?

thanks
S.

If you have a 1D uniform $M$-channel multirate filter bank (FB), and $N$ is the decimation ($M \ge N$),  rational oversampling ratios  ($M/N$) are indeed possible (see Fig. 2), or from Figure 1 in Optimization of Synthesis Oversampled Complex Filter Banks, J. Gautier et al., IEEE Trans. Signal Processing, 2009 (doi), which provides optimized synthesis FB with any ratio (like $7/4$).

More precisely: for an analysis filter bank, you can choose any decimation you want. But you may loose information, and the invertibility of the analysis/synthesis system. In a  lot of practical applications (for instance in speech processing), integer oversampling ratios are quite common, by choosing $N = M/k$, for instance $k=M/2$ or $k=M/4$. Such choices make the  computations and some calculations much simpler: for instance, an inverse synthesis filter bank may take a closed form, related to the  analysis window for DFT FBs.

However with redundancy, there are an infinity of inverses, and polyphase formulations help finding synthesis FBs.

A set of useful books (to me):

Subband compression of images: principles and examples, T. Ramstad et al., 1995
Wavelets and Subband Coding, M. Vetterli et al., 1995
Wavelets and filte banks, G. Strang et al., 1996


Echo cancelling using autocorrelation function I was given a problem, but I couldn't solve. I did some researches but I still didn't figure it out. Here is the problem:

An audio signal $ s(t) $ is generated by a speaker reflects in a wall with reflection coefficient $ \mu (\mu < 1) $. The signal $ x(t) $ is recorded by a microphone close to the speaker and far from the wall, after sampled, is given by:

$ x(t) = s(t) + \mu s(t - k) $

Where $ k $ is the delay given in samples due to the echo.

So, my goal here is to estimate $ \mu $ and $ k $, observing the autocorrelation function of $ x(t) $, $ r_{xx}(l) $, that can be writen in function of the autocorrelation function of $ s(t) $, $ r_{ss}(l) $. I was able to compute this autocorrelation function and I got:

$ r_{xx}(l) = (1 + \mu^{2})r_{ss}(l) + \mu r_{ss}(l-k) + \mu r_{ss}(l+k) $

So, it obvious has three peaks: when $ l = 0 $, when $ l = k $ and when $ l = -k $.

Assuming that $ r_{ss}(l) = 0 $ for $ |l| > k/10 $, I think that I can obtain the value of k observing the side peaks, and the value of $ \mu $ observing the central peak, but i wasn't able to formalize it mathematically. Could anyone help me to do so?

Thanks in advance.

Haven't you already got it there?

$$
\hat{k} = \arg \max_{l > \frac{k}{10}} r_{xx}(l)
$$
and
$$
\hat{\mu} = \sqrt{\frac{r_{xx}(0) }{ r_{ss}(0) } - 1}
$$

The R code below outputs the figure and:

  "k Estimate: 613  vs  613  mu Estimate :  0.747619585689531  vs  0.768768"

R Code

# 26617

T <- 10000
mu <- 0.768768
k <- 613

s <- filter(runif(T,-1,1), rep(10/k,k/10), circular = TRUE)

x <- s + mu*c(rep(0,k),s[1:(T-k)])

par(mfrow=c(2,1))
r_xx = acf(x, lag.max=1000, type="covariance")

min_idx <- floor(k/10)
mx <- which.max(r_xx$acf[(min_idx+1):10000])

khat <- min_idx+mx-1 # R indices start from 1, not 0

points(min_idx+mx,r_xx$acf[khat], col="red", pch=19)

r_ss = acf(s, lag.max=1000, type="covariance")

muhat <- sqrt( r_xx$acf[1] / r_ss$acf[1] - 1 )
print(paste("k Estimate:" , khat, " vs " , k, " mu Estimate : ", muhat, " vs ", mu ))


Using Kalman Filter to fuse two sensor readings of the same type I have always used Kalman Filter to smooth a signal comes from one sensor only.

I would like to know if Kalman Filter could be used to fuse data coming from two different sensors that provide the same type of reading. An example is fusing the position data that comes from the GPS with the position data that is calculated from the IMU. In the light of this example, I wish to have the following:

State vector: [x;y]
measurements vector: [x_imu;y_imu;x_gps;y_gps]

Let me know if you need more information/clarification.

The answer to the first question is yes, it is possible to fuse data coming from the same type of sensor. For example, the fusion of two different accelerometers, each with a measurement quality, gives a weighted average where the weights are the inverse of individual variances (see Kays's[1] book for details).

The answer for the second is yes as well, it is possible to fuse IMU with GPS measurements. But note the one returns a position itself, while the other measures the second derivative of the position. The fusion demands well elaborated models. There are plenty of papers and books with models and its theory. I would recommend, for example, Farrel's[2] and Titterton's[3] books.You can easily find implementations of the fusion filters in internet. Open source projects of drones, for example, have it.

[1] KAY, Steven; Fundamentals of Statistical Signal Processing, Estimation Theory.

[2] FARREL, Jay A.; The Global Positioning System and Inertial Navigation.

[3] TITTERTON, D. H.; Strapdown Inertial Navigation Technology.


Why dont sin and phase shifted cosine overlap? I plotted a sine wave(blue) and a 90 degrees phase shifted cosine wave (red) expecting them to overlap each other.

x = 0:0.001:1;

%ploting sine
y = sin(2*pi*x);
plot(x, y);

hold('on');

%ploting cosine
z = cos(2*pi*x-90);
plot(x, z);

The two graphs do not overlap.
Whats wrong with my assumption ?

You're confusing degrees and radians.

$$\sin(x)=\cos(x-\pi/2)$$

$\pi/2$ radians corresponds to 90 degrees.


Signal squared in time domain In some publication I have found a formula that for some reason squares signal in time domain.
I am analyzing it now.

What should be the Fourier transform of a signal that was squared in time domain? I do not understand why $\mathcal{F}[sin(x)^2]$ results in peak in zero and why side peaks are moved.

I fully understand the fact that multiplication in time domain is convolution in frequency domain and squaring is just $x\times x$.

If somebody squares a real signal in the time domain, he probably deals with the signal's energy. Since this is a pointwise product, this results in a convolution in the Fourier domain:

$$ \mathcal{F}(x\times x) = \mathcal{F}(x)\ast  \mathcal{F}(x) $$

In some cases, as stated by @MBaz, sometimes trigonometric identities can explain the frequency shift, as:

$\cos^2 x=\frac{1}{2}(1+\cos2x)$
$\sin^2 x=\frac{1}{2}(1-\cos2x)$

From the above formulae, one can see that:

they are positive, and the non-oscillatory term is $1/2 \cos(0x)$, explaining a peak at frequency zero,
the scalar before $x$ has doubled, explaining the peak moves.


Benefits of oversampling in a realtime noise-cancellation system When sampling a signal, it is standard practice to use a sampling rate that is greater than two times the bandwidth of the signal in order to avoid aliasing. The result is then low-pass filtered (and sometimes upsampled prior to this) in order to get rid of the harmonics.

Unfortunately, in control systems applications preventing aliasing in not enough. The system also needs to respond quickly enough to the input (as determined by a spec), and a low sampling rate will limit this performance. Specifically, if the system samples a signal at time t, then any event that occurs immediately afterwards will be ignored until one sampling period later.  Thus, there are clear performance advantages to having a sampling rate much greater than two times the signal bandwidth. I'm looking to quantify this performance improvement.

Consider an open-loop control system that attempts to cancel out an incoming signal by outputting an inverted version of this same signal. In other words, active compensation through destructive interference. If the system has infinite magnitude resolution (both the ADC and DAC have infinite bits) and an infinite sampling time, then the inverted output signal will be a perfect inverse of the original and the signals will cancel out completely. But if the system has a finite sampling rate, then the output signal will not be a perfect inverse of the original and the two will not cancel completely. The lower the sampling rate, the worse the cancellation will be (even if no aliasing occurs).

Q: Assuming the following:

The input signal is a  sinusoid of frequency $f_0$ and magnitude $A$
The system has an ADC with a sampling rate of $f_s$ where $f_s > 2 f_0$
The system uses sample-and-hold (ZOH) interpolation to output an inverted sinusoid and add it to the original signal
The system has zero time delay between an ADC sample and the corresponding DAC output (ADC settling time, DSP calculations, etc. are instantaneous)

How do I calculate the frequency and magnitude of the resulting frequency components?

The system has zero time delay (sampling, calculations, etc. are instantaneous)

This is an invalid assumption. Almost every digital system has at least one sample delay. DAC and ADC or clocked off the same master clock. So at sample time n, the ADC is acquiring sample n while at the same time the DAC is outputting sample n-1. The only way around this is using a much higher master clock, donwsampling that clock and clocking the A/D and D/A on different phases. However that's equivalent to running at a higher sample rate in the first place.

So your ability to cancel is given by how much the amplitude of the sine wave changes during 1 sample. The residual signal is simply given by $$e(t) = \sin(f_0 \cdot t) - \sin(f_0 \cdot (t-T_s)) = 2 \cdot \cos(f_0 \cdot (t- T_s/2)) \cdot \sin(f_0 \cdot T_s/2)$$
 where $T_s = 1/f_s$ is the sample time. The error is simply a function of the ratio of signal frequency and sample rate. For $f_0 << f_s$ the error is very small.

The only way to get around this is to crank up the sample rate or use some a priori knowledge of the signal to do some prediction and speculative cancellation.