fft - spectral leakage from noise into target frequency? I have an acoustic setup. A sender, emitting a clean sine signal. The receiver is far away. In the time domain, the sine is hardly visible because the noise is nearly as big as the signal. In the frequency domain, the noise is distributed across all frequencies.

I need to extract the amplitude at the target frequency, I do not care about the phase.
Currently, I calculate the length of my fft-window to contain exactly N periods of my target frequency (n_bins=N * samplerate / frequency) to avoid a loss in magnitude due to spectral leakage.
To remove the noise at this frequency, I do the same with a measurement w/o signal to estimate the noise background and subtract it.

My question now is, should I still use a window-function like blackman? From my understanding, I could now have leakage from noise in adjacent bins into my signal?

I have tried using a window-function in a python-example for this, with simulated data. Here, when using a window, the amplitude at my target frequency was actually higher than expected, while the results without the window were closer to reality.

EDIT:
Python source code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import blackman as w_function

target_freq = 12000
samplerate = 1000000.
window = True

# generate data
x = np.array(range(0,17166)) / samplerate
y = 290*np.sin(2*np.pi*target_freq*x) # signal
y += 150*np.sin(2*np.pi*(target_freq - 451)*x) # noise near the signal
y += 232*np.sin(2*np.pi*200185*x) # noise far away

# find closest multiple
reduced_len = len(y)
while not float(reduced_len * target_freq / samplerate).is_integer():
    reduced_len -= 1

def calc_fft(window):
    # prepare window
    if window:
        wndw = w_function(len(y))
    else:
        wndw = np.array([1] * len(y))
    wndw_correction = float(len(y)) / float(np.sum(wndw))  # amplitude     correction factor for window

    # calculate fft, apply corrections for correct amplitude     representation
    rft = np.abs(np.fft.fft(y * wndw, reduced_len, norm=None))
    rft *= 2 * wndw_correction / float(len(rft))

    # calculate corresponding frequencies, show only positive frequencies
    fft_freq = np.fft.fftfreq(len(rft), 1 / samplerate)
    rft = rft[np.argsort(fft_freq)]
    fft_freq = fft_freq[np.argsort(fft_freq)]
    low_ind = (np.abs(fft_freq-0)).argmin()
    fft_freq = fft_freq[low_ind:]
    rft = rft[low_ind:]
    return fft_freq, rft

fft_freq, rft_window = calc_fft(True)
fft_freq, rft_no_window = calc_fft(False)

plt.plot(fft_freq, rft_window, label="window")
plt.plot(fft_freq, rft_no_window, "-.", label="no window")
plt.legend()
plt.show()

Without any noise (just comment lines in python out), the fft without window gives the exact amplitude (290), with window the amplitude is to high:

With noise, the amplitude with window is still to high, and the amplitude without window is somewhat to low:

I am assuming you are wanting to measure the phase of the signal, perhaps the amplitude, since you already know the frequency.

I don't think a window will help you at all for the reasons you have cited.  In addition, since you are only interested in a single bin value, that is all you have to calculate instead of doing a whole FFT.

A much better noise reduction approach is to smooth your signal in the time domain.  There are many ways to do this.  My favorite way is to run a simple exponential IIR in the forward direction, then in the backward direction, and average the results of the two passes.  This process leaves the phase unaffected and attenuates the amplitude by a known frequency dependent factor.  This takes out a lot of the noise since a DFT inherently does a least squares fit and outlier data points have a disproportionate effect.

Others may recommend other filters, such as symmetric FIRs, but they take a little bit more calculation.

You may also want to turn up the volume on your emitter.

Hope this helps.

Ced

=========================================

Followup:

I've looked your code over and there are some issues I see:

First, your calculation for reduced_len is very fragile/ineffecient and depends on good values for samplerate and target_freq.  You are better off using this code:

samples_per_cycle = samplerate / target_freq
cycles_per_frame  = 17166 / samples_per_cycle
whole_cycles      = int( cycles_per_frame )
new_length        = whole_cycles * samples_per_cycle
reduced_length    = int( new_length + .5 )

Second, it appears that you are constructing your window function for the full length of y, not reduced_len.

Third, your window correction calculation makes no sense to me.

Now, turning to your problem in general.  There are basically two ways to get rid of the influence of nearby tones in a DFT.  Windowing or estimating its parameters and removing its contribution from the DFT.  Obviously, the latter approach is more complicated.  I have three blog articles that give the math on how to do this and if you are interested I will elaborate further.

For windowing, the one I use for display purposes, and the simplest to understand is the VonHann window.  The window correction factor is a factor of 2 as $ Y[k] = -.25 X[k-1] + .5 X[k] - .25 X[k+1] $ is what happens in the DFT.  For calculation purposes, I don't do windows, they muddle values.

For your application, it doesn't seem to me that important that you calculate the actual amplitude.  If you are measuring changes in amplitude on a relative basis any proportional value will work.

If you want to contact me directly, you can find my email address in my bio in my blog articles.  You can find my blog link on my profile page.

Ced


Inverse DFT of $-i\pi \text{ sgn}(f)$ the following MATLAB code should plot the iFFT of $-i \pi sgn(f)$ which has to be $\frac{1}{t}$ in continuous time. As you can see, the output kinda looks like a 1/x function, but is modulated.
How can this phenomena be explained and is there a way to get the desired $\frac{1}{t}$ function?

clc;
clear;
close all;

L=1000;
fStart=-100;
fEnd=100;
df=(fEnd-fStart)/L;
f=linspace(fStart,fEnd-df,L);

X=-1i*pi*sign(f);
x=ifft(fftshift(X));

figure()
ax1=subplot(3,1,1);
ax2=subplot(3,1,2);
ax3=subplot(3,1,3);

stem(ax1,imag(X));
plot(ax2,fftshift(real(x)));
plot(ax3,fftshift(imag(x)));

title(ax1,'spectrum');
title(ax2,'Re(iFFT)');
title(ax3,'Im(iFFT)');

Edit: for clarification: I especially wonder, why every other sample in the result is zero:

There are two issues here. First of all, you go from continuous time to discrete time, and second, you go from a continuous spectrum to a sampled spectrum. Note that the impulse response of an ideal discrete-time Hilbert transformer with frequency response

$$H(e^{j\omega})=-j\text{ sgn}(\omega),\qquad \omega\in[-\pi,\pi]\tag{1}$$

is given by

$$h[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}H(e^{j\omega})e^{jn\omega}d\omega=\begin{cases}\frac{2}{n\pi},& n\textrm{ odd}\\0,&n\textrm{ even}\end{cases}\tag{2}$$

So the zero samples at even indices don't come as a surprise.

However, with your method you don't get $(2)$ because you sample the spectrum, which results in aliasing in the time domain. You also get a non-zero imaginary part because your discrete spectrum does not satisfy the condition

$$X[k]=X^*[N-k]\tag{3}$$

where $^*$ denotes complex conjugation. Using the following definition for $X[k]$ will give you a purely real-valued time-domain sequence (up to numerical errors):

$$X[k]=\begin{cases}0,&n\in \{0,N/2\}\\
-j,& n=1,\ldots,N/2-1\\
j,&n=N/2+1,\ldots,N-1\end{cases}\tag{4}$$

where I've assumed that $N$ is even.

The inverse discrete Fourier transform (IDFT) of $(4)$ is given by

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi kn/N}=\begin{cases}\frac{2}{N}\cot\left(\frac{\pi n}{N}\right),&n\text{ odd}\\0,& n\text{ even}\end{cases}\tag{5}$$

Consequently, it is expression $(5)$ which you should expect instead of a sampled version of $1/\pi t$. Of course, $(5)$ approximates $(2)$, but the approximation gets worse for increasing magnitude of the index $n$ ($n\in [-N/2,N/2-1])$:


a-weighting signal without digital filter I essentially want to a-weight a signal by using the fft in Python.
But I guess I'm making some logical faults.

My Steps:

Initialize Sine: y = np.array(np.rint(32768 * np.sin(1000 * 2.0 * np.pi * np.linspace(0,1,44100))))
take RMS: print(20.0 *np.log10(rms_flat(y))) - gives me $87\textrm{ dB}$
fft with hanning window:yf = np.abs(np.fft.fft(samples * np.hanning(len(samples))))
take RMS: gives me around $130\textrm{ dB}$ - why is it not at $87\textrm{ dB}$ anymore?
finally a-weight spectrum - see below - get Inf RMS...

import matplotlib.pyplot as plt
import numpy as np

#samplerate
Fs = 44100

#samplecount
A = Fs

#Sample Interval
Ts = 1/Fs

#Zeitvektor
t = np.linspace(0.0, Fs*Ts, Fs)

#Sinus Freq
ff = 1000

y = np.array(np.rint(32768 * np.sin(ff * 2.0 * np.pi * t)))

fig = plt.figure(figsize=(10,10))
fig.subplots_adjust(hspace=0.7)
fig.subplots_adjust(wspace=0.25)
plt.subplot(321)
plt.title("Samples")
plt.ylabel("Amplitude")
plt.xlabel("Time")
plt.plot(t, y)

def rms_flat(a):  # from matplotlib.mlab
    """
    Return the root mean square of all the elements of *a*, flattened out.
    """
    return np.sqrt(np.mean(np.absolute(a)**2))

print(20.0 *np.log10(rms_flat(y)))

def fft(samples):
    yf = np.abs(np.fft.fft(samples * np.hanning(len(samples))))
    xf = np.linspace(0.0, 1 / (2.0 * Ts), Fs / 2)
    plt.subplot(322)
    plt.title("Abs Pos Spectrum")
    plt.ylabel("Amplitude")
    plt.xlabel("Frequency")
    plt.plot(xf, 2.0/A * yf[:A//2])
    print(20.0 * np.log10(rms_flat(yf)))
    return(yf, xf)

yf, xf = fft(y)

def aWeighting(yf, xf):
    f1 = 20
    f2 = 107
    f3 = 737
    f4 = 12194
    N = len(xf)
    i = 0
    yfa = []
    while i < N:
        a = 20.0 * np.log10((f4 ** 2 * xf[i] ** 4) / ((xf[i] ** 2 + f1 ** 2) * np.sqrt(xf[i] ** 2 + f2 ** 2) * np.sqrt(xf[i] ** 2 + f3 ** 2) * (xf[i] ** 2 + f4 ** 2))) + 2.0
        yfa.append(20.0 * np.log10(yf[i]) + a)
        i += 1

    yfa = np.array(yfa)
    print(rms_flat(yfa))
    plt.subplot(323)
    plt.title("A-weighted Spektrum in dB")
    plt.ylabel("dB")
    plt.xlabel("Frequency")
    plt.plot(xf, 2.0/A * yfa[:A//2])
    plt.show()
    return(yfa, weight)

yfa, weight = aWeighting(yf, xf)

The FFT has a processing gain of 10Log(N) that would apply to your measurements of the standard deviation, where N is the length of the signal, and the Hamming Window has a coherent gain of 0.54 (which is a loss). Refer to fred harris' famous paper On the Use of Windowing that details these gains and losses.

Let's see if this works out. You have 44100 samples, so the net processing gain is $10Log10(44100\times.54) = 43.8 \text{ dB}$

Your peak sinewave is 32768 so the rms should be $20log_{10}(32768/\sqrt{2})=87.3 \text { dB}$ which is consistent with your measurement.

After taking the window and the FFT, you saw the RMS signal increase to 130 dB, while you measured 87 dB before the FFT. 130-87 = 43 dB.

The difference of 43 dB is consistent with that predicted by the formula for the net processing gain given above.

The following is more details that may further help to understand how this occurs:

The coherent gain of your input signal is N. Your input signal in time is

$$A sin(\omega t) = \frac{A}{2}\frac{e^{j\omega t}}{j}-\frac{A}{2}\frac{e^{j\omega t}}{j}$$

The important item to see from Euler's expansion above is the sine is composed of two exponentials each with an amplitude of half of the sine wave itself. After the FFT, each of these is increased by a factor of N due to the coherent processing gain in the FFT.

Thus each of the two tones in the FFT would have a magnitude of $\frac{AN}{2}$ (if you sampled coherently at an exact integer number of samples this would be exact, otherwise we will see it slightly reduced from this due to "scalloping loss").

If sampled coherently, every other bin besides the two tones described will be zero (we could further show through Parseval's theorem the general case where all the energy in the signal is conserved, but to keep my explanation simple I will describe just the case of coherent sampling, meaning the sampling clock is an integer multiple of the sine wave).

Therefore when computing the rms, since there are only 2 non-zero values but we are using N samples, the rms calculation of the FFT would reduce to:

$$\sigma = \sqrt{\frac{\left(\frac{AN}{2}\right)^2+\left(\frac{AN}{2}\right)^2}{N}}$$

Which reduces to

$$\sigma = \sqrt{N}\frac{A}{\sqrt{2}}$$

Which we recognize as the $\sqrt{N}$ times the rms of our input sinusoid.  When we window the signal before taking the FFT, the coherent amplitude (that previously increased by N) would now be reduced by the coherent gain of the window, which for the case of the Hamming window was 0.54.

Thus the standard deviation of the FFT is expected to be $\sqrt{(N\times0.54)}$ times larger than the standard deviation of the input sinusoid for the case of a Hamming Window being used. So in dB:

$$20log_{10}(\sigma) = 10log_{10}(N\times0.54) $$


Fourier transform of a Fourier transform I generate a Gaussian noise and then I filter it with a passband FIR Kaiser window filter. When I perform the Fourier transform of the output of the filter and plot its magnitude spectrum, it is filtered in the correct frequency band. Now, I perform a second Fourier transform on the output of the first Fourier transform (I'm performing it on a complex vector this time and in the frequency domain), and when I plot the magnitude spectrum I see something very similar to the original filtered noise. I know this isn't time domain, I know this spectrum is not even because the input signal is complex. But as the time plot of the filtered noise and the absolute value plot of the Fourier transform (of the Fourier transform of the noise) are visually similar a doubt arised on me.

Why is the amplitude spectrum of the complex vector obtained after performing a Fourier tranform on the filtered noise once again occupying the whole vector (the whole frequency band) as if they were time samples spread in the time vector?
Is the negative symbol in the inverse Fourier transform of the Fourier transform of the noise having no effect so that it is visually similar performing a Fourier transform of Fourier transform of noise and performing an inverse Fourier transform of Fourier transform of noise?

This doubt arised to me when I was trying to filter on a matched filter implemented in Matlab using IFFT(FFT(signal_1).*FFT(signal_1_inverted)) and I accidentally filtered the output of an already Fourier transformed signal and coudn't notice the difference (not on the output of the filter but on the output of the first FFT stage). And this leads me to another question.

Suppose I would like to detect the likelihood of 2 signals by comparing their spectrums and I want to measure the likelihood by using a matched filter. Knowing that the convolution in time is the multiplication in frequency, is the convolution in frequency the multiplication in time and I get the same result if I use a matched filter (sort of convolution operation) on the spectrums than if I perform a sample by sample multiplication in time domain of the signals whose spectrums likelihood I want to determine?

Applying the discrete Fourier transform (DFT) twice to a time-domain signal, should produce the original time-domain signal in a circularly time-reversed, and linearly scaled (by N) form; as the inverse DFT is very similar to the forward DFT, except the sign of the complex exponential,a linear scale by 1/N.
The inverse DFT of $X[k]$ is:
$$ x[n] = \frac {1}{N} \sum_{k=0}^{N-1} {X[k]e^{j \frac{2\pi}{N} k n}}$$
By reverting the sign of the exponent $e^{j \frac{2\pi}{N} k n}$ into negative, by a change in the sign of the $n$ of $x[n]$ what you get is:
$$  N x[(-n)_N] = \sum_{k=0}^{N-1} {X[k]e^{-j \frac{2\pi}{N} k n}} $$
which is N times the circularly reversed (obtained via modulus operator denoted by $(.)_N$) version of the original signal $x[n]$, which is obtained by an apparently forward but effectively inverse DFT applied to $X[k]$ which is treated as a time-domain sequence, which is what you do at your second stage of applying forward FFT to $X[k]$ of $x[n]$.


Conceptualising the continuous time unit impulse function as derivative of unit step This a very newbie question.
I just watched Lecture 3 of Oppenheim's Signals course and he defines  here the continuous time function as the derivative of the unit step function like so:
$$ \delta_\Delta(t)=\frac {du_\Delta(t)}{dt}$$
and that $ \delta(t)= \delta_\Delta(t) $ as $\Delta \to 0$
He claims that the derivative is equal to 1 no matter the value of $\Delta$, because that derivative can be interpreted as the area of a rectangle with sides $\Delta$ and $\frac 1 \Delta$
I can't conceptualise this the way the function $u_\Delta(t)$ is drawn at all. If the function is linear, that is, $y = mx + b$ passes through the origin, meaning $b=0$, and we can see that it has the point $(\Delta, 1)$ we can easily tell that $m= \frac 1 \Delta$ and that should be the derivative.
Can someone explain to me the error in my line of thought? Why is the derivative the area of a square?

I didn't look at how the function $\delta_\Delta(t)$ is defined in that lecture, but the convention is that it's a rectangular pulse of area being $1$. The definition is :
$$\delta_\Delta(t) = \begin{cases}{ \frac{1}{\Delta} ~~~, ~~~ 0 < t < \Delta \\ ~0~ ~~~,~~~\text{otherwise} } \end{cases} $$
As you can see the area of that pulse, $A(\Delta) = \frac{1}{\Delta} \times \Delta = 1$ for any finite value of $\Delta$. According to calculus, that area is still $1$ even when the parameter $\Delta$ takes on inadmissible values of $0$ or $\infty$ through a limiting process. In particular the continuous Dirac delta function is defined as $\delta(t) = \lim_{\Delta \to 0} \delta_{\Delta}(t)$ ; when the limit of $\Delta$ goes to zero.
The rectangular pulse of unit area, $\delta_{\Delta}(t)$,  is only one of many other possible ways of defining $\delta(t)$ through a limit of functions process.


Z-transform of an affine function What is the transfer function of the system described by the following affine input ($x$)-output ($y$) relationship:
$$
y[n] = \alpha x[n] + \beta.
$$
Using the Z-transform we find:
$$
Y[z] = \alpha X[z]+ \beta
$$
and
$$
H[z] = \frac{Y[z]}{X[z]} = \alpha + \frac{\beta}{X[z]}
$$
But I don't know how to take it from here to compute the transfer function! I'm not sure if this can be done. Any help would be appreciated.

That system cannot be described by a transfer function because it is not linear. It neither satisfies the condition of additivity nor does it satisfy the condition of homogeneity.
If $y_1[n]$ is the response to an input $x_1[n]$, and $y_2[n]$ is the response to an input $x_2[n]$, then the response to $x_1[n]+x_2[n]$ doesn't equal $y_1[n]+y_2[n]$ (additivity is violated). Furthermore, the response to $\alpha x_1[n]$ doesn't equal $\alpha y_1[n]$ (homogeneity is violated).
This is unrelated to existence questions concerning the $\mathcal{Z}$-transform. The $\mathcal{Z}$-transform of a constant $\beta$ doesn't equal $\beta$, but it simply doesn't exist.


What happens when I try to resample a speech recording from 8kHz to 16kHz? I have a question regarding resampling in the context of speech. Given a speech recording sampled at 16kHz, downsampling to 8kHz will basically remove half of the samples (each 16000 samples become 8000). Now, I'm wondering about the inverse senario, given a telephony quality speech recording (300Hz -> 3400Hz) sampled at 8kHz, what happens if I try to resample the signal to 16kHz ? How will the 8000 additional samples for each second be computed ?

I tried this using sox and a new recording has been generated with no complaints or error messages, So my question is : how is this done ? Is there some kind of standard procedure, like an interpolation used to build the missing samples  ?

The new samples are generated by interpolating between the original ones. Exactly how this is done will vary by implementation, but the most typical way would be to use a linear interpolating filter. With this technique, you would interpolate by a factor of 2 by inserting zeros between each of your input samples. Assuming your input signal is $x[n]$, your expanded signal would look like:

$$
x_e[n] = [ x[0], 0, x[1], 0, x[2], 0, \ldots ]
$$

Due to the properties of the discrete time Fourier transform, this has the effect of compressing and duplicating the spectrum of the input signal (so all of the spectral content that were present between DC and the Nyquist rate of the original signal are compressed into a region half as wide, then duplicated twice).

Finally, you then apply a lowpass filter to remove the extra copy of the spectrum above $f_s/4$ in the expanded signal. The result is an interpolated-by-2 version of your input signal that is bandlimited to the same region of frequencies as the original.

Doing this for an audio signal isn't likely to do much for you, unless you have some equipment or algorithm that requires a particular input sample rate. The interpolated audio, when played back at the new rate, should sound the same as the original signal; the interpolation process doesn't create any new information.

Edit: I should also note that the above is just a conceptual description of the process. In practice, you wouldn't use the explicit step of zero-stuffing followed by linear filtering. Instead, you can use multirate techniques like a polyphase realization of the interpolation filter that allow you to achieve the same effect with fewer computations.


Why is bit error rate = symbol error rate / number of bits per symbol in QPSK I saw a code computing the error probability of QPSK, and I found that in this code the bit error rate is computed as symbol error rate divided by $2$.

I think that $2$ is because the QPSK symbol consists of two bits. But why can the bit error rate be equal to the symbol error rate divided by the number of bits per symbol?

In general, the bit error probability depends on the symbol error probability and the coder, i.e., the mapping between bits and symbols. The formula you refer to

$$P[\textrm{bit error}]\approx\frac{1}{M}P[\textrm{symbol error}]\tag{1}$$

with $M$ being the number of bits per symbol ($M=2$ for QPSK) is an approximation for a certain type of coder under the assumption that the SNR is relatively high.

In $(1)$ it is assumed that the system implements a Gray code, where nearest neighbor symbols only differ by one bit. So if the SNR is high, most symbol errors occur due to confusing a symbol with one of its nearest neighbors, which causes only one out of $M$ bits to be wrong.


Ratio of the power in the sidebands to the total power in conventional AM Message signal and carrier signal are given as

$m=\cos(2 \pi  \cdot f_{m} \cdot t) $

$c=A_{c}\cdot \cos(2 \pi  \cdot f_{c} \cdot t) $

I was asked to find modulation index $a$ so that ratio of the power in the sidebands to the total power in the transmitted signal will be $\frac{1}{3}$

My attempt

We know that, in conventional AM, total power of transmitted signal $P_{u}= \frac{Ac^2}{2} + \frac{Ac^2}{2}\cdot a^2 P_{m} $  where $P_{m} $ is message signal's power.

I think that power in sidebands is related to $\frac{Ac^2}{2}\cdot a^2 P_{m}$ term

$\frac{Ac^2}{2}\cdot a^2 P_{m} = \frac{1}{3}P_{u}$

Submitting $P_{u}$

$\frac{Ac^2}{2}\cdot a^2 P_{m} = \frac{1}{3}(\frac{Ac^2}{2}(1 +  a^2 P_{m}))$

$ a^2 P_{m} = \frac{1}{3}(1 +  a^2 P_{m})$

$ a^2 =  \frac{1}{2P_{m}}$

It seems false because modulation index might have been a numerical value.

this is John BG jgb2012@sky.com

In AM modulation, as shown in

https://en.wikipedia.org/wiki/Amplitude_modulation#Modulation_index

the modulation index is the ratio between the amplitude of the modulated and the amplitude of the carrier.

So the sought modulation index $a=\frac{1}{Ac}$

modulating signal power: $\frac{1}{2}$

carrier power: $\frac{Ac^2}{2}$

side bands to total power ratio:   $\frac{1}{(1+\frac{Ac^2}{2})}$

solving:  $\frac{1}{(1+\frac{Ac^2}{2})} =\frac{1}{3}$ gives $Ac=2$

so the answer is $a=\frac{1}{2}$


Several questions about power of periodic signal I am currently learning Fourier series (periodic signals) and there are few things I am not sure about.

(1)

$$
P = \frac{1}{T}\int_{-T/2}^{T/2}|s(t)|^2dt = \frac{1}{T}\int_{0}^{T}|s(t)|^2dt
$$

Let's say, that period is $T = 2$ and my signal $s(t) = t$:

$$
P = \frac{1}{T}\int_{-T/2}^{T/2}|s(t)|^2dt = \frac{1}{2}\int_{-1}^{1}t^2dt = \frac{1}{3}
$$

but

$$
P = \frac{1}{2}\int_{0}^{T}|s(t)|^2dt = \frac{1}{2}\int_{0}^{2}t^2dt = \frac{4}{3}
$$

What am I missing here? Are these equations wrong? I've checked several sources and all of them state it this way. What kind of power do I calculate this way? Is it total power of periodic signal?

(2) I also have a problem, I don't understand how to calculate.
I am given:

\begin{align}
x(t) &= t\\T &= \frac{1}{3}
\end{align}
what is the power on interval
$$
\left(\frac{1}{3};\frac{2}{3}\right)
$$
How do I calculate this?

(3) Am I right if I say, that I calculate power of first $N$ harmonics with Parseval's theorem?

$$
P_n = a_0^2 + \frac{1}{2}(a_n^2 + b_n^2)
$$

I am really trying to understand it, but the more I read, the more confused I am. So I really need to point me in the right direction.

You are missing the effect of the integration interval on the integrand function...

If you choose a period; T = 2 , and a periodic signal whose base period is $s_0(t) = t$  in the interval [0,2]... Then the integrand in the shifted interval [-1,1] will be different as given by

$$
s_1(t) = \begin{cases}{ t + 2  ~~~, -1 < t < 0 \\ t   ~~~~~~~~~~ , ~   0< t < 1}\end{cases}
$$

Then you should write the integrals as follows: :
$$
P = \frac{1}{2}\int_{-1}^{1}|s_1(t)|^2dt = \frac{1}{2}\int_{0}^{2}|s_0(t)|^2dt
$$

You should get the same result now...


How can I correctly plot an impulse_response() of a discrete transfer function? I posted this question on stack overflow and I think it was the wrong place so I moved it here.

Hello, I'm new to python, I just want to use the control library and I'm slowly learning, I did my research and read the docs but as a newbie programmer I don't really understand much. I know it's about the syntax but I don't know what exactly is wrong. I'm moving from matlab to python and I can't get this to work:

from control import *
from sympy import *
from numpy import *
from matplotlib.pyplot import *

Gp = tf(100,[1,0,100])
T = 0.05
Gz =sample_system(Gp, T)

x = linspace(0.0, 0.7, num=15)

T, yout = impulse_response(Gz, T=x)
plot(T, yout)

if I run this I get this output:

ValueError: Time steps ``T`` must match sampling time

and if I change x = linspace(0.0, 0.7, num=15) for x = linspace(0.05, 0.7, num=14) I get this:

ValueError: A value in x_new is below the interpolation range.

No idea what I'm doing wrong, I read the docs on impulse_response() and linspace(), I can't find any examples of similar problems or people plotting impulse responses using python. Thanks in advance

Edit1: Okay I've gotten here so far:

imp = signal.unit_impulse(100, 'mid')
T, yout = signal.butter(4, 0.2)
z = array([0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5])
T, yout = impulse_response(Gz, T=z, input=imp)
T, yout = signal.butter(4, 0.2)
plot(T, yout)
grid(True)
show()

now I get a plot that looks like this:

But I need something like this:

Edit2: Okay seems like the whole problem was that yout had these dimensions: (1,15) and all I had to do was .flatten(), on the other hand I was using wrong linspace so I changed it for arange, and the end result is this:

Gp = tf(100,[1,0,100])
ts = 0.05
Gz =sample_system(Gp, ts)
pprint(Gz)

t = arange(0.0, 15*ts, ts)

T, yout = impulse_response(Gz, T=t)
yout=yout.flatten()
plot(T, yout)
margins(0.1, 0.1)
xlabel('Time [sec]')
ylabel('Amplitude')
grid(True)
show()

and the output is:

Still seems to me that something is missing, but I don't know what, and that this is not a clean solution, anyone got a cleaner way of doing this? Having to flatten() seems wrong to me... in matlab you just had to impulse(Gz, t) and got the plot right away.

So finally this is the way I got to plot it right:

Gp = ct.tf(100,[1,0,100])
ts = 0.05
Gz =ct.sample_system(Gp, ts)
sym.pprint(Gz)

t = np.linspace(0,14*ts,num=15)

T, yout = ct.impulse_response(Gz, t)
yout=yout.flatten()

fig, ax = plt.subplots()
ax.step(T, yout)
ax.yaxis.set_minor_locator(MultipleLocator(0.05))
ax.yaxis.set_major_locator(MultipleLocator(0.1))
plt.margins(0.1, 0.1)
plt.title('Impulse Response c.1)')
plt.xlabel('Time [sec]')
plt.ylabel('Amplitude')
plt.grid(True)
plt.show()

and this is the output:

this would be the equivalent of matlab's impulse() in python


Mathematically Inclined Signal and Systems / Signal Processing Book Recommendations I'm an electronics engineering student with high inclination to analysis and pure mathematics.  I was just wondering if there was any book ( or any resource ) that treats signal and systems and signal processing with a lot of mathematics rigour ( actually doing proper complex analysis, using functional analysis and linear algebra rigorously to explain convolution, fourier, laplace and z transforms for example ).

I'm very disapointed with the books i've read ( Oppenhein, Lathi and related ) because it actually throws a lot of the beauty of analysis and algebra away, focusing on the computational side.

Thanks a lot

Oppenheim and Willsky's Signals and Systems or Lathi's Linear Systems and Signals are intended for Sophomores who have only a single semester of differential equations under their belts, so it is a bit unfair to criticize them for leaving out the functional analysis and the conformal mapping.  At the sophomore level my favorite book is Siebert's Circuits, Signals, and Systems.  It won't give you the mathematical rigor you desire either, but you can see his great love for the mathematics and he has these wonderful, witty, footnotes that provide a great historical perspective.

There is a great book (which I love, but do not recommend to you) by the (applied) mathematician Richard Hamming (of the "Hamming window", "Hamming code", "Hamming distance", "Hamming bound" and "Hamming problem") called Digital Filters.  In it he makes a number of snarky comments like:

  Since we are interested mainly in using mathematics, we are obliged in our turn to be ambiguous with respect to mathematical rigor.  Those who believe that mathematical rigor justifies the use of mathematics in applications are referred to Lighthill and Papoulis for rigor; those who believe that it is the usefulness in practice that justifies the mathematics are referred to the rest of this book.  (1998 Dover edition, page 72.)

So in addition to the book by Papoulis that @Matt_L recommended I will add Hamming's (and Siebert's!) recommendation of M.J. Lighthill, Fourier Analysis and Generalized Functions, Cambridge University Press, 1958.


How does Siri recognize me saying "Hey Siri"? I am trying to understand how my iPhone can continually listening for me saying Hey Siri, Alexa, Hey Cortana or Okay Google without quickly draining my battery down.

I imagined two kind of algorithm. One that record slice of time such as 10 ms wide slices each 200 ms and perform a synchronous detection on specific frequencies. However these parameters strongly depend on the characteristic of my voice. Moreover, it will still consume a lot of CPU power to continually try to match a Hey Siri in the middle of nowhere.

What kind of efficient low-power algorithm/implementation (hardware or software) can perform such task?

This is somehow related to this patent I think: https://www.google.com/patents/US20160253997

I read some articles that talk about hidden Markov models, but I doubt it is very low power approach.

"Ok Google" is described in many publications by Google

Automatic Gain Control and Multi-style Training for Robust Small-Footprint Keyword Spotting with Deep Neural Networks

Convolutional Neural Networks for Small-Footprint Keyword Spotting

It is based on DNN specifically trained for keyphrase and runs really fast. It does not consume a lot of power even on mobile.

Alexa spotting is implemented by the same guys and available as Snowboy

UPDATE: Apple describes their implementation here.


Welch's overlapped method for amplitude by frequency estimation pwelch  is a useful matlab function for computing real/magnitudes for power using WOSA method, but I'd like to apply the same thinking in order to output complex pairs for amplitude (not power) spectra.

Does anyone know if this is possible? Or if there is an existing script that I might be able to modify?

Thanks!

I've found this:
Implementing Welch's method for Power Spectral Density

That seems to answer the question. However, the author encountered a problem, which I believe may be addressed by modifying the code that was identified as problematic in the thread by:

dataA(istart:iend,j) = data1(istart:iend,j).*window;
dataB(istart:iend,j) = data2(istart:iend,j).*window;


Decoding of binary signal from a temperature sensor I have the following signal from a temperature sensor (broadcasts on 433.92 MHz, apparently as ASK).

I cannot for the life of me extract the temperature data.

I've tried reading off the bits in so many different ways (mostly by looking at the pulse widths) then searching the data for the temperature info.

For example:

This corresponds to a temp of 21.0.

I tried to read each pulse after the first 4 big pulses - that gives 37 bits:

1000111110010010110001100011111111000

Then I tried searching just for the 21 part. So I tried searching for 21, 61, 71 (in case the designer was adding 40 or 50 degrees to avoid negative numbers). I tried searching for the 2 and then the 1 bit shifted by 4. I tried negating all the bits.  I tried searching in reverse.

Is it possible this is not pulse width modulation?  Is it possible there are not 37 bits here?

It is impossible to extract the correct data from this string with only one sample (as is mentioned by Laurent Duval in the comments). Here is something that caught my eye:

Markers

You have these four blocks at the beginning and the end of the signal that are some kind of start/end markers, I suppose. Are you sure that you caught the whole markers?

1000111110010010110001100011111111000 could be split up in the following way:

1000
1111

100100101100011000111

1111
1000

so the 1000 and 1111 might actually be part of the markers.

The temperature itself

You say the signal corresponds to 21.0°, i.e. a number with an accuracy of 0.1°. It could be that this is encoded as the integer 210. 210 in binary is 11010010 and you can find this string (reversed) in there:

100100101100011000111
  01001011

The prefixing 10 is 2 in decimal and could be a 2 bit indicator of the number of 4 bit words that encode the temperature, however this would be strange since the temperature is reversed and the prefix is not...

Please try the following:

Measure the same temperature multiple times in independent experiments. Is the resulting data string the same in each experiment? Are parts of the string the same in each experiment?
Measure multiple temperatures that only differ by 0.1°. Is there a corresponding pattern in the data string?
Do you have a datasheet of your sensor?


what is the relationship between the reliability of a monte carlo test and the length of the samples tested and the number of errors found? Considering this setup, single tone transmission over a channel impaired by AWGN, doppler spread, and inter symbols interference.
We are to compute the BER as a function of the SNR.
We want to find the number of samples that need to be tested and the number of errors that need to be found.
What number of samples/errors need to be simulated to observe a bit-error rate values ($10^{-3}$ for instance) that have 90% confidence level of being within a factor of two of the actual values?

For each sample, we note $1$ if error and $0$ otherwise. Then, for bit error rate $q$, the tests are independent Bernoulli random variables $\{X_i\}$ with probability $q$.
The estimate is $\hat{p} = \frac{S}{n}$  where
$$S = \sum^n_i X_i \tag{1}$$
and $\mu = \mathbb{E}[X] = np$.
Chernoff bound (see Corollary 5) tells us that for every $0 < t < \mu$,
$$\mathbb{P}\left( \mid S-np\mid \leq \frac{t}{\mu}\mu\right) \geq 1 - 2 e^{\frac{-t^2}{3\mu}} \tag{2}$$
Choose $t=\epsilon np$ with $0 < \epsilon < 1$,
$$\mathbb{P}\left( \mid S-np\mid \leq \epsilon np\right) \geq 1 - 2 e^{\frac{-\epsilon^2 np}{3}} \tag{3}$$
$$\mathbb{P}\left(p(1-\epsilon) \leq  S/n \leq p(1+\epsilon)\right) \geq 1 - 2 e^{\frac{-\epsilon^2 np}{3}} \tag{4}$$
Therefore, $\epsilon$ is relative estimation error and the right hand side of $(4)$ is confidence level which depends on sample size $n$.
For example, if you want estimate a BER within $10^{-3}\times(1 \pm 10\%)$, using sample size $n=8.99 \times 10^5$ will guarantee a confidence level greater than 90%.
Note that choosing another $t$ and/or using another bound may result in a tighter required sample size.


RADAR Video Signal and Matched Filtering In order to make my question more clear and give more details
I'm trying to understand the processing signal procedure in radars,i understood till this moment that a transmitter generates a signal called transmitted signal with a center frequency and then when this signal meets a target an echo signal reflected to radar and captured with a received stage including many RF and non RF components.After amplification and noise rejection with LNA a mixer is used to melange received and transmitted signal to produce a signal called video signal and then the signal processing stage take place.
My question are : can anyone explain me what's the role of the matched filter in radar signal processing stage what are the alternatives methods?

The simple answer is the matched filter increases the signal to noise ratio of the reflected return signal.

As the signal propagates from transmitter, the power that impinges on the target is proportional to $1/r_{target}^2$.  The reflected energy back to the transmitter is also proportional to $1/r^2$, so the received energy is proportional to $1/r^4$.

The received back echo power is very small and the matched filter maximizes the deflection of the echo with respect to the noise in the receiver.

This is a simple explanation. There are many types of RADAR, and there are additional ways to make the echo stand out with respect to the noise of the system and the clutter echos present in all real systems

Those additional ways that are specific to the different kinds of RADAR systems to make the echo discernible. There are various ways to modify a matched filter for a specific kind of applications but it would be hard to say they were alternative to a matched filter. The underlying optimality characteristic of a matched filter justifies few alternatives.


What does this block operator mean? What does this block do? It's part of a block diagram for a RMS/Peak detector.

It symbolizes the input-output characteristic of a half-wave rectifier, which only passes positive input signals:

$$y(t)=\begin{cases}x(t),&x(t)>0\\0,&x(t)\le 0\end{cases}$$


Defining bound on input signal to test accumulator for BIBO stability For an accumulator, defined as shown in the image below, why would I define $B_x=1$? $u[n]$ is defined at zero so my (possibly misguided intuition) is telling me that I'd choose $B_x = 0$ to not eliminate that point.

$B_x$ is a bound on the magnitude of the input signal. So if $x[n]=u[n]$ you have $|x[n]|\le1$ and, consequently, $B_x=1$ is a valid bound on $|x[n]|$. $B_x=0$ would mean that the input signal is zero.


Maximum Magnitude Deviation between DFT and DTFT Let $x[n]$ be a finite-length discrete-time signal with length $N$. The continuous DTFT $X(\omega)$ is then
$$
X(\omega) = \sum_{n = 0}^{N-1} x[n] e^{-j \omega n}.
$$

The length-$N$ DFT of $x[n]$ is
$$
X[k] = \sum_{n = 0}^{N-1} x[n] e^{-j 2 \pi \frac{k n}{N}}.
$$
For this, the DFT is a sampled version of the DTFT, i.e.,
$$
X[k] = X(2\pi k / N).
$$
Please also see posts here and here. Now, we consider the maximum magnitudes $m_{\textrm{d}} = \max_k |X[k]|$ and $m_{\textrm{c}} = \max_\omega |X(\omega)|$. Because of above, we have $m_{\textrm{d}} \leq m_{\textrm{c}}$.

What is the relation between these maximum values? Is there $\gamma > 1$ such that $\gamma m_{\textrm{d}} \geq m_{\textrm{c}}$?

okay, i'm gonna steal a diagram from this publication:

now this is an approximation that is based on the assumption that the discrete spectral points $|X[k]|$ are close enough to each other that only the discrete peak:

$$m_\textrm{d} = \max_k |X[k]| = |X[m]|$$

and the two adjacent neighbors $|X[m-1]$ and $|X[m+1]|$ need be considered. (these are shown as $k$ in the diagram above, not $m$, but i don't know if i should change the convention of discrete variables we're using and i don't wanna draw out another diagram.)

the assumption is, within this very tight range of $\omega$, that the peak is adequately approximated by a quadratic:

$$ |X(\omega)| = m_\textrm{c} - a(\omega - \omega_\textrm{m})^2 $$

where $a > 0$, and the "spectral peak":

$$\star \quad |X(\omega_\textrm{m})| = m_\textrm{c} = \max_\omega |X(\omega)|$$

so we have three unknowns: $|X(\omega_\textrm{m})| = m_\textrm{c}$, $a$, and $\omega_\textrm{m}$ and three equations by evaluating this quadratic at the discrete points $|X[m-1]|$, $|X[m]|$, $|X[m+1]|$:

$$\begin{align}
|X[m-1]| &= |X\left(2\pi \tfrac{m-1}{N}\right)| &= m_\textrm{c} - a (2\pi \tfrac{m-1}{N} - \omega_\textrm{m})^2 \\
|X[m]| &= |X\left(2\pi \tfrac{m}{N}\right)| &= m_\textrm{c} - a (2\pi \tfrac{m}{N} - \omega_\textrm{m})^2 \\
|X[m+1]| &= |X\left(2\pi \tfrac{m+1}{N}\right)| &= m_\textrm{c} - a (2\pi \tfrac{m+1}{N} - \omega_\textrm{m})^2 \\
\end{align}$$

so you can solve those three equations for the three unknowns.  we're not so interested in $a$, but we are interested in the frequency location of the peak $\omega_\textrm{m}$ and the height of the peak, $m_\textrm{c}$.  if you slug out solving this set of three equations and three unknowns, you get for those two parameters:

$$\begin{align}
\omega_\textrm{m} &= \frac{2 \pi}{N}\left( m + \tfrac12 \frac{|X[m+1]| - |X[m-1]|}{2|X[m]| - |X[m+1]| - |X[m-1]|} \right) \\
\\
|X(\omega_\textrm{m})|= m_\textrm{c} &= |X[m]| + \frac{\tfrac18 \big( |X[m+1]| - |X[m-1]| \big)^2}{2|X[m]| - |X[m+1]| - |X[m-1]|}   \\
\end{align}$$

Now this makes no assumptions about the nature of the spectrum except that the discrete spectral peak is at $|X[m]|$ and the neighboring spectral samples, $|X[m-1]|$ and $|X[m+1]|$ are no larger.  With no other information, this is about the best or only guess you can make.  It's not always perfect, but it should get you to within one-half FFT bin of the true DTFT peak.

It is perfect if you're looking at the peak of an isolated sinusoid, that was windowed with a gaussian window before going to the FFT, and you're looking at the dB or log magnitude.  Then this quadratic behavior in the log magnitude scale is mathematically the correct model and this quadratic peak interpolation is exactly correct.  So then you replace $|X[m]|$ and $|X[m-1]|$ and $|X[m+1]|$ with the logarithms (or dB values) of the same before using the above formulae to get the continuous-frequency peak.


Confused about Hamming window and its relation with frequency response specification $w(n)= 0.54 + 0.46\cos(2n \pi/(M-1))$ is Hamming window and we have say $f_{s_1}$,$f_{p_1}$,$f_{p_2}$,$f_{s_2}$ with pass and stopband tolerance specifications of the desired filter. From this specification can we get mathematical corresponding frequency response?

Windowing has nothing much to do with your filter design/response spec's, however, it has impact on your transition bandwidth that also depends on the window length.Your filters frequency response is mainly dependent on design parameters such as sampling rate, maximum frequency, transition bandwidth, cut-off frequency and impulse response length or kernel length.Now coming to windowing, these are some techniques basically used for obtaining smoother and ripple free impulse response h[n], so before your multiply your impulse response to the input/incoming samples, you must window your impulse response. For e.g. If you have a sinc impulse response, Ideally we know that its amplitudes run infinite on both negative and positive ends, hence in practice we just taper or cut off the ends of the sinc impulse in order to make it a finite length, therefore, this causes ripples in the impulse response due to the Gibbs Phenomenon (for more info on Gibbs Phenomenon, check online), So to eliminate these ripples and get a smoother impulse response, we use these windowing techniques.And each windows vary, for e.g. Blackman windows have better stopband attenuation than hamming window. Between why have you specified two stopband and passband frequencies?. Anyways,hope this gives you some insight to windowing!


Deriving the DFT magnitude of $A\cos(2\pi nk/N)$ Given that $$x(n) = A\cos(2\pi nk/N),$$ the $N$-point DFT of $x(n)$ can be expressed as follows—the derivation can be found in here: $$X(m) = \color{red}{\frac{A}{2}\sum_{n=1}^{N-1}e^{-j2\pi n(m-k)/N}}+\frac{A}{2}\sum_{n=1}^{N-1}e^{-j2\pi n(m+k)/N}.$$ It is said that the left summation above represents a positive-frequency spectral component. Does it mean to say that, when $m$ is negative, the summation goes down to zero?

So here's a plot of the summation.  It is true that $X(m)$ is zero for some values of $m$ negative. However, due to the $N$-periodic nature of the summation, any $m$ for which $k\ \mbox{mod}\ N = m$ will be non-zero.

EDIT

  I still don't understand how the negative portion of bins goes down to zero as a result of m−k

So, rewrite the summation as:

$$
\sum_{n=1}^{N-1}e^{-j2\pi nm/N} \cdot e^{j2\pi nk/N}
$$

so that all you are doing is modulating a complex exponential of frequency $m$ with one of frequency $k$ and summing the result.

This is just taking the inner product of two complex exponentials.  The inner product just projects one exponential onto another. In this special case (for discrete time) where the frequencies are harmonically related, the only time when this is non-zero is when $m = k \ \mbox{mod}\ N$.

R Code Below
    #30344

N <- 100
k <- 5
m_vals <- seq(-N,N)

Xm <- rep(0,2*N+1)

for (m in m_vals)
{
  Xm[m+N+1] <- sum(exp(-1i*2*pi*seq(0,N)*(m-k)/N))
}

plot(m_vals, Re(Xm), type="l", col="blue", xlab="m", ylab="X(m)")

title("Plot of X(m) vs m for m -100 to 100")


Circular convolution of a non causal signal I know how we compute the $N$ point circular convolution of a two causal signals, but what about a signal such as $\{1,-1,2,1\}$ where, the position of 2 is the $0^{th}$ index and the other sequence is $\{2, -1\}$ which we can assume  to be causal, what about the 4 point circular convolution. According to me it is

$$\begin{bmatrix}1&1&2&-1 \\-1&1&1&2\\2&-1&1&1\\1&2&-1&1 \end{bmatrix} \begin{bmatrix} 2\\-1\\0\\0\end{bmatrix} =\begin{bmatrix} 1\\-3\\5\\0\end{bmatrix}$$
With the position of 5 being the zeroth index because only then the 2 from the first signal got multiplied with the 2 of the second signal, giving off the zero position. But now I am confused, as to how to arrange the other indices. Can anyone help me out?

For an $N$-point circular convolution you can think of each signal as being periodically extended with period $N$. For your example with $N=4$ that would mean that the two sequences are

2 1 1 -1 and 2 -1 0 0

where both now start at index $n=0$. The result of the cyclic convolution is

5 0 1 -3

which is just a cyclic shift (by $2$) of the (correct) result that you obtained.


Average value of power spectral density I have a spectra of power spectral density (PSD) (vector (1:1e6)) from experimental data (blue color).

However I want to plot only average value of PSD (red curve).

How can I do it by using MATLAB?

Thank you so much !

You should remove low frequency content of the diagram. Simply use moving average filter. find the optimum filter size that produces your desired result. Here is a matlab example code.
t=0:0.001:1;
PSD=2*sin(2*pi*1*t)+10*randn(size(t))+5; %assume as the power spectrum
filtered=zeros(size(t));
M=20 ;% Moving average filter size
for i=1:length(t)-M
    filtered(i)=mean(PSD(i:i+M));
end
plot(PSD);
hold on
plot(filtered,'r')


How to remove lines from captcha image which are acting as noise using java? The images are attached here. I first converted the sample image to grayscale and then binarized it.

Finally i applied median filter to remove the dots from the image. However i am unable to remove the lines which are acting as noise in the image. Kindly help me out.

Original Image:

Processed Image:

I worked on this problem some time ago, my solution was to compress the image with wavelet transform with an opportune coefficient that depends on captcha type that I choose with some expermients.
Then I saturate the compressed image from mean value to 0 and 1 in order to create a mask.
Next I multiply mask to the original image and I get the image without lines or noise. I post a picture example where I use also a derivative filter on the result but in reality you don't need it.

Your lines are a bit different but I think that you can start from this approach because your lines are more thin than text so wavelet compression can works but you have to improve it for your problem.

Finally your lines have different color and you can use also this information.


breaking the input signal into segments while calculating DFT Im reading a book about DSP and there is an example about investigating the sound that travel through the ocean. The sound was recorded as the time-domain signal and they tried to analyze it in frequency domain.

The input signal was broke up into 256 point segments. Each of these segments is multiplied by the Hamming window, then run through a 256 point DFT. The resulting frequency spectra are then averaged to form a single 129 point frequency spectrum. Please take a look at the figure

I don't understand:

How they break up the input signal into 256 segments?
How each segment was multiplied with the Hamming window?
How they calculate 256 DFT points?
How the frequency spectra was averaged?

How they break up the input signal into 256 segments?

This is usually done by a sliding a window of a given width (here 256). Depending on the step-size this results in overlapping (if smaller than 256) or non-overlapping (if equal or larger than 256) data segments.

  How each segment was multiplied with the Hamming window?

Usually, by an element-wise multiplication.

  How they calculate 256 DFT points?

The number of DFT points is, if I am not mistaken, determined by the length of the input signal.

  How the frequency spectra was averaged?

All resulting spectra were summed up and divided by the number of spectra.


Energy Detection Algorithm Uncertainty I'm currently working on a project to perform energy detection of an IEEE 802.15.4 using the block diagram illustrated below. However, I do have a little uncertainty about the algorithm. I'm wondering after finding the magnitude squared values at each FFT bin, do you then add each bin value up and divide by the number of bins (N) to get the average power? And also, do you just repeat this averaging M times? And why exactly do you repeat this averaging M times??

I'm wondering after finding the magnitude squared values at each FFT bin, do you then add each bin value up and divide by the number of bins (N) to get the average power?

Yes.

  And also, do you just repeat this averaging M times?

Yes, you would have to...

  And why exactly do you repeat this averaging M times??

...to reduce the impact of noise on your measurement and also, depending on the length of the FFT, to make sure that you have captured the whole dynamics of your signal.

For more information, please see Power Spectrum Estimation Using the FFT.

Hope this helps.


Calculating cutoff frequency for Butterworth filter I have a problem while calculating cutoff frequency, suppose we have these specs.

Firstly, I calculated the order of the filter and got $N=5.8858$ and round it up to get $N=6$.

Now I'm supposed to get $\Omega_c$. Using these equations:

\begin{cases}1+\left(\frac{0.2\pi}{\Omega_c}\right)^{2N} = \left(\frac{1}{0.89125}\right)^2 \quad&\quad (1)\\
1+\left(\frac{0.3\pi}{\Omega_c}\right)^{2N} = \left(\frac{1}{0.17783}\right)^2\quad&\quad (2)
\end{cases}

Now with $N=6$ and $T=1$, substituting in $(1)$

\begin{align}
\left(\frac{0.2\pi}{\Omega_c}\right)^{12} &=\left(\frac{1}{0.89125}\right)^2 - 1 =0.25893\\
\implies (\Omega_c)^{12} &= \frac{{(0.2\pi)}^{12}}{0.25893}\\
 &= 40.29
\end{align}
But in the textbook it says $\Omega_c = 0.7032$, what I did wrong? Any help would be appreciated.

Everythng is alright. You made a mistake in the final calculation.

$$\Omega_c^{12} = \frac{{(0.2\pi)}^{12}}{0.25893} \implies \Omega_c = \sqrt[12]{\frac{{(0.2\pi)}^{12}}{0.25893}} =\frac{0.2\pi}{\sqrt[12]{0.25893}} = 0.7032$$

which is the desired result.


Confused in Bit rate, Symbol rate, Sample rate and Data rate I'm using GMSK modulation to send $300 \ \rm bits$ and assuming $1 \ \rm bit/symbol$. $BT$ product is $0.5$. Sampling frequency is $100\ \rm kHz$ and bandwidth of the signal is $28 \ \rm kHz$.

(a) Is bit rate calculated as $BT/B = 0.5/28\rm k$ ?
(b) Here, symbol rate = bit rate ?
(c) The sample rate of input will be equal sampling frequency i.e. $100\ \rm kHz$?
(d) How can I calculate data rate ?

The $BT$ product is linked to the single-side $-3\rm dB$ down bandwidth and the bit rate $R_b$ with the relation below:
$$BT = \frac{f_{-\rm 3 dB}}{R_b}$$
Assuming you have a single-side $-3\rm dB$ bandwidth, you can extract your bit rate as:
$$R_b = \frac{f_{-\rm 3 dB}}{BT}$$
Also see the answer to this question. For more on GMSK see this paper $^1$ and this paper $^2$.

You have one bit per symbol, meaning for every bit you have one symbol. And their rate should be in that proportion as well. Or see it this way,
\begin{align}
\text{For} \ \ &n \ \ \text{bits you have} \ \ n \ \ \text{symbols}\\
\implies \text{for} \ \ &\Delta n \ \ \text{bits you have} \ \ \Delta n \ \ \text{symbols}\\
\implies \text{for} \ \ &\frac{\Delta n \ \rm bits }{\Delta t} \ \ \text{you have} \ \ \frac{\Delta n \ \rm symbols}{\Delta t}
\end{align}

See the sampling rate as the rate at which you are having your samples, and this value is in $\rm Hz$, or basically
$$
\frac{\Delta n \ \rm Samples}{\Delta t}
$$
The $\Delta t$s above are normally given as $1 \ \rm second$. The sampling rate is the number of samples per second, also known as sampling frequency (in $\rm Hz$).

Generally speaking, your bit rate is your data rate. The answer to this question details more on this.

You can use the maths and plug in your numbers.

$[1]$: A. Linz and A. Hendrickson, "Efficient implementation of an I-Q GMSK modulator," in IEEE Transactions on Circuits and Systems II: Analog and Digital Signal Processing, vol. 43, no. 1, pp. 14-23, Jan 1996.

$[2]$: K. Murota and K. Hirade, "GMSK Modulation for Digital Mobile Radio Telephony," in IEEE Transactions on Communications, vol. 29, no. 7, pp. 1044-1050, Jul 1981.


Applying Photoshop's "Shadow / Highlight" Correction Using Standard Image Processing Algorithms Adobe's Photoshop has an excelent filter called "Shadow/Highlight" used to correct, in one step, dark areas AND overexposed areas of a photo. By simple analysis, I suppose it uses some radius-dependent adaptive contrast algorithm (halo effect is typical), possibly combined with interpolation (even extra-dark or extra/white, clipped areas get very improved with this filter.

The parameters used can be seen in the dialog box below

What I would be looking for is a way to start exclusively with the parameter values given in the dialog box, apply a sequence of parameterized operations "blindly", and get back the result non-interactively (well, perhaps previewing the image while adjusting parameter sliders in some GUI environment).

Any suggestion about how to perform this with standard, well-known image-processing algorithms?

I think this can be achieved through histogram stretching. Input pixel values (levels) are mapped using linear transform from $<\textbf{min}_{in},\textbf{max}_{in}>$ to $<\textbf{min}_{out},\textbf{max}_{out}>$. The intensity transform for input and output pixel values look like this:

$$I_{out}=\frac{I_{in}-\textbf{min}_{in}}{\textbf{max}_{in}-\textbf{min}_{in}}\cdot(\textbf{max}_{out}-\textbf{min}_{out})+\textbf{min}_{out}$$

The formula can be simpler, but this one is more clear, at is transform input values from one interval to the other.

By default, $\textbf{min}$ is 0 and $\textbf{max}$ is 255 (depending on bit depth). Increasing value of $\textbf{min}_{in}$ emphasizes shadows, while decreasing $\textbf{max}$ emphasizes lights.

There are also "radius" track bars, so I suppose the histogram operation can be done adaptively / locally, i.e. histogram is computed for a window at every pixel instead of whole image.

Gaussian blur on work image prior to histogram computation may be convenient when the window is small to avoid sharp transitions.

Similarly, you can compute difference map which you smooth and then apply.

The adjustment will likely disturb color balance, so applying to luma channel only or computing adjustment of gray value and applying to R,G and B channels should suffice.


Can we get the same result from between partial FFT? I'm trying to implement fft and gaussian low pass filter multiplication to 16x16 image.

As I know fft can be separable. What if there are 16x16 one image, then is this the same ? between 4x4 part of the whole 16x16 fft result and  just only 4x4 fft result?

If they are same then I think that the latter method is more resource efficienciable.

I don't think "separable" means what you think it means – it means that you can implement the 2 DFT as the concatenation of a $N_{rows}$ amount of $N_{columns}$ long DFTs followed by $N_{columns}$ DFTs of length $N_{rows}$.

So, no, your approach is not correct.


Relation between original points and 1st/2nd derivative points I have points $\{p_0, p_1, ... p_n\}$, I create a discrete derivative consisting of $\{d_0, d_1, ..., d_{n-1}\}$ like this: $d_k=p_{k+1}-p_k$. I'd like to choose a point in the original signal by exploring the derivative. But after I find a specific point $d_k$, which point from the original set should I choose: $p_k$ or $p_{k+1}$? It looks like I can't really choose between them.

Similarly I explore 2nd derivative for inflection points. But in 2nd derivative I'm guessing that $p_{k+1}$ is a sensible choice because it's a midpoint out of 3 participants ($p_k, p_{k+1}, p_{k+2}$). Though this is still not perfect because of the 1st and last points in the set.

What are possible strategies to resolve this ambiguity?

Numerical signals are discrete and finite. So you are likely to hit issues either on end-points or discrete locations. I won't treat end-points here.

If you use only one operator, many signal processing folks care relatively little about a shift, as long as it is known, and constant. Here, the two-point derivative yields a $0.5$-point shift. You can consider that the derivative is valid at some mid-point. And you can interpolate it an integer locations if needed. Smoothing first can help limit a too-high fluctuations of differences.

A classical option is to choose an odd-length derivative operator. The three-point derivative writes:

$$ d_3[k] = (p[k+1]-p[k-1])/2$$

and gives you an information at the center of the schemes. Indeed, it coincides with the average of the left- and the right two-points skewed differences:  $ d_2[k] = p[k]-p[k-1]$ and $ d_2[k+1] = p[k+1]-p[k]$, since:

$$d_3[k] = (d_2[k] +d_2[k+1])/2$$

Many domains uses such schemes, and name them differently. To better find resources, some uses other terms, like generic ones, like "Finite difference coefficients", or specific ones, like "five-point stencil" or $n$-point gradient/Laplacian. One example is:

$$ d_5[k] = (p[k+2]- 8p[k+1] + 8p[k-1]-p[k-2])/12$$

The above thinking extends to higher order derivatives. One can further constraints methods using data properties, knowledge of noise, additional penalties, etc. To start with,  some literature on numerical derivatives:

CSERD Numerical Differentiation
Finite difference coefficients


Performance of Viterbi detector over non-minimum phase channels For sequences that are transmitted over channels with memory $\mu=n$ and response H=$[h_0 h_1 \ldots h_n]$,Viterbi algorithm implements Maximum Likelihood (ML) detection and BCJR implements Maximum A-Posteriori (MAP) Detector. Trellis over which these detectors operate have total states = $2^n$.

While applying any trellis based detector either ML or MAP detectors, channels are always considered causal i.e. $h_i=0 \quad \forall i<0$. More over equalizers/detectors are tested for channels that have $h_0=\max(h_i)\quad \forall i=0, \ldots,N$.

My question is related BER vs SNR performance of these trellis based detectors (that jointly equalize and detect) when applied to channels that are not causal and have a peak at $t=k,k\neq 0$. Since a delay to shift response to right has no effect on the structure of trellis, one can easily achieve this. But now the peak is no more at $t=0$. Following is an example to make it more clear.

Let $h_{ij}$ represent the $j^{th}$ coefficient for channel $i$. Consider two channels,

\begin{align}H_1 &=\left[h_{10} \space h_{11} \space h_{12}\right]=[4 \space 2 \space 1]\\H_2&=\left[h_{20} \space h_{21} \space h_{22}\right]=[2 \space 4 \space 1]\end{align}

The performance of ML and MAP detector on $H_1$ is better than that of $H_2$ which is obviously because of high ISI in $H_2$ as compared to $H_1$.

What transformation on $H_2$ can be done so that the resulting channel has a BER vs SNR that is close to $H_1$?

The Viterbi algorithm for optimum detection in the presence of intersymbol interference (ISI) is always preceded by a whitened matched filter (cf. this article by G.D. Forney). So if you use - as mentioned in a comment - the equivalent discrete-time channel, which then also includes the whitened matched filter (WMF), then, due to the properties of the WMF, this model must always be causal and minimum-phase. So, regardless of the actual impulse response of the channel, the complete equivalent discrete-time model (including the WMF) can never be non-causal or non-minimum-phase if implemented correctly.


estimating spectral optimization I'm relatively new to DSP so excuse my simplified words, and my detailed explanation.

if the signal have non-coherent sinusiod, it will induce energy spreading into the frequency domain.

One of the solutions is :

how to estimate the signal component correctly?

And I'm stuck here, I don't really know how to optimize such problem and setting a stopping criterion in Matlab,R..etc

How to optimize such objective function?

What are you thoughts about the drawbacks about such solution to separate signal components?

Thanks in advance :) ,

If you know the number of sinusoidal components, parametric estimators such as MUSIC may be suitable. ESPRIT might be another.  Both seem similar to you proposed optimization/minimization search.

But parametric estimators are reported to produce bogus decompositions if given the wrong number of components.  Possible reason: For a decomposition into an arbitrary number of arbitrary frequency sinusoids, there can be a near infinite number of solutions that might be arbitrarily close (e.g. you can get a bunch of cosines to converge to one sine wave of yet another frequency.)


Doubt on behaviour of frequency responce of a system I am doing loop shaping using $H_{\infty }$ approch, and I have a plant with an unstable pole:

G = 10/((s+10)*(s-1));
K1 = 4*(s+1)/s;

where G is the plant, and K1 is the controller. I am trying to do a confrontation between the approch with $H_{\infty }$ and the one without, in particular I am placing a weight only for the sensitivity function, and designing a controller using the mixed sensitivity design.

What I get if I plot the frequency response of the system in the two cases is the following:

in which the yellow line is the frequency response of the system in which I applied the $H_{\infty }$ and the red line is the one in which I used G and K1 as they are, so without mixed sensitivity and weights.

What I don't understand is why the yellow line has this strange behaviour, while the red line decrease smoothly.

I see that the bandwidth of the yellow line is lower the the other one, so maybe this fact is due to the fact that I have a RHP pole which imposes a lower limitation on the bandwidth of the system. But, if I plot the sensitivity functions of these two systems, I have:

so the yellow line behaves better with respect to disturbance rejection with respect to the red line, as expected.

I don't understand how to interpret the behaviour of the yellow line of the complementary sensitivity function.

Can somebody help me?

It appears to have -40 dB/decade roll-off in magnitude starting at 1 rad/sec and transitioning to -120 dB/decade roll-off in magnitude at approximately 300 rad/sec (this suggests you have 2 dominant poles that begin to be dominant when f is at 1 rad/sec and 4 additional poles that begin to contribute at f> 300 rad/sec since each pole would contribute -20 dB/decade to frequency roll-off.

In comparison, red is rolling off at -80 dB/decade consistent with 4 poles.  The poles would actually be complex given the frequency peaking in all cases.


How are Object Recognition and Content-based Image Retrieval interlinked? I need help on understanding the linkage between Object Recognition and Content-based Image Retrieval. Does OR use CBIR to find similar objects with the help of algorithms like SIFT etc or does CBIR use OR to improve its results on finding images showing the same object?

CBIR could use color histograms to find similar images but that has not too much in common with or so CBIR can stand alone i guess. OR needs to retrieve information which could come from learned images thus I suppose OR needs CBIR in some way. That would lead to the conclusion that OR is dependent and builds on CBIR, is that correct?

The difference is not so much in the algorithms used, but in the objectives.  Object recognition is a problem of naming an object or objects depicted in an image.  Content-based Image Retrieval is the problem of finding images in a database which match a user's query, which may be in the form of an image or in form of text.

I would say that CBIR is a more general problem than object recognition.  If anything, it is CBIR that builds on object recognition, and not the other way around.  If you have an image retrieval system, and the query is an image of a cat, then your system would need to recognize that this is an image of a cat (object recognition), and then find other images of cats in its database.  Presumably to be able to do that it would have to have performed object recognition on all of its images in advance, and tagged the ones that depict cats accordingly.

Note, that you can also have a query that may not require object recognition.  For example you may ask the system for images of urban landscapes, or images of indoor scenes.  This is a problem of "scene recognition" and it is typically approached  by using  some sort of a global image descriptor (e. g. scene gist), rather than trying to find all the object in the image.

Also note that an object recognition system must give you a single answer, i. e. the name of the object in the given image, and that answer is either right or wrong.  On the other hand, an image retrieval system gives you multiple possible matches for your query, and you would consider it successful even if only some of the top matches are good.  So there is a big difference in how you might evaluate an object recognition system vs. an image retrieval system.


How range of the signal increases when modulation is done? I am very confused because we know by increasing frequency we decrease the wavelength hence the distance covered will be very less but when we modulate a signal we are basically increasing the signal frequency which increases the range of the signal  , so how is that possible??

I think that you are interpreting the Friis equation to mean that RF signals carry for less distance than lower frequency RF signals. The “path loss” term “predicts this”.

The first point which is still confused of why we modulate should be fairly obvious, because we don’t want a 5000 meter antenna. One also uses modulation with amplification.  One also needs to be aware that frequency spectrum is regulated. In the US the FCC is not tolerant of unlicensed transmission so no signal regardless of carrier frequency will carry if you are in jail.

The Friis equation covers lossless spherical spreading just like the homogeneous wave equation.

https://en.m.wikipedia.org/wiki/Wave_equation

The wave equation predicts that the signal amplitude from an ideal point source drops off as $1/r$ and if that ideal point source has managed to radiate a signal there is no frequency dependence on the electromagnetic field.

The Friis equation doesn’t predict electromagnetic field behavior. It predicts the voltage at the output of an antenna. Path loss refers to the path loss with respect  to an antenna with a fixed aperture where voltage does drop off with the carrier frequency.

For someone like myself who started out as a SONAR engineer, the term “path loss” is an unfortunate choice.

Terrestrial RF propagation almost never conforms to the physical assumptions which the Friis equation is based on. It does have `relevance because it is actually useful as a first order heuristic.  Higher frequencies are more susceptible to scattering and shading than low frequencies.
Low frequencies penetrate buildings deeper. High frequencies tend be favored in line-of-sight applications.

There are other reasons why we want to use a modulated carrier. There isn’t a deep theoretical contradiction in play.

Essentially your conception of Electromagnetic attenuation is a bit too simple. Waveguides for instance have a characteristic known as the cut off frequency. A field has to have a frequency above the cut off frequency to propagate. Photons are electromagnetic waves and they travel billions of miles. The wavelength is not in itself the cause of attenuation, it is the interactions with matter which have the various effects


Signal amplitude to FFT amplitude I am using Python to generate a sine wave in order to cancel out part of a signal I'm analyzing. I would like the magnitude of this signal in the FFT to be equal to the magnitude of the signal that I want to get rid of. My questions are:

A. How does the amplitude of a signal translate to the amplitude of the resultant FFT?

B. Is there a way I can easily control the amplitude of an FFT using the amplitude of an input signal?

I can post my code if it helps.

This is the FFT of the signal I'm analyzing. I want the large amplitude at 1.4-ish(from the signal I generated) to be equal in magnitude to the large amplitude at 2.9-ish (the frequency I want to get rid of)

Consider a continuous-time signal $$x_a(t) = A \cos(\Omega_0 t + \theta) $$ where $\Omega_0$ is the frequency, and $\theta$ is the phase.
To perform spectral analysis of $x_a(t)$, sample it with a period of $T_s$ resulting in the discrete-time sequence $$x[n] = A \cos( \omega_0 n + \theta)$$ where $\omega_0 = \Omega T_s$ is the discrete-time frequency in radians.
Consider the relationship between the CTFT $X_a(\Omega)$ of $x_a(t)$ and the DTFT $X(e^{j\omega})$ of $x[n]$ for $\omega \in [-\pi,\pi]$
$$X(e^{j\omega}) = \frac{1}{T_s} X_a(\frac{\omega}{T_s})$$
Computation of $X(e^{j\omega})$ of $x[n]$ may require indefinetely long duration of $x[n]$, such as infinite duration if $x_a(t)$ is an ideal sine wave. And also $X(e^{j\omega})$ is a function with continuous domain of $\omega$. Thus it's not practical to compute the DTFT.
However, one can compute samples $X[k]$ of $X(e^{j\omega})$, known as the DFT of $x[n]$, defined over a finite length of $x[n]$ obtained by applying a window: $$v[n] = w[n]x[n]$$ where $w[n]$ is the window of length $L$.
The consequence of this windowing on the result of the computed spectrum is: $$V(e^{j\omega}) = \frac{1}{2\pi} X(e^{j\omega}) \star W(e^{j\omega}) $$
Computing the $N$-point DFT of the sequence $v[n]$ therefore provides the samples $V[k]$ of the result of the periodical convolution between the true spectrum $X(e^{j\omega})$ of $x[n]$ and the DTFT $W(e^{j\omega})$ of the window $w[n]$. Note that $N \geq L$.
For the specific example we have considered, $$x[n]=A\cos(\omega_0 n + \theta)$$ whose DTFT is $$X(e^{j\omega})= A\pi e^{j\theta} \delta(\omega-\omega_0) + A\pi e^{-j\theta} \delta(\omega+\omega_0)$$ a pair impulses weighted by $A\pi$.
And assuming a rectangular window of length $L$ for $w[n]$, its DTFT is: $$W(e^{j\omega})= \sum_{n=0}^{L-1} w[n] e^{-j\omega n} = e^{-j{\omega (L-1)/2}} \frac{\sin(\omega L / 2)}{\sin(\omega / 2)} $$ the result is $$V(e^{j\omega}) = 0.5 A e^{j\theta} W(e^{j(\omega - \omega_0)}) + 0.5 A e^{-j\theta} W(e^{j(\omega + \omega_0)})$$
Finally the DFT samples $V[k]$ are given by the following:
$$V[k] = V(e^{j\frac{2\pi k}{N} }) = 0.5 A e^{j\theta} W(e^{j(\frac{2\pi k}{N} - \omega_0)}) + 0.5 A e^{-j\theta} W(e^{j(\frac{2\pi k}{N} + \omega_0)}) $$
One peak peak of these DFT samples would occur for $\omega = \omega_0$ for some $k$ if $\frac{2\pi k}{N} = \omega_0$ is satisfied for an integer $N$ and $k$. Assuming $T_s$ ,$L$,$N$,and $\Omega_0$ are selected such that the equality is satisfied for some integer $k_0$, then the magnitude of the peak of the observed spectrum is:
$$ |V[k_0]| = | 0.5 A e^{j\theta} W(e^{j(0)}) + 0.5 A e^{-j\theta} W(e^{j(2\omega_0)}) | $$
Solving for $A$, the amplitude of the original analog signal whose spectrum analysis are after, yields:
$$A = \frac{2 |V[k_0]|}{|e^{j\theta} W(e^{j(0)}) + e^{-j\theta} W(e^{j(2\omega_0)}) |} $$
An approximation to this can be obtained by assuming, for the rectangular window, $W(e^{j 2\omega_0}) \approx 0$ and noting that $W(e^{j0}) = L$ : $$A \approx \frac{2 |V[k_0]|}{L}$$
If a window other than rectangular is used, formulation should be adjusted accordingly, specifically $W(0)$ being the sum of the window samples, which is $L$ for the rectangular window.


Efficient implementation of 2-d circularly symmetric low-pass filter How can a two-dimensional ideal circularly symmetric low-pass filter or its approximation be efficiently implemented on data sampled on a square grid? I'm referring to an ideal filter with a spatial frequency response that equals $1$ inside radius $\omega_c$ (the cutoff frequency) and zero outside it, and with ideal impulse response:

$$h[x,y] = \frac{\omega_c}{2\pi \sqrt{x^2 + y^2} } J_1 \big( \omega_c \sqrt{x^2 + y^2} \big).$$

A "not efficient" approach is to convolve with a 2-d filter kernel that is the product of $h[x, y]$ and a circularly symmetric window function.

Approximation by the real part of a weighted sum of separable complex Gaussian component kernels

Figure 1. The proposed scheme illustrated as 1-d real convolutions ($*$) and additions ($+$), for cut-off frequency $\omega_c = \pi/4$ and kernel width $N=41$. Each the upper and the lower half of the diagram is equivalent to taking the real part of a 1-d horizontal and a 1-d vertical pass separable 2-d complex convolution. Except for truncation/cropping, the two final components summed and the resulting kernel are guaranteed circularly symmetrical due to use of complex-number-weighted complex-variance Gaussian 1-d kernels.

I have suggested this approach for this problem, but not tried it before. A discrete-space 2-d convolution kernel (impulse response) $h[x, y]$ that is circularly symmetric (see Eq. 4) about $[0, 0]$, can be approximated to arbitrary accuracy by the real part of a separable (see Eq. 3) complex number weighted sum of Gaussian functions of complex variance:

$$h[x, y] \approx \tilde h[x, y],\tag{1}$$

$$\begin{align}\tilde h[x, y] = &\sum_{m=0}^{M-1}\operatorname{Re}\left(a_m^2e^{b_m(x^2+y^2)}\right)\tag{2}\\
= &\sum_{m=0}^{M-1}\operatorname{Re}\left(a_m e^{b_m x^2}\times a_m e^{b_m y^2}\right)\tag{3}\\
= &\sum_{m=0}^{M-1}\operatorname{Re}\left(a_m^2 e^{b_m r^2}\right)\tag{4}\\
= &\sum_{m=0}^{M-1}e^{\operatorname{Re}(b_m)r^2}\Big({\big(\operatorname{Re}(a_m)^2 - \operatorname{Im}(a_m)^2\big)\cos\big(\operatorname{Im}(b_m)r^2\big)\\
- 2\operatorname{Re}(a_m)\operatorname{Im}(a_m)\sin\big(\operatorname{Im}(b_m)r^2\big)\Big),}\tag{5}\end{align}$$

where $x$ and $y$ are integer horizontal and vertical coordinates, $r = \sqrt{x^2 + y^2}$ is the radius or distance between $[0, 0]$ and $[x, y]$, $M$ is the approximation order, and $a_m$ and $b_m$ are complex amplitude and variance related coefficients.

Optimized approximation, $\omega_c = \pi$

Eq. 5 with $r^2 = x^2 + y^2$ can be numerically optimized by global optimization to approximate a given $h[x, y]$ in spatial domain within the square $-(N-1)/2 \le x \le (N-1)/2$, $-(N-1)/2 \le y \le (N-1)/2$ with odd width $N$ of the square kernel. It suffices to do the optimization over a triangle $0 \le y \le (N-1)/2$, $y \le x \le (N-1)/2$. If the cost function is a sum of squares of errors or another sum of errors at each $(x, y)$, then cumulative $2\times$ error multipliers should be applied for each satisfied condition $x > 0$, $y > 0$, $x > y$.

The following C++ program optimizes the approximate kernel $\tilde h$ in spatial domain least squares sense using Differential Evolution from an optimization library. The target kernel $h$ is windowed using a 2-d rotated cosine window. The source code includes pre-optimized approximations for cut-off frequency $\omega_c = \pi$, kernel size $N=41$, and approximation orders $1 \le m \le 5$.

// -*- compile-command: "g++ -Wno-unused-result -march=native -O3 -ffast-math optitest.cpp opti.cpp" -*-
// This file is optitest.cpp
#include <stdio.h>
#define _USE_MATH_DEFINES
#include <math.h>
#include "opti.hpp"
#include "keyboard.h"

double optimal0[4] = { // Dummy starting point for optimization
  0,0,0,0
};

double optimal1[1*4] = {//bestcost=0.03593696209385285195 Least squares
  0.90954095944946312,0.21237789645057972,-0.72663606661376057,0.76147736559107793
};

double optimal2[2*4] = {//bestcost=0.01072419992771576382 Least squares
  -0.14935266179235707,0.28729650415029556,-0.08227954155725942,0.36475170816661134,1.03225137134955114,0.47158356759095016,-0.60039465413238513,0.47344404338750434
};

double optimal3[3*4] = {//bestcost=0.00322191268986909149 Least squares
  -0.00646998371362690,-0.15983965263134517,-0.02321635125852370,0.20289085505437962,-0.30298121364781033,0.34395623806288950,-0.10323939264492392,0.32564190139614663,1.13148495464950427,0.66639187966234981,-0.54567361845465989,0.37519342198291905
};

double optimal4[4*4] = {//bestcost=0.00089888657916611925 Least squares
  0.04481214216545543,0.08998210660870602,-0.01130470880820374,0.14003687599487286,0.02506113778123746,-0.24720956687380186,-0.03419309947921907,0.19643391040650712,1.27604032320787875,0.91159589755622883,-0.49585695598458995,0.29480666144138823,0.48259638667856242,-0.42280882153371496,-0.12730997070336811,0.29988730534029784
};

double optimal5[5*4] = {//bestcost=0.00017259634849055045 Least squares
  -0.06413980110992069,-0.03272558679644168,-0.00736621171073370,0.10943881870260203,0.08213936888117918,0.16226536153011967,-0.01975126456281387,0.13958603391531316,-0.07270617897425770,0.37473199916953354,-0.04702984945995840,0.18925309227383197,1.43523098843984531,1.17056740170289952,-0.45755771004055446,0.24332707322808175,0.69700596730108921,-0.48854609666142051,-0.14144974773647198,0.27672159791886242
};

const int numOptimalKernels = 5;
double *optimalKernels[numOptimalKernels + 1] = {optimal0, optimal1, optimal2, optimal3, optimal4, optimal5};

class CircularLPProblem : public Opti::Problem {
private:
  double *minimum;
  double *maximum;
  int numComponents;
  int numSamples;
  int *r2s;
  int *multipliers;
  double *target;
public:

  int getNumDimensions() {
    return numComponents*4;
  }

  double *getMin() {
    return minimum;
  }

  double *getMax() {
    return maximum;
  }

  double costFunction(double *params, double compare) {
    for (int m = 0; m < numComponents; m++) {
      params[4*m + 2] = -fabs(params[4*m + 2]);
      params[4*m + 3] = fmod(fabs(params[4*m + 3]), M_PI);
    }
    for (int m = 0; m < numComponents - 1; m++) {
      if (params[4*m + 3] > params[4*(m + 1) + 3]) {
        for (int k = 0; k < 4; k++) {
          double temp = params[4*m + k];
          params[4*m + k] = params[4*(m + 1) + k];
          params[4*(m + 1) + k] = temp;
        }
      }
    }
    double cost = 0;
    for (int k = 0; k < numSamples; k++) {
      double r2 = r2s[k];
      double trial = 0;
      for (int m = 0; m < numComponents; m++) {
        //  trial += exp(params[4*m + 2]*r2)*(params[4*m + 0]*cos(params[4*m + 3]*r2) + params[4*m + 1]*sin(params[4*m + 3]*r2)); // Max absolute error
        trial += exp(params[4*m + 2]*r2)*((params[4*m + 0]*params[4*m + 0] - params[4*m + 1]*params[4*m + 1])*cos(params[4*m + 3]*r2) - 2*params[4*m + 0]*params[4*m + 1]*sin(params[4*m + 3]*r2)); // Least squares
      }
      /*      if (fabs(trial - target[k]) > cost) { // Max absolute error
              cost = fabs(trial - target[k]);
              }*/
      cost += (trial - target[k])*(trial - target[k])*multipliers[k]; // Least squares
      if (cost > compare) {
        return cost;
      }
    }
    return cost;
  }

  // numComponents = number of components
  // N = kernel width
  // omega_c = cutoff frequency (radians)
  // Parameter vector: Re(a_0), Im(a_0), Re(b_0), Im(b_0), Re(a_1), ...
  CircularLPProblem(int numComponents, int N, double omega_c = M_PI): numComponents(numComponents) {
    numSamples = 0;
    for (int y = 0; y < (N-1)/2 + 1; y++) {
      numSamples += (N-1)/2 + 1 - y;
    }
    r2s = new int[numSamples];
    multipliers = new int[numSamples];
    target = new double[numSamples];
    int k = 0;
    for (int y = 0; y < (N-1)/2 + 1; y++) {
      for (int x = y; x < (N-1)/2 + 1; x++) {
        r2s[k] = x*x + y*y;
        target[k] = omega_c*j1(omega_c*sqrt(x*x + y*y))/(2*M_PI*sqrt(x*x + y*y));
        double window = cos(M_PI/2*sqrt(pow(x/((N - 1)/2 + 1.0), 2) + pow(y/((N - 1)/2 + 1.0), 2)));
        if (window < 0) {
          target[k] = 0;
        } else {
          target[k] *= window;
        }
        multipliers[k] = ((x > 0) ? 2 : 1) * ((y > 0) ? 2 : 1) * ((x > y) ? 2 : 1);
        k++;
      }
    }
    target[0] = omega_c*omega_c/(4*M_PI);
    minimum = new double[4*numComponents];
    maximum = new double[4*numComponents];
    k = 0;
    for (int i = 0; i < 4*numComponents; i++) {
      minimum[i] = optimalKernels[numComponents - 1][k]-pow(0.1, numComponents - 1);
      maximum[i] = optimalKernels[numComponents - 1][k]+pow(0.1, numComponents - 1);
      k++;
      if (k >= (numComponents-1)*4) {
        k -= 4;
        if (k < 0) {
          k = 0;
        }
      }
    }
  }

  ~CircularLPProblem() {
    delete[] minimum;
    delete[] maximum;
    delete[] r2s;
    delete[] multipliers;
    delete[] target;
  }
};

int main()
{
  INITKEYBOARD;

  CircularLPProblem problem(1, 41, M_PI); // Parameterize this!

  Opti::Strategy *optimizer;

  //optimizer=new Opti::G3(&problem, 2000);
  optimizer=new Opti::DE(&problem, 2000);

  printf("\nOptimizing...\n\n");

  for(int t = 0;; t++) {
    double bestcost = optimizer->evolve();
    if (!(t % 1000)) {
      printf("gen=%d, bestcost=%.20f, average=%.20f\n", t, bestcost, optimizer->averageCost());
      if (kbhit()) {
        printf("Parameter vector printout:\n");
        problem.print(optimizer->best());
        if (getch() == 27) break;
        getch();
      }
    }
  }
  delete optimizer;

  DEINITKEYBOARD;
  return 0;
}

Figure 2. Sum-of-squares cost of the best found approximate kernel $\tilde h$ as function of approximation order $M$, for $\omega_c = \pi$, $\omega_c = \pi/2$, and $\omega_c = \pi/4$, and $N=41$. A good indicator that these are the globally optimal parameterizations is the steady decrease of error as $M$ is incremented, except for $\omega_c = \pi/4$, $M=3$ for which we might (or not) have a sub-optimal solution.

This Python script generates the "exact" (windowed) kernel for comparison and implements the approximations:

import matplotlib.pyplot as plt
from scipy import special
import numpy as np
import scipy.ndimage
import skimage
import time

def circularLowpassKernel(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.fromfunction(lambda x, y: omega_c*special.j1(omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  kernel[(N - 1)//2, (N - 1)//2] = omega_c**2/(4*np.pi)
  return kernel

def rotatedCosineWindow(N):  # N = horizontal size of the targeted kernel, also its vertical size, must be odd.
  return np.fromfunction(lambda y, x: np.maximum(np.cos(np.pi/2*np.sqrt(((x - (N - 1)/2)/((N - 1)/2 + 1))**2 + ((y - (N - 1)/2)/((N - 1)/2 + 1))**2)), 0), [N, N])

N = 41  # Kernel width
M = 1  # Number of components
omega_c = np.pi  # Cutoff frequency <= np.pi
kernel = circularLowpassKernel(omega_c, N)*rotatedCosineWindow(N)

def saveKernel(name, kernel):
  plt.imsave(name+'.png', plt.cm.bwr(plt.Normalize(vmin=-kernel.max(), vmax=kernel.max())(skimage.transform.rescale(kernel, 4, 0))))
  absF = np.abs(np.fft.fftshift(np.fft.fft2(kernel)))
  plt.imsave(name+'_dft.png', plt.cm.Greys(plt.Normalize(vmin=0, vmax=absF.max())(skimage.transform.rescale(absF, 4, 0))))

saveKernel('exactpi', kernel)

plt.imsave('bwrkey.png', plt.cm.bwr(np.repeat([(np.arange(41*4)/(41*4-1))], 16, 0)))
plt.imsave('Greyskey.png', plt.cm.Greys(np.repeat([(np.arange(41*4)/(41*4-1))], 16, 0)))

def reComponentKernel(N, param):
  x = np.arange(N)-(N - 1)//2
  return np.exp(param[2]*x**2)*(param[0]*np.cos(param[3]*x**2) - param[1]*np.sin(param[3]*x**2))

def imComponentKernel(N, param):
  x = np.arange(N)-(N - 1)//2
  return np.exp(param[2]*x**2)*(param[1]*np.cos(param[3]*x**2) + param[0]*np.sin(param[3]*x**2))

optimal1pi = [0.90954095944946312,0.21237789645057972,-0.72663606661376057,0.76147736559107793]
optimal2pi = [-0.14935266179235707,0.28729650415029556,-0.08227954155725942,0.36475170816661134,1.03225137134955114,0.47158356759095016,-0.60039465413238513,0.47344404338750434]
optimal3pi = [-0.00646998371362690,-0.15983965263134517,-0.02321635125852370,0.20289085505437962,-0.30298121364781033,0.34395623806288950,-0.10323939264492392,0.32564190139614663,1.13148495464950427,0.66639187966234981,-0.54567361845465989,0.37519342198291905]
optimal4pi = [0.04481214216545543,0.08998210660870602,-0.01130470880820374,0.14003687599487286,0.02506113778123746,-0.24720956687380186,-0.03419309947921907,0.19643391040650712,1.27604032320787875,0.91159589755622883,-0.49585695598458995,0.29480666144138823,0.48259638667856242,-0.42280882153371496,-0.12730997070336811,0.29988730534029784]
optimal5pi = [-0.06413980110992069,-0.03272558679644168,-0.00736621171073370,0.10943881870260203,0.08213936888117918,0.16226536153011967,-0.01975126456281387,0.13958603391531316,-0.07270617897425770,0.37473199916953354,-0.04702984945995840,0.18925309227383197,1.43523098843984531,1.17056740170289952,-0.45755771004055446,0.24332707322808175,0.69700596730108921,-0.48854609666142051,-0.14144974773647198,0.27672159791886242]

class SeparableCircularLowpassFilter:
  def __init__(self, N, coefs):
    self.N = N
    self.reKernels = []
    self.imKernels = []
    for i in range(len(coefs)//4):
      self.reKernels.append(np.array([reComponentKernel(N, coefs[i*4:])]))
      self.imKernels.append(np.array([imComponentKernel(N, coefs[i*4:])]))
  def filter(self, x):
    reZ = scipy.ndimage.convolve(scipy.ndimage.convolve(x, self.reKernels[0].transpose()), self.reKernels[0]) - scipy.ndimage.convolve(scipy.ndimage.convolve(x, self.imKernels[0].transpose()), self.imKernels[0])
    for i in range(1, len(self.reKernels)):
      reZ += scipy.ndimage.convolve(scipy.ndimage.convolve(x, self.reKernels[i].transpose()), self.reKernels[i]) - scipy.ndimage.convolve(scipy.ndimage.convolve(x, self.imKernels[i].transpose()), self.imKernels[i])
    return reZ

filter1pi = SeparableCircularLowpassFilter(N, optimal1pi) # 1 component kernel
filter2pi = SeparableCircularLowpassFilter(N, optimal2pi) # 2 component kernel
filter3pi = SeparableCircularLowpassFilter(N, optimal3pi) # 3 component kernel
filter4pi = SeparableCircularLowpassFilter(N, optimal4pi) # 4 component kernel
filter5pi = SeparableCircularLowpassFilter(N, optimal5pi) # 5 component kernel

x = np.zeros([N, N])  # Input image
x[N//2, N//2] = 1     # (unit impulse)
# x = plt.imread('sample.tif').astype(float)[:,:,1]  # (green channel of some image from file)

t0 = time.time()
Z = scipy.ndimage.convolve(x, kernel)  # Exact
t1 = time.time()
print(t1-t0)

t0 = time.time()
reZ1pi = filter1pi.filter(x)  # 1 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ1pi', reZ1pi)

t0 = time.time()
reZ2pi = filter2pi.filter(x)  # 2 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ2pi', reZ2pi)

t0 = time.time()
reZ3pi = filter3pi.filter(x)  # 3 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ3pi', reZ3pi)

t0 = time.time()
reZ4pi = filter4pi.filter(x)  # 4 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ4pi', reZ4pi)

t0 = time.time()
reZ5pi = filter5pi.filter(x)  # 5 component kernel
t1 = time.time()
print(t1-t0)
saveKernel('reZ5pi', reZ5pi)

$M=1$
$M=2$
$M=3$
$M=4$
$M=5$

exactFigure 3. Kernels and the absolute value of their discrete Fourier transform (DFT), enlarged by a factor of 4 to make individual pixels visible. Top to bottom: Optimal approximating kernels $\tilde h$ with $M = 1$, $2$, $3$, $4$ and $5$ complex separable components, ideal kernel $h$ for $\omega_c = \pi$ and $N = 41$. Color keys: kernel: blue negative, white zero, red positive (normalized); abs DFT: white zero, black maximum.

Optimized approximation, $\omega_c = \pi/2$

Approximation of lower-cutoff kernels works better with this approach. I think this is because the number of ripples that fit in the windowed target kernel are reduced. For a high-cutoff filter then perhaps the kernel width $N$ could be reduced, as this would give the same number of ripples and frequency-domain approximation accuracy. But that would then favor a direct 2-d implementation (see Fig. 6). The C++ source code in the above was modified (not shown) for $\omega_c = \pi/2$ (and later for $\omega_c = \pi/4$) and the obtained parameters were used in a Python implementations of the approximate isotropic low-pass filters: (continued from the previous script)

omega_c = np.pi/2
kernelpi2 = circularLowpassKernel(omega_c, N)*rotatedCosineWindow(N)
saveKernel('exactpi2', kernelpi2)

optimal1pi2 = [0.44103810622146067,0.08998875769710178,-0.17010258583392401,0.19960767673288432]
optimal2pi2 = [-0.07233719880423649,0.14289630144713414,-0.02205699413927855,0.09165233018125875,0.51013245392078410,0.22632987351129516,-0.15215724343836151,0.12131467270512424]
optimal3pi2 = [-0.03867959516035375,-0.07754823299868645,-0.00905504286350691,0.05573851697983074,0.60460672854618647,0.36747347219951876,-0.13550005930232881,0.08061493799161984,-0.15528642640407436,0.23307152837452039,-0.03561076440525033,0.08259420496099962]
optimal4pi2 = [0.06323093028956613,0.02993129946018375,-0.00937470404526276,0.03826932634049100,0.78321346985185014,0.58823521191007977,-0.11611698422922974,0.05181454568143690,-0.07541120991444331,-0.16653642294319324,-0.01587217217852400,0.05400812761661748,-0.30454011473801174,0.38069299325271988,-0.04600795167685135,0.07149320869042795]

filter1pi2 = SeparableCircularLowpassFilter(N, optimal1pi2) # 1 component kernel
filter2pi2 = SeparableCircularLowpassFilter(N, optimal2pi2) # 2 component kernel
filter3pi2 = SeparableCircularLowpassFilter(N, optimal3pi2) # 3 component kernel
filter4pi2 = SeparableCircularLowpassFilter(N, optimal4pi2) # 4 component kernel

reZ1pi2 = filter1pi2.filter(x)
saveKernel('reZ1pi2', reZ1pi2)
reZ2pi2 = filter2pi2.filter(x)
saveKernel('reZ2pi2', reZ2pi2)
reZ3pi2 = filter3pi2.filter(x)
saveKernel('reZ3pi2', reZ3pi2)
reZ4pi2 = filter4pi2.filter(x)
saveKernel('reZ4pi2', reZ4pi2)

$M=1$
$M=2$
$M=3$
$M=4$

exactFigure 4. Kernels and the absolute value of their discrete Fourier transform (DFT), enlarged by a factor of 4 to make individual pixels visible. Top to bottom: Optimal approximating kernels $\tilde h$ with $M = 1$, $2$, $3$, and $4$ complex separable components, ideal kernel $h$ for $\omega_c = \pi/2$ and $N = 41$. Color keys: kernel: blue negative, white zero, red positive (normalized); abs DFT: white zero, black maximum.

Optimized approximation, $\omega_c = \pi/4$

In Python: (continued from the first Python script)

omega_c = np.pi/4
kernelpi4 = circularLowpassKernel(omega_c, N)*rotatedCosineWindow(N)
saveKernel('exactpi4', kernelpi4)

optimal1pi4 = [0.22210250507512549,0.04738076096442791,-0.04565916667208288,0.04924268271256500]
optimal2pi4 = [0.28445762468654334,0.14050386040493884,-0.03988634664177378,0.02123143861749057,-0.01717623293554685,0.11158938041615812,-0.01191733840600387,0.02614013272136841]
optimal3pi4 = [-0.00332428227606697,0.01722626235358831,-0.00516799783527921,0.00727825629264402,0.30296628929626013,0.16328004086794692,-0.03813290200362767,0.01805103148816284,-0.02296006483369628,0.12768470633404483,-0.01299438181561595,0.02558273151301529]

filter1pi4 = SeparableCircularLowpassFilter(N, optimal1pi4) # 1 component kernel
filter2pi4 = SeparableCircularLowpassFilter(N, optimal2pi4) # 2 component kernel
filter3pi4 = SeparableCircularLowpassFilter(N, optimal3pi4) # 3 component kernel

reZ1pi4 = filter1pi4.filter(x)
saveKernel('reZ1pi4', reZ1pi4)
reZ2pi4 = filter2pi4.filter(x)
saveKernel('reZ2pi4', reZ2pi4)
reZ3pi4 = filter3pi4.filter(x)
saveKernel('reZ3pi4', reZ3pi4)

$M=1$
$M=2$
$M=3$

exactFigure 5. Kernels and the absolute value of their discrete Fourier transform (DFT), enlarged by a factor of 4 to make individual pixels visible. Top to bottom: Optimal approximating kernels $\tilde h$ with $M = 1$, $2$ and $3$ complex separable components, ideal kernel $h$ for $\omega_c = \pi/4$ and $N = 41$. Color keys: kernel: blue negative, white zero, red positive (normalized); abs DFT: white zero, black maximum.

Computational complexity

Complexity as the number of 1-d dot products of length $N$ per pixel is $4M$ for the approximations and $N$ for the naive 2-d convolution implementation. The approximations suffer from an additional performance overhead due to the use of additional memory for storing intermediate results.

Figure 6. Complexity as the number of real multiplications per pixel vs. kernel width $N$. In reality the approximations are not quite as efficient, see Fig. 7. Possible advantages from symmetry of the kernels has not been taken into account.

Figure 7. Actual measured run time in seconds for 2-d convolution of a 5424×3636-pixel single-channel image by the ideal kernel $h$ and by separable approximating kernels $\tilde h$ with $M = 1$, $2$, $3$, $4$ and $5$ complex separable components, for $N = 41$, using Python's scipy.ndimage.convolve.

Further ideas

In order to adjust $\omega_c$ without kernel parameter re-optimization, a nearly continuous approximating kernel could be optimized by using a small $\omega_c$ and a large $N$. The result could then be sampled using interpolation.

It would be interesting to try direct optimization of the complex kernels without a parameterization that enforces circular symmetry. Perhaps the approximation could even be optimized as a sum of $M$ real separable kernels, similar to this Mathematics Stack Exchange question. Also minimization of (least squares) error in the frequency domain is something to try.


Why is bandwidth always limited in a real (physical) channel? I'm talking about a continuous analog channel.

Why can't it support infinite bandwidth?
Is there a physics reasons for it say for electrical signals?

Information bandwidth is dependent on signal to noise ratios.  At absolute zero, quantum level signal quantization and quantum noise will limit the lower bound on the noise floor.  At higher temperatures, thermal noise creates a higher noise floor in any information receiving equipment.  There may also be an upper limit on power density in a signal before its mass equivalent creates a black hole, which may be problematical in terms of information conservation.  In terms of transmitters, at some point you will melt the wires.

But in terms of pure math, there is no (aleph zero) finite limit on the amount of information that a single real number (one sample) can contain.


When sampling theory condition fails, what should I do? I met with a problem: if I sample audio data at 11kHz (yes, that is slow!) for 10 seconds and apply a Fourier Transform to the data, at what frequency will a 15 kHz tone appear?

So the original signal is 15kHz? which is greater than half of sampling rate 5.5kHz, right? and thus alias appears. And then, in this case, how should I find the appears one?

Thank you!

Because you've sampled at 11kHz the Frequency spectrum becomes periodic with a period of 11 KHz. So you're 15kHz will not only appear at 15 kHz but also at every frequency that is a multiple of 11 kHz away from 15 kHz. So the frequencies at which the 15 KHz tone are given by $ f= 15\times10^3 +N \cdot 11 \times10^3 $, where $N$ is an integer.

By convention you are looking at the frequency bins from -5.5 kHz to 5.5 kHz. Therefore you're 15 kHz tone will appear at 4 KHz.

Note that if your original signal was real then there is also a freqency component at -15kHz. This also gets extended periodically every 11 kHz - so this will produce a corresponding peak at -4 kHz.