Perform Convolution in Frequency Domain Using FFTW? I'm trying to convolve two signals $x(n)$ and $h(n)$ in C by using the FFTW library's functions to perform a Fourier transform on each, multiply the appropriate complex components together, and take the IFFT of the resultant product. Unfortunately, my IFFT'd result isn't correct, and I'm not sure what I'm doing incorrectly.

I'm testing using a modified version of code I found online (.c and Makefile here), which populates two arrays with random numbers between 0 and 1. My process is as follows:

Create double buffers for x (audio signal) and h (impulse response) of lengths n and nIR respectively.
Populate each buffer with random numbers between 0 and 1.
Create two complex buffers X and H, each of length nOut = n + nIR - 1
Create and then execute two real-to-complex FFT plans, one for each of the signals, both using nOut elements (i.e. plan_forward = fftw_plan_dft_r2c_1d ( nOut, x, X, FFTW_ESTIMATE ); and plan_forwardIR = fftw_plan_dft_r2c_1d (nOut, h, H, FFTW_ESTIMATE);
Create complex array fftMulti to hold the frequency products of the two arrays. Its size is nc = nOut / 2 + 1 because FFTW doesn't store the redundant half of the FFT
Loop for ( i = 0; i < nc; i++ ) and perform complex multiplication, i.e. fftMulti[i][0] = X[i][0] * H[i][0] - X[i][1] * H[i][1]; for the real component and fftMulti[i][1] = X[i][0] * H[i][1] + X[i][1] * H[i][0]; for the imaginary component
Create double array convolvedSig of size nOut. This array will hold the IFFT of fftMulti
Create and execute a complex-to-real FFT plan of size nOut, i.e. plan_backwardConv = fftw_plan_dft_c2r_1d(nOut, fftMulti, convolvedSig, FFTW_ESTIMATE);
Normalize convolvedSig by dividing each element by nOut

The values in convolvedSig are not correct (they are way too high), but I'm not sure what I'm doing wrong. I also created arrays and complex-to-real plans for my unmodified signals X and H (also of size nOut), and those IFFTs worked just fine (i.e. values were exactly the same as before I performed FFTs on them).

Using my process and/or my code, can someone please help me identify what I'm doing incorrectly?

You need to zero pad x and h.

nZeros = length(x) + length(h) -1

Example in MATLAB:

clear all; close all; home

xSamples = 10;
hSamples = 7;

x = randn(1,xSamples);

h = randn(1,hSamples);

nZeros = xSamples + hSamples - 1;

X = fft(x,nZeros);
H = fft(h,nZeros);

x_conv_h = ifft(X.*H);

figure(1)
plot(real(x_conv_h))
hold on; grid on
plot(conv(x,h),'.r')
legend('fft based convolution','convolution')


How to show that y[n] = x[n] * h[n] turns into the Y(z) = X(z).H(z)? I'm trying to show that $y[n]=x[n]*h[n]$  turns into $Y(z) = X(z)H(z)$ in Z-domain by first applying convolution then by taking the inverse Z-transform of the $Y(z)$, stating that it's the same sequence after all. However, in power series expansion an extra term appears which kills the equation. You can find the question below; I added it as an image.

I'd be grateful if you could help.
Thank you

The correct sequence $y[n]$ is the one you obtained via the $\mathcal{Z}$-transform:

$$y[n]=\frac{1}{a-b}\left[a^{n+1}-b^{n+1}\right]u[n]\tag{1}$$

I'm not sure where that $n+1$ term is coming from in your time-domain derivation, so I'll write out the convolution sum:

$$y[n]=\sum_{k=-\infty}^{\infty}a^ku[k]b^{n-k}u[n-k]=u[n]\sum_{k=0}^{n}a^kb^{n-k}=u[n]b^n\sum_{k=0}^{n}\left(\frac{a}{b}\right)^k=\\=u[n]b^n\frac{1-\left(\frac{a}{b}\right)^{n+1}}{1-\frac{a}{b}}=u[n]b^{n+1}\frac{1-\left(\frac{a}{b}\right)^{n+1}}{b-a}=u[n]\frac{b^{n+1}-a^{n+1}}{b-a}$$

which equals (1), as it should.


How to transform a signal to go through specific points? I have a 1d signal obtained using a Fourier based resample method (TDIFDZP) for which the resampled points don't necessarily go through the original samples. I want to transform the upsampled signal to go through these points to reduce the error between the resampled and the original signal as shown in the image below. I'm using scipy.sigal.resample implementation in python.
I'm looking for a transformation of the resampled signal so that it goes through the original samples.

Edit: A little of background, I'm exploring different interpolation methods for sub-Nyquist sampling. My original signal's Nyquist rate is ~20Hz but the number of samples available in my application is <6.

To have all the samples in the interpolated waveform land on the original samples implies using an integer interpolation rate (otherwise time domain distortion would need to be introduced, and I can't think of an application where that would make sense to do).
For integer interpolation, a very simple approach is to use the DFT with zero padding since the samples chosen will be exact, but this will have more deviation in between then other optimized interpolation methods further detailed below. This approach is very simple, so if the deviation between samples is of no concern then this approach would be favored and would be as follows:
Compute the FFT of the samples without zero padding.
Scale the result by $N$, where $N$ is the total number of original samples, zero pad the result and then compute the inverse FFT.  (Due to reciprocity the fft and ifft operations can be swapped).
In MATLAB or Octave this is:
N = length(test);
interp = 5;
result = ifft(N*fft(sequence), interp*N);

Other optimized interpolation methods using zero-insert and subsequent image filtering would also preserve the values, and with post processing the image filtering can be done with a zero-phase filter using filtfilt() and for an optimized solution in the least squares sense, together with multiband filter coefficients using firls(). The zero-insert and filtering approach is further detailed here: What is the impulse response used in an interpolation filter when upsampling?
If your samples are non-uniform to begin with, first resample to a uniform rate, where this post already addresses that approach: What is an algorithm to re-sample from a variable rate to a fixed rate?
But if the new samples do not fall on the same time locations (which will occur if there isn't an integer relationship in the duration between the original samples) then the interpolated samples should not be expected to match the original samples when they occur at new time samples that were not in the original sequence.


Looking for a analytical formula to compute the central frequency of a signal analyzed by a discrete wavelet at a given scale I am reading this paper by Han et al. (2014). In this article, the authors extract detailed information from geomagnetic data sampled at 1Hz using a daubechies 5 (db5) wavelet: they reconstruct the signal using the coefficients from the 6th level only. The authors explain that the central frequency of the resulting reconstructed signal is 0.01Hz.

I have reprocessed similar data following the same steps and have convinced myself that the central frequency of the signal reconstructed using only coefficients from the 6th level is indeed around 0.01Hz. Below are the PSD plots for both the original and reconstructed signals, generated using Python:

But is there a way to analytically determine the central frequency of the reconstructed signal?

On one hand, I have found explanations stating that, considering a SR of 1Hz, the central frequency of the signals reconstructed using the 6th level should lie within [0.5, 1]/2**6, i.e. [0.0078125, 0.015625]Hz, which is the case. Is this the correct answer? If this was the case, it would mean that the central frequency of the analyzed signal would be "wavelet-independent"...
On the other hand, I know about the formula from Abry, 1997, stating that the pseudo-frequency of a wavelet at a given scale Fa is the central frequency of the wavelet Fc divided by the scale a: Fa=Fc/a. But my undestanding is that this formula yields the pseudo frequency of the wavelet, not the central frequency of the analyzed signal!

Any clarification would be appreciated. Thanks!

They do mean the pseudo-frequency of the wavelet which is not dependent on the signal being analyzed. The misleading terminology that they use seems to come from from one of the references, Han, P. (2013), Investigation of ULF seismo-magnetic phenomena in Kanto, Japan during 2000–2010, PhD thesis, Chiba University, Chiba, Japan. Quoting an article of a similar name and where Han is one of the authors, Hattori et al., Investigation of ULF Seismo-Magnetic Phenomena in Kanto, Japan During 2000–2010: Case Studies
and Statistical Studies, Surv Geophys (2013) 34:293–316,
DOI 10.1007/s10712-012-9215-x, directly answers your question:

  In this paper, we yield discrete wavelet transform (DWT) and decompose the original
  signature into six levels. The central frequency ($F_a$) of signals in each level is given as

  $$F_a = \frac{F_c}{a\Delta} \tag{4}$$

  where $F_c$ is the frequency of basic wavelet; $a$ is the scale of the
  level, and $\Delta$ is the sampling rate of original data. After evaluating
  different kinds of wavelets, we adopted Daubechies 5 (db5) as the
  mother wavelet. Moreover, the db5 wavelet has been proved to be
  effective in the previous study (Jach et al. 2006). According to Eq.
  4, for db5 wavelet, the central frequency of signals in the sixth
  level is 0.01 Hz.


Algo Name: Detecting Elongated Shape in 2D Image by Convolution with Directed Lines I am detecting a shape with roughly known orientation (obtained by using Hough line transform) in thresholded binary image by convolving the image with line with that orientation, and counting the ratio of pixels which match; this line is swept across the image, and lines which have hit ratio higher than some threshold (like 0.9) are considered belonging to the object completely.

The motivation is not only to remove holes but also discontinuities across the object; how a hole is removed is shown in the image, where grey line parts are not-matching and red are matching; blue lines are from Hough transform:

This routine was found by intuition and I wonder whether there is some name for it, or perhaps a more sophisticated algorithm based on similar idea. In that case I could use an optimized implementation instead of the hand-written code (I am using python, OpenCV and scikit-image).

What you are essentially doing is a matched filter. However, thanks to Hough transform, your filter (line) is oriented and therefore I would call it an oriented matched filter. For generating the Bresenham line and sampling the pixels you might want the use the OpenCV line iterator. The simple usage would be similar to:

cv::LineIterator it(image, pt1, pt2, 8);
std::vector<cv::Vec3b> buf(it.count);
std::vector<cv::Point> points(it.count);

for(int i = 0; i < it.count; i++, ++it)
{
    buf[i] = *(const cv::Vec3b)*it;
    points[i] = it.pos();
}

The hole detection part can also be implemented using convolutions as follows:

Convolve the image with a horizontal convolution kernel composed of the oriented line. Use filter2D.
Threhsold the result to identify the areas with holes (use threshold).
Perform a blob-filtering and select the desired hole.

Even though I think that this approach is really valid and would work fine in many cases, the typical blob-person would approach the problem differently: First identify the large connected components in the image and filter out the rest. Within that single large component, perform a morphological hole filling. This is a pretty standard pipeline, and with that method, you will be able to find more tools (blob analysis) to help you. For example OpenCV methods such as findContours, drawContours or morphologyEx come handy.

Finally, if you know the shape of the hole, you might benefit from e.g. ellipse detection to directly find it.

Hope these help.


Matlab HIgh frequency sine wave generation I need to generated a high frequency sinusoidal signal for modulation in MATLAB. But it turns out to be something weird. This is the simple code snippet I used...

t = 0:0.001:100;
A = 1;
s = A*sin(2*pi*1e9.*t);  %1 GHz signal
plot(t, s);
title('Modulating Signal');
xlabel('Time');
ylabel('Amplitude');

Instead of the sine wave I am getting something weird like this..

The function you need to evaluate is:

A*sin(2*pi*FREQUENCY*(1:NUM_SAMPLES)/SAMPLE_RATE)

Use a reasonable sample rate (at least twice the frequency) and a reasonable interval (remember that at a sample rate of 2Ghz, 1 second of signal will fill up the memory of your computer) ; and it'll work! From your code it seems that you are trying to generating 100s of signal at a sample rate of 1kHz, which doesn't make sense.

At the moment, in your code, the argument of the sin function is an integer multiple of $2\pi$ so in theory your code should be plotting a flat line of 0s. You are seeing something a bit different due to the finite resolution of floating point numbers.


Is it possible to encrypt an audio wave? NFL Example In the NFL, quarterbacks receive playcalls from coaches over their helmet's wireless interface.

Someone mentioned how occasionally, these radios will pick up interference from airplanes, and that frequently the helmets lose communication.

One thing I wonder about is if an opposing team's coaches could pick up what was being communicated between the offensive coordinator and the quarterback.

And then I wonder, well, why wouldn't there just be an encryption scheme between the coach and the quarterback.

A chip on the coach's side would take his voice, encrypt it, send it to the quarterback's helmet, in which a matching chip would decrypt it.

Obviously there are tons of other applications, military, etc.

But then I wonder if the encrypted signal would be distorted with all this noise, and if that distortion would completely disrupt the decrypted signal.

Because if you change even a bit in an encrypted signal, it would obfuscate the resulting signal when it was decrypted by the quarterback's helmet.

Are there any methods for this kind of "loose encryption"?

There are "human" encryption schemes. Like NFL coaches could change the name of their plays from game to game. This would keep it known to only the QB and the coach.

But I guess I was wondering if there was a mathematical approach.

Someone mentioned how occasionally, these radios will pick up interference from airplanes, and that frequently the helmets lose communication.

Is this the 1970s?

Modern headsets should be pretty robust against interference, and they have a dedicated spectrum allocation. If there's interference from airplanes, either device malfunctions.

So, no, this sounds like rumors, or long obsolete claims, especially since for superposition reasons and 6 decades of compatibility, aircraft voice comms is partly still AM; modern (as in: post-1970) headsets are either digital or at the very least FM.

  One thing I wonder about is if an opposing team's coaches could pick up what was being communicated between the offensive coordinator and the quarterback.

Sure. It's on the air. Radio waves can't tell friend from foe. And since a quarterback can be all over the field, there's certainly no high-gain antenna focusing the radio beam on the player alone.

  And then I wonder, well, why wouldn't there just be an encryption scheme between the coach and the quarterback.

Just because you can receive the mobile phone call of your neighbor, or his WiFi, or pay TV, doesn't mean you can decrypt stuff.

So, yes, there's plenty encryption systems in the technological fields of voice comms, some civil, some military... plenty to pick from.

  But then I wonder if the encrypted signal would be distorted with all this noise, and if that distortion would completely disrupt the decrypted signal.

Digitally, whether or not the signal is encrypted makes little difference to the signal quality in case of errors. You typically break one full "frame" worth of channel coding / encryption.

All wireless digital communication schemes incorporate channel coding that allows for forward error correction (FEC) so that this makes no difference, and the audio codecs (specifically: the VoCoder) for such communication channels are robust against "surviving" errors.

  Are there any methods for this kind of "loose encryption"?

I'll not go into crypto here, because a) that's not my forte and b) because this isn't the perfect platform, but:

No, that's something a cryptographer must never implement; every encryption must have zero correlation between encrypted and decrypted bits. If that wasn't the case, the cryptography would be inherently insecure, because you could, through statistics, relatively easily reconstruct (parts of) the unencrypted communication.

There are "human" encryption schemes. Like NFL coaches could change the name of their plays from game to game. This would keep it known to only the QB and the coach.

  There are "human" encryption schemes. Like NFL coaches could change the name of their plays from game to game. This would keep it known to only the QB and the coach.

That would be called a "substitution cipher", and these are extremely weak. Especially, you call a specific move "rabbit", and then the quarterback does a move, for the rest of the game it's clear what "rabbit" means.

Commentary on the situation: Obviously, the solution to this problem is organizational, not technological: It's a sport. If you cheat, you get disqualified. So, there's strong incentive not to take performance-enhancing drugs. And there's just as much incentive to not use information you mustn't get.

Just that the supervision of the former requires invasive, expensive testing, and the question what kind of knowledge a coach had can be relatively easy tested within similar margins of error as drug tests with basic statistics.


Implementing FIR Hilbert Transformer in Python I want to implement FIR Hilbert Transformer in Python using window method. I'm new to digital signal processing so please point out to me all mistakes you notice.
Ideal Hilbert transformer is described as:
$H(e^{j\omega}) =
    \begin{cases}
        -j & \text{, $0 < \omega < \pi$} \\
        j & \text{, $-\pi < \omega < 0$}    \\
    \end{cases}
$
from which we conclude that ideal Hilbert Transformer is all pass filter $|H(e^{j\omega})| = 1$.
I know that FIR filter can be implemented using scipy.signal.firwin function.
What I don't know is which parameters to pass to function in order to get Hilbert filter aproximation.

You can't implement an ideal Hilbert transformer: its impulse response is non causal and infinite in time.
So you can only implement an "approximate" one and the best way to do this depends on the specific requirement of your application: what's your frequency range of interest, how much magnitude and/or phase deviation can you tolerate, are you constrained in latency or computational resources, do you need to maintain causality and transients, etc.
Python has a Hilbert transform function https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hilbert.html but without understanding your requirements first, it will be difficult to control the parameters properly.
Very good article on Hilbert transform and analytical signals in the discrete domain: http://andrewduncan.net/air/


Separate Signal Values from Noise let us consider following graph of singular values

i want to make  some kind of  clustering of these data,namely to  seperate  main components from non main components,let say signal components from noise components,i would like to do like this(any software matlab,etc is great,better matlab)let us start like this  ,first do linear regression  with first point,then do regression analysis with first two point and so on,point is that  one want to  create two group,in first group   coefficients of  regression lines should be  close to each other, also in another group  regression coefficients should be close to each other,but regression  coefficients  in one group must be different then coefficients in another group,that means that we should find such point which separate this group optimally,  in other word i should stop regression analysis at some point which  well separate  two group,and another regression will start from this stop point till other rest points,please help me how to do it programaticaly

I would use the SVD (Singular Value Decomposition).
By looking at the Singular Values I'd determine which vectors spread the data and which spread the noise.

Practically, they both do both, but if we speak which are dominant, this would be a great starting point.


How calculate channel gains in NOMA? I am developing a user group algorithm in NOMA, but I need channel gain for that.
I know that the channel gain consists of some components, but I have not found a way to calculate the gain.
Does anybody know?

You won't know the exact channel, instead you'd use the channel estimates.
The basic idea behind Non Orthogonal Multiple Access (NOMA) is that you are receiving multiple overlapping signals, decoding the strongest signal, subtracting it, decode the next strongest signal, and so on. With that being said, NOMA and the associated successive interference cancellation (SIC) receiver requires a high enough SNR so as to be able to reliably decode all of the signals. For this high SNR case, it is reasonable to assume that the channel estimates you have will be good enough for the processing.


How does an equalizer prevent a sample and hold effect? How does an equalizer prevent a sample and hold effect?

Does a reconstruction filter only smooth out the output of an DAC (+S&H) or does it have another function?

Equalizers do not prevent the sample and hold effect, though if you are talking about compensation filters, see the edit below.

There are two types of reconstruction filters (also known as low-pass filters): analog and digital.  An analog low-pass filter can smooth out the output of a DAC.  A digital low-pass filter cannot, since it comes before the DAC.  What it can do is reduce the effect of the sample and hold by increasing the sample rate (assuming it is an interpolation filter) and then getting rid of the aliases via the low-pass filtering.  The effect of the sample and hold is reduced due to the increased sample rate.

EDIT: As Robert points out, the frequency droop that happens as a result of the zero-order hold can be partially compensated for with digital filters whose frequency domain gain is the inverse of the zero-order hold effect.  While this will not eliminate the sample and hold effect, it will mitigate it.


Is there a public repository of labeled sound files coming from an industrial area? I know about the UCI Machine Learning Repository for ML datasets; Where could one find sound files e.g. in wav format, which are partially or completely labelled?

Freesound is a repository of sound files categorised by user-defined tags. Through those, you can spot what other users have labeled industrial sounds but be prepared that some of those might be artificially created (so, not really a recording of an engine but something that sounds like an engine).

In the same category is also Soundsnap.

The other popular choice would of course be Soundcloud but there, you will have to do a bit of digging to find sounds that suit your use case because the majority of sound clips there are about music (so, "industrial rock" for example).

In all of these, you can curate a subset of samples that seem to have a strong relationship between their content and the tags, to then use the tag as your class variable if you are looking to download a large dataset.

Hope this helps.


Discrete time aperiodic signal: Concept of frequency? Is concept of frequency to a-periodic signal whether it is continuous or discrete time is valid?

I am bit confused, since discrete time sinusoids have frequency $-\pi<\omega<\pi$ or $-\frac{1}{2}<f<\frac{1}{2}$, how can I relate these to discrete time a-periodic signals?

Also how is this relation $f=\frac{F}{F_s}$ applied to a-periodic signals?

TL;DR

Your statements are inaccurate. The concept of frequency is $\frac{1}{T}$ where $T$ is the period time. Be careful with the continuous-time and discrete-time definition for a perioud.

  discrete time sinusoids have frequency −π<ω<π or $−\frac{1}{2}<f<\frac{1}{2}$

So, this statement is not accurate.

A frequency is defined as the number of times a signal repeats itself in a second and is measured in $[Hz]=[\frac{1}{sec}]$. This is, by definition, not limited to the interval you stated and, infect, not defined well for negative values at all.

An angular frequency, $\omega=2\pi f$ is the rate of change of angular displacement, θ, (during rotation), or the rate of change of the phase of a sinusoidal waveform. It is measured in angle per second or $\frac{radians}{sec}$. As can be seen from the units, this is also not limited to an interval.

When we start to analyze using the Fourier and Laplace transforms, we analyze in the frequency domain. First, the Laplace transform analyzes the signal and if we want to consider frequencies we convert it to continuous-time Fourier transform (CTFT) by using $H^F(\omega)=H^L(s=j\omega)$, which exist iff $0\in ROC$ for the Laplace transform. Otherwise, CTFT is not defined.

CTFT is a function of $\omega$ and is not limited to an interval. However, let us look at a $T$ periodical signal $x(t)=x(t+T)$ with a CTFT $X^F(\omega)$.
$$x(t)=x(t+T)$$
$$\mathcal{F}\{x(t)\}=\mathcal{F}\{x(t+T)\}$$
$$X^F(\omega)=X^F(\omega)e^{-j\omega T}$$
$$X^F(\omega)-X^F(\omega)e^{-j\omega T}=0$$
$$X^F(\omega)(1-e^{-j\omega T})=0$$

Either $X^F(\omega)=0$, or $(1-e^{-j\omega T})=0$ and sinc $(1-e^{-j\omega T})=0\implies\omega T=2\pi k$ we know that $X^F(\omega\neq\frac{2\pi k}{T})=0$. This means that the CTFT of a periodic signal is discrete.

  Is the concept of frequency to a-periodic signal whether it is continuous or discrete-time is valid?

You have to be careful here. Is the discrete signal $\sin[n]$ a periodic signal? The surprising answer here is no. A descrit-time periodic signal is defined as a signal for which exist $N\in \mathbb{N}$ S.T. $x[n]=x[n+N]\forall n$. Let us consider $0=\sin[0]=\sin[N]\implies N=2\pi k \notin \mathbb{N}$. on the other hand $\sin(t)$ is most definitely periodical. The concept of frequency in both cases is $\frac{1}{T}$ where $T$ is the time of a single period. For the continuous-time, it is the smallest $T$ for which $x(t)=x(t+T)$ and for a periodic discrete signal, the smallest $N$ for which $x[n]=x[n+N]$. The frequency of both is $\frac{1}{T},\frac{1}{N}$ respectively.

the discrete-time Fourier transform (DTFT) is the analogy of CTFT for sampled\descrete signals. The difference between both is that it is periodical with respect to the sampling rate. The time distance between every two samples is $T_s$ hance the sampling frequency is $F_s=\frac{1}{T_s}$ and $\omega_s=\frac{2\pi}{F_s}=2\pi T_s$. we now take the original Fourier transform and replicate it around $k\omega_s$ where $k\in\mathbb{N}$. After this operation, we normalize (divide) by $\omega_s$. Now, the DTFT is a function of $\theta$ which is limited to $-\pi \leq \theta\leq \pi$ and is periodic around $2\pi$.
All this can be written as: $$X^f(\theta)=\frac{1}{T_s}\sum_{k=\infty}^{\infty}X^F(\frac{\theta-2\pi k}{T_s})$$

Maybe this is what you meant when you wrote the limitation on $\omega$?

I am not sure where did you get the final relation. Maybe a reference will help and I will edit the answer.


How do you send signals through geological formations Does anyone have suggestions for performing the above?  Presently, in my job we send signals through geological formations; however, if the rock is too dense, then we are unable to receive any signals.

The signals must be sent from a battery-powered transmitter.  Here is an article I've found sending neutrinos, but I am more interested in sending electromagnetic radiation:

http://www.gizmag.com/neutrino-messaging/21842/

For digital communications, bit error rate performance is all about $\frac{E_b}{N_0}$, the energy per bit divided by the noise power spectral density. The part that you can control is the energy per bit. Energy per bit is equal to $\text{(average power)(bit duration)}$. You can increase $E_b$ (and therefore decrease your bit error rate) in two ways:

Raising your transmitter power. This may or may not be possible based on the hardware you're using, the frequency band in use, and any governmental regulations for your locale.
Decreasing the bit rate (i.e. give more time per bit). In cases of extremely high path loss, like communication into caves, you can get reliable operation at a very low bit rate (although other issues such as long-term oscillator stability can come into play and diminish the gain that you might expect).

Communication into deep geological formations is something of a niche application, so you may have trouble finding a lot of literature on suggested approaches. There are companies who specialize in making this type of equipment, so it may serve you better to consult the marketplace for existing solutions that might be off the shelf instead of trying to engineer one from scratch (if your budget permits, of course; I expect such systems to be expensive due to the rather small market).


DFT product of sinusoids (From Shaums DSP outline, 2nd edition, page 248, problem 6.21)

Book says, evaluate the Sum:

$$
S = \sum^{N-1}_{n=0} \Bigg( x_1[n] \ x^{*}_2[n] \Bigg)
$$

when:

$$
\begin{aligned}
x_1[n] = \cos\left( \frac{2\pi n k_1}{N} \right) \\
\\
x_2[n] = \cos\left( \frac{2\pi n k_2}{N} \right)
\end{aligned}
$$

using Property:

$$
\sum^{N-1}_{n=0} \Bigg( x_1[n]\  x^{*}_2[n] \Bigg) = \frac{1}{N} \sum^{N-1}_{k=0}  \Bigg(  X_1[k]\ X_2^{*}[k]\Bigg)
$$

I start by converting cosines into complex exponentials:

$$
x_1[n] = \cos\left( \frac{2\pi n k_1}{N} \right)
$$

$$
x_1[n] =  0.5e^{j2\pi n k_1 / N}  + 0.5e^{-j 2\pi n k_1 / N}
$$

Applying Definition of DFT to find:
$$
X_1[k] = \sum_{n=0}^{N-1} \Bigg( 0.5e^{j2\pi n k_1 / N}  + 0.5e^{-j 2\pi n k_1 / N} \Bigg) e^{-j2\pi nk/N}
$$

$$
X_1[k] = \sum_{n=0}^{N-1} \Bigg( 0.5e^{-j(2\pi n/N)(k- k_1)}  + 0.5e^{-j(2\pi n/N)(k+ k_1)} \Bigg)
$$

From this I determine that sum over one period of a complex exponential is zero, except when k is selected to cancel out exponential function input to zero, such as when $k=k_1$ and $k=-k_1$.  Thus:

$$
X_1[k] = \begin{cases}
\frac{N}{2}&(k=k_1)\ or\ (k=-k_1) \\
0 & else
\end{cases}
$$

By similar means:
$$
X_2[k] = \begin{cases}
\frac{N}{2}&(k=k_2)\ or\ (k=-k_2) \\
0 & else
\end{cases}
$$

Here's where I get lost:

Books says, therefore:

$$
\sum^{N-1}_{n=0} x_{1}[n] x^{*}_{2}[n] = \frac{1}{N}\Bigg[ \frac{N^2}{4} + \frac{N^2}{4} \Bigg] = \frac{N}{2}
$$

Where does:
$$
 \frac{N^2}{4} + \frac{N^2}{4}
 $$
come from?

Actually you have quite done it, but lets indicate.

When the frequencies of both sinusoidal signals $x_1[n]$ and $x_2[n]$ are the same; i.e., $k_1 = k_2 = m$, then the impulses in the corresponding N-point DFT sequences $X_1[k]$ and $X_2[k]$ will occur at the same index $k=m$ and $k=-m$ (or $N-m$) with a weight of $N/2$.

Hence multiplication of $X_1[k]$ with $X_2[k]$  will be $N/2 * N/2 + N/2 * N/2 = N^2/2$, as there are two impulses in the DFT of a cosine signal...

Finally dividing by $N$ yields the result as $N/2$ when the frequencies are the same and $0$ when not.


Bode plot of discrete-time transfer function $H(z)$ $H(z)$ is the transfer function of a biquad filter as described here.

I would like to plot the Bode plot of the magnitude response of $H(z)$.

Scipy has a bode method (scipy.signal.bode) for continuous-time transfer functions. Is there an option I didn't see for discrete-time functions?

I found in this nice file formula 18. I tried plotting it with gnuplot directly, but that didn't look correct.

set samples 100000, 100000
set logscale x

b0= 0.2514
b1= 0.5028
b2= 0.2514
a1=-0.1712
a2= 0.1768

H(x) = sqrt((b0**2 + b1**2 + b2**2 + 2*(b0*b1+b1*b2)*cos(x) + 2*b0*b2*cos(2*x))/(1 + a1**2 + a2**2 + 2*(a1+a1*a2)*cos(x) + 2*a2*cos(2*x)))

plot [1:22050] H(x)

So, how do I do this? Is  it even possible to Bode plot the magnitude response of a discrete-time transfer function?

You can use the freqz command, available in Matlab/Octave and SciPy. In SciPy it is available as scipy.signal.freqz.

For your case it would be scipy.signal.freqz([b0, b1, b2], [a1, a2])

Here is the reference for more info:

http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.freqz.html


FMCW Radar using FSK Superimposed LFM ramp While reading the datasheet of the ADF4159 I came upon a section describing how to create an FMCW ramp with FSK superimposed to "enable unambiguous (distance and velocity) multitarget detection". While this sounds like a wonderful radar technique, I do not understand how it works. There does not seem to be a lot of online documentation on it.
Datasheet: https://www.analog.com/media/en/technical-documentation/data-sheets/ADF4159.pdf
Page of the section: 29
I am especially confused with the picture they provide of the frequency domain:

What are the dashed lines supposed represent? Why does each step seem so "uneven" with each other? What exactly is the smaller dimension arrows/line measuring?
I already understand the concepts that drive FMCW radar using a linear sawtooth or triangular ramp, so could anybody explain how LFM-FSK radar works, possibly in a visual way? How exactly does it interact with the target? What does it look like on the frequency and/or time domain? And is this method of measurement worthwhile (in terms of complexity vs functionality) compared to other alternatives?
I simply seek knowledge here, I don't plan to try to detect missiles anytime soon :)

First thing first: Once you start talking about stepped FM waveforms with the goal of achieving LFM-like range resolutions you are introducing a more complex receiver along with other associated disadvantages. Stepped-LFM waveforms require that you process each pulse individually, and thus the receiver must reconfigure itself for every pulse which are at different frequencies. Having said that, you must transmit multiple pulses (think multiple PRIs) in order to reap the range resolution benefits. Contrast this to an LFM system that has the large enough instantaneous bandwidth where it can achieve high range resolution with only a single pulse.
This documentation implies that you use the stepped-FM portions of the waveform to achieve high range resolution and the FSK processing to perform Doppler filtering.
Specifically, with this waveform you can process the strictly stepped-FM portions (the solid black lines) as mentioned before to achieve the range resolution required. You can then pass the FSK portions (dotted lines) through a Doppler filter bank, after appropriate mixing and integration, to yield the target velocity.
As you can see, the receiver for such a system sounds pretty complex! In the past, hardware was not up to par to provide the instantaneous bandwidth and computational power that we have today, so great effort was made to design these complex receivers to achieve the desired performance.
Thankfully, hardware has come a long way even since this datasheet was published and we now have access to high instantaneous bandwidth systems which allow for simpler receiver designs. Modern DAC/ADC and embedded processor technology has enabled multi-GHz bandwith waveform generation  and processing. We can easily generate a multi-hundred-MHz LFM pulse and process it digitally via a matched filter using a rather simple receiver design.


How to Check Separability of 2D Filter / Signal / Matrix Given:
    x(n1,n2) = {1 ,n1=0,n2=0 ;
                2 ,n1=1,n2=0 ;
                3 ,n1=0,n2=1 ;
                6 ,n1=1,n2=1 }

How could one prove it is separable?

Nilesh Padhi, Welcome to the DSP Community.

The classic definition of separable means the data (2D) given by $ X \in \mathbb{R}^{m \times n} $ can be written as:

$$ X = \sigma u {v}^{T} $$

Where $ \sigma \in \mathbb{R} $, $ u \in \mathbb{R}^{m} $ and $ v \in \mathbb{R}^{n} $.
This is called Rank 1 Matrix.

How can you get those parameters and vectors given $ X $?
Well, the Singular Value Decomposition (SVD) is here to save the data.

The SVD of $ X $ is given by:

$$ X = U \Sigma {V}^{T} = \sum {\sigma}_{i} {u}_{i} {v}_{i}^{T} $$

You can see those match when $ {\sigma}_{j} = 0 $ for $ j \geq 2 $.
So what you should do is the following:

epsThr = 1e-7;
[mU, mD, mV] = svd(mX);
vD = diag(mD);
if(all(vD(2:end) < epsThr))
    vU = mU(:, 1);
    vV = mV(:, 1);
end

We checked if indeed the singular value of 2 and onward are small.
If they do (You can decide to what degree of small by epsThr) then it is separable and the vectors are vU and vV.

In your case:

mX = [1, 3; 2, 6];
[mU, mD, mV] = svd(mX);

vD = diag(mD);
disp(vD);

The result is:

vD =

    7.0711
    0.0000

Since vD values are zero besides the first element (Only single non vanishing Singular Value) it is separable.

Indeed you can see that:

mD(1) * mU(:, 1) * mV(:, 1).'

ans =

    1.0000    3.0000
    2.0000    6.0000

As expected.

This method is really useful in Image Processing when we want to convolve with 2D kernel and we find it is separable and hence we can apply the 2D convolution using 2 1D convolutions (Along Columns / Rows).

In that case we define $ \hat{u} = \sqrt{{\sigma}_{1}} u $ and $ \hat{v} = \sqrt{{\sigma}_{1}} v $ where $ {\sigma}_{1} $ is the Singular Value.
Then we convolve $ \hat{u} $ along columns and $ \hat{v}^{T} $ along rows.


How do I convert a timeseries to a different frequency band? I have a real time series sampled at 32 MHz. So, when I channelized it and plot it against time, I get a Frequency vs. Time image plot where the frequency domain spans 16 MHz. To elaborate, this basically entails taking, say, 2048-sample chunks of data and doing a Fourier Transform on each to get a 1024-channel '16 MHz'-wide Fourier Transform over time (because we only care about positive frequencies). It looks something like this:

Every vertical line in that plot is a Fourier Transform of 1024 channels. Which means that the time step between adjacent vertical lines is 2048 samples wide or 64 microseconds wide.

Now, what I want to do is convert this time series so that it looks like it was sampled at 48000 Hz, but not actually. Basically, I want the 16 Mhz frequency channel to correspond to a 24 kHz channel. So, I still want to preserve the structure of the data. So, a naive low-pass filter doesn't work, because I don't actually want to remove the higher frequencies, I just want to "redshift" it all. I want to go from a 0 - 16 MHz band to a 0 - 48 kHz band. But those pulses that pulse at 4 Hz, I still want them to pulse at 4 Hz.

The way I would try to do this, is by binning the channels. So, if I have 1024 channels (as from above) - I'd try to bin them into, say, 32 channels (and compress the spectrum by 32 times) - which basically means that I would take blocks of 32 channels and sum them or average them somehow to get one channel. Now, I really want to compress this spectrum 666 times but it doesn't matter. I really just need to understand how to meaningfully combine multiple Fourier Transforms into one.

Unfortunately, the only way I know how to do this is by taking the amplitude of the Fourier Transforms and then summing up adjacent channels to bin the channels. But this LOSES phase information - and I lose the pulses that I want to survive. So, how do I those these channels into a bin meaningful without losing both phase or amplitude information?

More concisely, my question is the following:

  Given N Fourier Transforms of N real functions, is there a meaningful
  way to get a sum of them without losing phase and amplitude
  information?

Of course, this doesn't have to be the question if someone knows a better way to handle my overall problem without binning channels.

Mathematically, shifting the frequency of a signal is pretty easy:

Following @OlliNiemitalo's answer, the 0.003 frequency shift can either be done in

time domain, or
frequency domain.

I recommend doing it in time domain, by multiplying the signal with a complex sinusoid,

$$e^{j 2 \pi f_\text{shift} n},\, n \in\{1,2,\ldots\}$$

that way you can get arbitrary frequency shifts.

Still, you're not right about needing to add up the amplitudes -- what you're doing is probably coherent sampling (i.e., the sampling clock doesn't stop after each 2048 sample chunk and starts with a random delay for the next, but continously runs at a stable frequency), so adding up multiple FFTs without reducing them to their amplitude works -- that's the beauty of the DFT being a linear operation, and the bases of the base transform it is being orthogonal.

Think of it like this:

Assume you add up many ($k$) consecutive FFTs. Because the DFT is linear, it'd be the same as adding up many consecutive 2048-vectors of time samples. What happens then is that things that repeat every 2048 samples get multiplied by $k$, whereas things that don't get attenuated or cancel themselves out. That sounds like a bad thing™, however, you have to realize one thing about the frequencies that the DFT can represent in each bin:

Every bin in the 2048-FFT represents one frequency that is an integer multiple of the fundamental frequency (ie. $\frac{f_\text{sample}}{N_\text{FFT}}$). So things that perfectly fall into each bin already inherently repeat every $N_\text{FFT} = 2048$. Hence, doing that will usually (assuming sufficiently stable sampling clock!) simply improve your SNR by $k$.

Concluding, just add up 2048-vectors of your time signal, and transform the result. You will get a 2048-DFT with an increased SNR, which definitely is meaningful.


Understanding FFT size in OFDM modulation I know that the FFT size is the the number of subcarriers OFDM system have. But I also know that the FFT has to be $ 2^N $ big (in order to process N samples from time domain to frequency domain).

But here in Matlab I can set my N to 333 for example, and The program will run without any problem (of course, the constellation diagrams will be different if I use for example N = 256, but that's just it). Why is that? Did I misunderstand something?

1) Also, If I have 1024-OFDM modulation, and the constellation on subchannels is 64 QAM, my FFT size is in the name of the modulation, right? Is it 1024?. Because 64 is the number of constellation points (size of the QAM symbols), so it can't be it.

An FFT does not need to be of length 2^N.  However, FFTs of lengths that are a multiple of only tiny primes run faster.  If you don't mind slower performance (or higher CPU utilization), you can use any length.

Note that multiple modulations can be combined, for instance, modulation of each subcarrier, and modulation of a collection of carriers.  Each can be different, thus differently named.  For instance, each subcarrier in an OFDM system can have a different size of QAM constellation.


How can one know the polynome equation based on the output of a system? I came across this situation in my textbook:

However I have no clue about how you can (starting from the stepresponse on the left), get the polynome equation on the right.

Could someone please explain?

EDIT: what about this case? http://postimg.org/image/726noo2k3/b8bf257

The right-hand polynomial is the Z-transform $X(z)$ of the discrete-time function $x[n]$, and is given by

$$X(z)=\sum_{n=0}^{\infty}x[n]z^{-n}$$

(unilateral or bilateral does not matter here, since $x[n]$ is causal anyway). So, if you now insert $x[n]$ you get

$$X(z)=\sum_{n=0}^{\infty}z^{-n}=\sum_{n=0}^\infty\left(z^{-1}\right)^n$$

You can use formula for the geometric series, as long as $|z|>1$ (i.e. $|z|>1$ is the region of convergence of the Z-transform):

$$X(z)=\sum_{n=0}^{\infty}\left(z^{-1}\right)^n=\frac{1}{1-z^{-1}}=\frac{z}{z-1}.$$

For the second case, i.e. $x[n]=a^nu[n]$, you can write the following:

$$X(z)=\sum_{n=0}^{\infty}a^nz^{-n}=\sum_{n=0}^{\infty}(az^{-1})^n$$

Then, again applying the geometric series formula you get (given that the series converges, i.e. $|az^{-1}|<1$)
$$X(z)=\frac{1}{1-az^{-1}}$$


how to find the frequency of the downsampled signal? So, I am trying to solve this question.
Consider the signal, $$x=\cos\left(\left(\dfrac{2\pi}{3} \right)n\right)$$
The signal is downsampled by a factor of two. Indicate the frequency of the resulting output, normalized by $2\pi$
(E.g., if the frequency is $\dfrac{\pi}{2}$, write $\dfrac{1}{4}$)
So, what I understood is,
from the signal, $f = \dfrac{1}{N} \implies \dfrac{1}{3}$  or $\dfrac{2\pi}{3}$ if we consider angular frequency.
So after downsampling,
$$f_d = \dfrac{f}{M} = \frac{1}{3}\cdot \frac{1}{2}  \implies \frac{1}{6}$$.
or, if taking $\dfrac{2\pi}{3}$ and downsampling by factor $2$ then, $\dfrac{2\pi}{3}\cdot\dfrac{1}{2} = \dfrac{\pi}{3}$ or $\dfrac{1}{6}$
But the answer is wrong. How should I solve it?

If the term downsampling is used in its correct meaning as anti-aliasing lowpass filtering followed by decimation, then your output should be zero, as the lowpass filter with a cutoff frequency of $\omega = \pi/2$ should remove the input $x[n] = \cos( \frac{2\pi}{3} n) $ ; i.e. , frequency of the cosine is out of the passband of the lowpass filter, and will be blocked by the filter.
However, if you simply mean decimation (or compression) by the term downsampling, then, the output will be: $y[n] = x[2n] = \cos(\frac{2\pi}{3} 2n) = \cos(\frac{4\pi}{3} n) = \cos(\frac{2\pi}{3} n) = x[n]$. This equivalence is due to the aliasing that occurs for omitting the anti-aliasing lowpass filter.
Hence, in this case the frequency of the output is the same as frequency of the input.


If your sampling frequency is 5Hz (Which I don't think is!) and the signal contains a known frequency of 5Hz. So, the sampled signal will have aliasing for sure.

You are correctly mapping your positive FFT spectrum from 0 to 2.5Hz, because digital frequency of sampled signal can only map from $[-f_s/2, f_s/2]$ which is $[-2.5hz, 2.5hz]$.

The question is where do you want to see your 5 Hz in this. You can only see a high DC which is actually the negative folded image of the 5Hz tone at $f_s$. And, you can actually see a peak very close to DC.

  Based on the Plots you have shown :

The plot of time domain data and assuming you have injected a 5Hz tone, I think the sampling frequency is close to 500Hz and not 5Hz. Simply, because you have ~2.5 cycles of 5Hz sinusoidal in 250 samples. 5Hz sinusoidal wave will take 1 second for 5 cycles, hence ~2.5 cycles of such wave will take 0.5 second. 250 samples/0.5 second means 500 samples/second. Hence, your sampling frequency is 500 Hz.

And, if you are taking 250 length FFT, frequency resolution will be of :
$$\frac{2\pi}{250} = \frac{2}{250}.\frac{f_s}{2} = 2Hz$$
Hence, your 5Hz tone will be seen as a peak between $2^{nd}$ and $3^{rd}$ coefficient of the 250-Point DFT, which you are able to see.

MATLAB Simulation showing your situation:

fd = 5;             % data tone of 0.25Hz
fs = 500;           % sampling frequency of 5Hz
Ts = 1/fs;
n_samps = 250;      % Number of samples
time = n_samps/fs;  % duration of data collection (0.5 secs)
n = linspace(0,time-Ts,n_samps); % sample points in the duration [0sec, 0.5sec)

% Carrier Frequency fc = 5Hz + a slow linear ramp as it seems from the
% time-domain plot you have attached.
signal = (n/250 + sin(2*pi*fd*n));

% AWGN Noise added
signal = awgn(signal, 15);
figure(1)
plot(signal)

Y = fft(signal);
figure(2)
plot(5/n_samps * (0:1:n_samps/2-1), (abs(Y(1:length(Y)/2))));

  The other small peaks are most probably being generated due the
  harmonics which are present as a consequence of clipping in sampled
  data due to fixed bit-width ADC. If you lower your amplitude a little
  to accommodate your sinusoidal wave within the dynamic range of ADC,
  these small peaks might go away.

Adding analysis of the time-domain samples you have pasted :

Assumption : Sampling Frequency = 500Hz

I have plotted FFT of the samples centered around 0 Hz(DC). And, based on number of samples (256) and assumed $f_s = 500Hz$, frequency spacing between 2 FFT samples is $\frac{f_s}{n_samps} = \frac{500}{256} = 1.9531$. Since 5Hz is not an integral multiple of this spacing, therefore, you will see the tone energy spread around into 1.9531*2Hz and 1.9531*3Hz i.e. 2nd and 3rd FFT coefficient.

You are also seeing a strong peak at DC (0Hz) because your sinusoid of 5Hz is having an amplitude shift of a non-zero value.

You can see that your tone is not missing at all. It is present right where is should be.

MATLAB Code to plot the signal and it's DFT with correct $f_spacing$ is as follows:

plot(500/256*(-128:127), abs(fftshift(fft(x))))


Starting off an IIR with plain averages to accelerate slow rise time A 1-term IIR can easily be started off with the first data point, instead of 0.
For data that's all > 0,
this ought to be better than starting off with y = 0 :

def update( x, t, a=0.1 ):  # data x (scalar or vector), t 0 1 2 ...
    global y  # or self.
    if t == 0:
        y = x
    else:
        y += a * (x - y)  # fraction `a` of the way towards x
    return y

It's also easy to start off with plain averages (FIR) for the first few data points,
then switch to IIR:

def update( x, t, a=0.1, naverage=10 ):
    if t <= naverage:
        a_t = 1. / (t + 1)  # start off with plain averages, 1 1/2 1/3 ...
    else:
        a_t = a       # then IIR
    y += a_t * (x - y)

The goal here is to accelerate the slow rise time of y when a is small.
Is this reasonable ?
Are there less ad-hoc ways of doing this ?

Added 28 Nov, a plot showing the difference:

You are omitting the initial condition of the filter:

$$y(t)=\alpha x(t-1)+(1-\alpha)y(t-1)$$

which is

$$y(0)=x(0)=x_0$$

which is implemented by correcting this line before the cycle:

if t == 1:
    y=a*x;

You dont need to switch anything. Only provide that value and an initial evaluation with that index, and you will not have a slow rise problem.

2016-11-22
The following image shows the result of putting the proper IC to the filter. In red there is the standard "zero" IC, while the green shows the correct IC. The $x(t)$ is a step periodic signal. Note that the IC is perfectly matched with the expected result.

2016-11-23

Another depict for a biased Laplacian Noise, and different parameters for a linear filter, with IC settled as indicated. Note that regardless of the IC setting, all the signals are around the bias, which is the correct result.


What's Logic Behind the Construction of Sobel's Filter in Image Processing? I am basically very new to this image processing field. I am presently working on edge detection on colour images. While learning the basics of edges and edge detection in images, I encountered image derivatives and spatial masks for the corresponding operations. It's where I happened to learn about Prewitt operator and Sobel operators. I cannot understand the logic behind the construction of these masks and how does it detect lines. Can someone help me please?

A first rationale is to be very short, as there was a time when computing on images was expensive. Then, a contour or an edge often present a fast variation in image intensities, that can be enhanced by derivatives. Sobel filters emulate such derivatives in one direction, and slightly average pixels in the complementary direction, to smooth small variations or noise.

One direction implements the shortest possible centered 1D discrete derivative:

$$\begin{bmatrix}
   -1 &0 &1
\end{bmatrix}
$$
to detect variations across lines, the other the shortest non-trivial Pascal/Gaussian smoothing
$$\begin{bmatrix}
   1&2&1
\end{bmatrix}^T
$$
to smooth along columns, resulting in, for instance:
$$
\begin{bmatrix}
   1&2&1
\end{bmatrix}^T\cdot
\begin{bmatrix}
   -1 &0 &1
\end{bmatrix}
$$
or
$$
\begin{bmatrix}
   -1 &0 &1 \\
   -2 &0 &2 \\
   -1 &0 &1 \\
\end{bmatrix}
$$

As you can see, this only involves dyadic numbers, so it can be implemented with adds and binary shifts.

Of course, the 3-point derivative often has an additional $1/2$ factor: $$\begin{bmatrix}
   -1/2 &0 &1/2
\end{bmatrix}
$$
to get the appropriate scale factor, and the Pascal smoother has a $1/4$ factor to have its coefficients sum to one
$$\begin{bmatrix}
   1/4&1/2&1/4
\end{bmatrix}
$$
but the resulting global scaling of $1/2\times 1/4$ does not change the edge detection power for such linear filters.


Do negative frequencies count towards null-to-null bandwidth? I am aware that negative frequencies do not normally count towards the bandwidth when transmitting in the baseband. However when specifically asked to look at the null-to-null or zero-crossing bandwidth, should the negative frequencies count? For example, in the sinc plot below, is the bandwidth 1 or 2?

It is null-null when the signal is bandpass and zero-null when the signal is baseband. So you will always take the positive frequencies.


Proving that FFT (and by extension any correlation based estimator) is unbiased I have a pretty simple question which i should have been able to answer. Just wanted to check if people here have a better solution to it.
I am trying to establish that FFT as a frequency estimator (imagine a single sinusoid + AWGN model) is unbiased. I guess we can agree that the answer to that depends on the number of grid points on the frequency axis or put in a different way it depends on how many zeros you pad to your signal. And no matter what, you will always be limited by your grid resolution. There are many papers out there which sort of interpolate using samples around the detected peak to refine the estimate and that pretty much is unbiased. This question can then be extended to other sort of correlation based estimators such as range estimation for a radar wavefrom etc.
I have a similar setting that i am working on, skipping the details for brevity but the signal model and algorithm are essentially what i described above. I plotted the resulting errors and they do come out to be unbiased. I was just wondering if i can somehow prove mathematically that the estimator is indeed unbiased. The reason i want to do that is because i eventually want to compare with the CRB.
If any papers/lectures or ideas come to mind, do share them. Thanks.

A had a long email discussion about the topic of bias in frequency formulas, aka estimators with a really sharp Grad student a while back.  It makes me hesitant to speak authoritatively on statistical matters.
I don't like the term estimator as it is too broad.
Formulas come in two basic flavors, either they are exact or they are not.  Exact means mathematically exact.  So, if you feed it the DFT of a noiseless pure tone (real or complex, different formulas for each), the calculation will result in an answer limited only by the precision of the calculation.
Other formulas are approximations.  Approximations are often times cheaper to calculate and good enough for the job.  In heavier noise, the precision advantage of a more precise formula may be rendered useless, even a disadvantage.  Approximation formulas are inherently biased.  Run a set of tones through them and maybe, like a broken clock, it will appear unbiased for certain tones.
Now run add your run-of-the-mill AWGN to the pure tone signal and test the frequency formulas.  Does the Gaussian shape pass through unskewed and still centered?
The (excellent, Thanks Julien) paper I cited in the comment didn't take quite that approach.
I would be quite interested in another, up to date, comparison with emphasis on CRLB calculation (especially on the real tone case) of all my formulas against the known standards (Start with the cited list, I think he missed Macleods.)
Here is another similar page:
http://www.ericjacobsen.org/fe2/fe2.htm
If you do endeavor to do a comparitive analysis, and there no better way to learn the alternatives, and post it somewhere, please cite it in a comment on this answer.
Frequency on bin (which seems to be what you are asking about in terms of a native bin) cases are sort of a different beast than non-whole cycle per frame tones.
The technique cited in my articles can be extended to any arbitrary number of bins.


Understanding index transformation in derivation of Fourier transform for sampling rate reduction Was going over some notes regarding deriving fourier transform equation for Sampling Rate Reduction. Reference to Notes from below link https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-341-discrete-time-signal-processing-fall-2005/lecture-notes/lec05.pdf or from Book Discrete-Time Signal Processing by Alan V. Oppenheim (2nd Edition), equation 4.75.
$$r = i + kM$$
I am lost as to how this is obtained. I understand that every $M-1$ samples are dropped from original sampling results; but still cannot understand how this expression for $r$ is derived.
Could someone help me understand this?

It's not derived, it's just chosen in a smart way such that the relationship between the decimated and the original sequences becomes obvious.
It's just a rearrangement of the terms of the sum. As a simple example, take an infinite sum of numbers $a_r$:
$$S=\sum_{r=-\infty}^{\infty}a_r\tag{1}$$
Under certain conditions that we don't need to bother with now we can rearrange the sum and write it like
$$\begin{align}S&=\ldots+a_{-2M}+a_{-M}+a_0+a_M+a_{2M}+\ldots\\&\ldots+a_{-2M+1}+a_{-M+1}+a_1+a_{M+1}+a_{2M+1}+\ldots\\&\vdots\\&\ldots+a_{-2M+(M-1)}+a_{-M+(M-1)}+a_{0+(M-1)}+a_{M+(M-1)}+a_{2M+(M-1)}+\ldots\tag{2}\end{align}$$
with some arbitrarily chosen integer $M$. So we just start with element $a_0$ and add every $M^{th}$ element, then we move to element $a_1$ and add again every $M^{th}$ element, etc. If we do this $M$ times, we've added all elements, just like in the originial sum $(1)$.
Eq. $(2)$ can be written as
$$S=\sum_{i=0}^{M-1}\sum_{k=-\infty}^{\infty}a_{i+kM}\tag{3}$$
which means that we expressed the index $r$ as $r=i+kM$.


Restriction of Fourier Transform I am currently reading Candes et. al.'s 2006 paper[1] on recovery of sparse signals from incomplete frequency samples. I am having trouble figuring out what is the form of the Fourier transform restricted to $\Omega$, $\hat{f}|_{\Omega}$, which appears in Lemma 1.2.

This restriction appears to be a bijection $\mathcal{F}_{T \to \Omega}: l_2(T) \to l_2(\Omega)$, when $|T| = |\Omega|$ and $T, \Omega$ are subsets of $\mathbb{Z}_p$, $p$ prime. Here, $l_2(C)$ denotes the space of signals with $C$ as their support.

Is there a concrete example of such a $\mathcal{F}_{T \to \Omega}$, or are we only able to prove its existence? In Terence Tao's paper[2] containing the proof of this Lemma, the proof mentions that the matrix of $\mathcal{F}_{T \to \Omega}$'s coefficients has a form similar to that of a Vandermonde matrix, i.e. it will contain elements $(e^{- 2 \pi x_i \xi_j / p})$, which makes me suspect that we can obtain concrete examples of this linear transformation.

Yes, we can explicitely construct such an example: The DFT matrix of dimension $N$ is given by
$$F=\{\exp(-j2\pi ik/N\}_{i=0,\ldots,N-1,k=0,\ldots,N-1}$$

Now, if $T,\Omega\subset\{0,1,\ldots,N-1\}$ the corresponding DFT matrix is given by

$$F_{T\rightarrow\Omega}=\{\exp(-j2\pi ik/N\}_{i\in\Omega,k\in T}$$

We can directly implement this in e.g. Python:

N = 17  # the order/dimension (P in your notation)
L = 8   # the dimension of T, Omega
# choose some random subsets of 0,...,N-1 for T and Omega
T = np.random.choice(np.arange(N), L, replace=False)
Omega = np.random.choice(np.arange(N), L, replace=False)

# Calculate the DFT matrix
t, o = np.meshgrid(T, Omega)
F = np.exp(2j*np.pi*t*o/N)

import numpy.linalg as lg
# Print its rank. It should be always equal
# to L to make F a bijection (i.e. invertible)
print (lg.matrix_rank(F))

When running above code with N=17, the rank is always 8. I've tried setting N=16, and it turns out that it can happen that the rank goes to 7 for some random configuration, i.e. F would not be a bijection.


How does a digital filter $H(z)=1/(1-z^{-1})$ change a continuous rectangular pulse? I inputted the following signal

in this discrete time filter

and got this output

Could anyone explain how ? Also how does a discrete filter process a continuous time output ?

Unless otherwise specified, your transfer function (with ROC $|z| > 1$)

$$ H(z) = \frac{1}{1 - z^{-1}}$$ defines an discrete-time accumulator with impulse response;
$$h[n] = u[n]$$ and the corresponding output
$$y[n] = \sum_{k=-\infty}^{n} x[k] $$ which reduces to
$$y[n] = \sum_{k=0}^{n} x[k] $$ for a causal input.

So practically the accumulator behaves like a continuous-time integrator and therefore for an input of square wave, you would expect a triangular wave output plus some DC term whose value depends on the initial sample timing of the square wave being applied.

Please note that this is for a AC-coupled input square wave, if the input square wave has nonzero DC component, then there will also be a ramp output term in addition to the triangular wave. However in that case you shoudl clearly indicate the time span of the input, otherwise an indefinitely long input would mean ramp output going to infinity.


How to calculate performance of your object recognition algorithm and compare results with others? I am working on a novel object recognition algorithm with BoF implementation for my MSc project and it performs fairly well. I am using it on Caltech101 dataset. I want to compare the results with other techniques, like SIFT, SURF, PHOW and so on. The problem I'm facing is that none of the papers that list out their performance in terms of accuracy mention the number of categories taken into consideration.

E.g. for the Vl_feat PHOW, the script says in the comment section that it should give you about 64% accuracy. But when I run it without any changes, I'm getting 92%. When I changed it to test on all 102 categories, it's giving me an accuracy in the 30's range.

I am fairly new to this, so if I'm missing something that's obvious, please accept my apologies.

The problem seems to be that the authors of the other papers did not include enough information for you (or anyone) to be able to reproduce their work. This could have been oversight, or the editor (or reviewers) may have removed it prior to publication. I'd suggest contacting the authors of the problematic paper and see if they'll answer your question. Be polite (it's always nice to know someone read your paper, so let them know that you have!)


matched filter of chirp signals In simulating linear frequency modulated (LFM) chirp signals, with the parameters:
Pulse width, bandwidth.

I find the matched filter response to be this:

In particular, there is a increase in the energy over the right side of the plot.

My question is, is this normal, or should the tapering be going down?

I found that the reason for the spike on the right is due to the inclusion of pulse width of time along the time axis.


The main issue is that the filters start with an initial position of (0,0). Then when you filter your actual data, the filter output has to ramp up to the actual position. In doing so, you can also observe some overshoot since your filter has a very short time response.

It would help if you were to start with a better estimate of the initial position by using for example the first position with:

dataOut_x = xcor_i(1)+filter(b,a,xcor_i-xcor_i(1));
dataOut_y = ycor_i(1)+filter(b,a,ycor_i-ycor_i(1));

The resulting plot would then look like:


Proving the upper bound of cross correlation I am reading about cross-correlation from this document and equation (5) states that

The maximum value of the crosscorrelation is not always when the shift
equals zero; however, we can prove the following property revealing to
us what value the maximum cannot exceed
$$R_{xy}(\tau) \le \sqrt{R_{xx}(0)R_{yy}(0)}$$

However, I could not find any proof for this equation. Does anyone have a resource regarding this property?

The  document cited by the OP actually states that
$$|R_{X,Y}(\tau)| \leq \sqrt{R_X(0)R_Y(0)}\tag{1}$$
which implies what the OP wants to prove, namely that
$R_{X,Y}(\tau) \leq \sqrt{R_X(0)R_Y(0)}$ but also tells us
that $R_{X,Y}(\tau) \geq -\sqrt{R_X(0)R_Y(0)}$.
Recall that for jointly wide-sense stationary processes (which are, of course, individually wide-sense stationary too and thus have constant (possibly nonzero) means $\mu_X$ and $\mu_Y$), the definitions of the functions in $(1)$ are
$$ R_X(\tau)\stackrel{\Delta}{=} E[X(t)X(t+\tau)],\quad R_Y(\tau)\stackrel{\Delta}{=} E[Y(t)Y(t+\tau)],\quad R_{X,Y}(\tau)\stackrel{\Delta}{=} E[X(t)Y(t+\tau)]$$
where, because of the wide-sense stationarity assumption, the functions depend only on $\tau$, and not upon $t$. Note also that the means $\mu_X$ and $\mu_Y$ of the random processes are not involved at all, much to the sorrow of those who prefer to demean their processes before calculating correlations and thus end up with the autocovariance  and cross-covariance functions instead of autocorrelation functions and cross-correlation functions as they are defined and understood in the communications and signal processing literature.
Now consider the quantity $$E\left[(\alpha X(t)\pm\beta Y(t+\tau))^2\right]
= \alpha^2 R_X(0)+\beta^2 R_Y(0) \pm 2\alpha\beta R_{X,Y}(\tau)$$ which, being the expected value of a nonnegative random variable, must be nonnegative for all choices of real numbers $\alpha$ and $\beta$. Thus,
$$\pm 2\alpha\beta R_{X,Y}(\tau) \leq \alpha^2 R_X(0)+\beta^2 R_Y(0)\quad \forall ~\alpha, \beta$$ or equivalently, for every possible choice of positive real numbers $\alpha$ and $\beta$, it must be that
$$2\alpha\beta |R_{X,Y}(\tau)| \leq \alpha^2 R_X(0)+\beta^2 R_Y(0). \tag{2}$$
Choosing $$\alpha^2 = \frac 12 \sqrt{\frac{R_Y(0)}{R_X(0)}},
\quad \beta^2 = \frac 12 \sqrt{\frac{R_X(0)}{R_Y(0)}}$$
so that $\alpha^2\beta^2=\frac 14$, we see that Eq.$(2)$ reduces to
$$|R_{X.Y}(\tau)| \leq \frac 12 \sqrt{R_X(0)R_Y(0)} + \frac 12 \sqrt{R_X(0)R_Y(0)}$$
which is $(1)$.


Essential Bandwidth of rect(t/T) Here is a question I have been trying to solve:

Estimate the "essential bandwidth" of a rectangular pulse

$$ g(t) = \operatorname{rect}\left(\frac{t}{T}\right), $$ with $T>0$, where this "essential" bandwidth contains 90% of the rectangular pulse energy.

What I have so far is that the Fourier Transform of $\operatorname{rect}\left(\frac{t}{T}\right)$ is

$$ G(f) = \mathcal{F}\{g(t)\} = \mathcal{F}\left\{ \operatorname{rect}\left(\frac{t}{T}\right) \right\} = T \operatorname{sinc}(fT) $$

where $$\operatorname{rect}(u) \triangleq \begin{cases}
0 & \text{if } |u| > \frac{1}{2} \\
\frac{1}{2} & \text{if } |u| = \frac{1}{2} \\
1 & \text{if } |u| < \frac{1}{2} \\
\end{cases}$$

$$\operatorname{sinc}(u) \triangleq \begin{cases}
\frac{\sin(\pi u)}{\pi u} & \text{if } u \ne 0 \\
1 & \text{if } u = 0 \\
\end{cases}$$

$$ X(f) = \mathcal{F}\{x(t)\} \triangleq \int\limits_{-\infty}^{+\infty} x(t) \ e^{-i 2 \pi f t} \ dt $$
and
$$ x(t) = \int\limits_{-\infty}^{+\infty} X(f) \ e^{+i 2 \pi f t} \ df. $$

Integrating $G(f)$ over $\pm \infty$ results in $1$. Also, integrating $|g(t)|^2$ over $\pm \infty$ results in $T$. This is about where I am lost.

Any help is appreciated.

I'll add another answer, even though MBaz's answer is correct, because I think that it doesn't actually address the problem you have with arriving at the final solution. Summarizing, you have to solve

$$\frac{1}{2\pi}\int_{-W}^W|P(\omega)|^2d\omega=0.9\int_{-\infty}^{\infty}\text{rect}^2(t/T)dt=0.9\cdot T\tag{1}$$

for $W$ (which is the bandwidth in radians). If you like, rewrite (1) using $f$ instead of $\omega$. $P(\omega)$ is the Fourier transform of the rectangular pulse:

$$P(\omega)=2\frac{\sin(\omega T/2)}{\omega}\tag{2}$$

If I understood your question and your comments correctly, then you understand all of this, but your problem is solving (1) for $W$. Note that the integral has no closed-form solution, unless you consider an expression including the sine integral $\text{Si}(x)=\int_0^x\sin(t)/t\;dt$ closed form. With a little help from WolframAlpha, Equation (1) can be written as

$$\frac{2}{\pi W}\left[TW\cdot \text{Si}(TW)+\cos(TW)-1\right]=0.9\cdot T\tag{3}$$

Equation (3) has no closed-form solution for $W$, so it is not surprising that you got confused. You could of course try to solve (3) numerically but I do not think that this is the idea of the exercise. What you need to know is that obviously the largest portion of the energy is contained in the main lobe of the spectrum, i.e. between $-W_0$ and $W_0$, where $W_0$ is the first zero crossing of $P(\omega)$. From (2) it is clear that

$$W_0T/2=\pi\quad\Longrightarrow W_0=\frac{2\pi}{T}\tag{4}$$

So this is the essential bandwidth in radians.
Checking the exact percentage using the left-hand side of (3) we get

$$\frac{2}{\pi W_0}\left[TW_0\cdot \text{Si}(TW_0)+\cos(TW_0)-1\right]=
\frac{2\;\text{Si}(2\pi)}{\pi}\cdot T=0.9028\cdot T$$

which seems close enough.

In sum, the essential bandwidth of a rectangular pulse is given by the width of  the mainlobe of its spectrum, so you only need to be able to calculate the first zero of the spectrum and you're done.


How to do this inverse Z-transform? Is this correct? It was explained to me here that I need to do an inverse z-transform on a z-based equation in order to get something I can use. That example was quite simple. But I'm not sure how to do the same process on a more complex equation.

I am starting with:

$H(s) = \frac{-ms}{ms+2R}$

Where m=mass and R=impedance. I can do a backward Euler substitution $s=\frac{1−z^{−1}}{T}$:

$H(z) = \frac{m (z^{-1} - 1)}{-m z^{-1} + m + 2 R T}$

How do I then get this in terms of $[n]$ and $[n-1]$? I understand $H(z) = H[n]$ and $H(z)*z^{-1} = H[n-1]$, but how do I mathematically accomplish this and get rid of the $z^{-1}$ terms here?

Here is my best try:

$-H(z) m z^{-1} + H(z) m + 2 R T H(z) = m z^{-1} - m$

$-H[n-1] m - H[n] m + 2 R T H[n] = m z^{-1} - m$

$H[n] (-m + 2 R T) = m z^{-1} - m + H[n-1] m$

$H[n] = \frac{m z^{-1} - m + H[n-1] m}{2 R T - m}$

And I'm stuck there with one remaining $z^{-1}$ term. Multiplying by $z^{-1}$ is supposed to be the same as delaying by a sample right? So if mass is constant, then can I just say for the remaining term $m z^{-1} = m$ and get:

$H[n] = \frac{H[n-1] m}{2 R T - m}$

Is this the right idea or if I screwed it up what was I supposed to do differently?

Thanks.

A rational transfer function

$$H(z)=\frac{Y(z)}{X(z)}=\frac{\sum_{n=0}^{N}b[n]z^{-n}}{1+\sum_{n=1}^{N}a[n]z^{-n}}\tag{1}$$

corresponds to the following difference equation:

$$y[n]=b[0]x[n]+b[1]x[n-1]+\ldots +b[N]x[n-N]-\\-a[1]y[n-1]-\ldots -a[N]y[n-N]\tag{2}$$

So the current output sample $y[n]$ can be computed from the current input sample $x[n]$ and $N$ past input and output samples. $N$ is the order of the difference equation $(2)$.

In your example you have

$$H(z)=\frac{Y(z)}{X(z)}=-c\frac{1-z^{-1}}{1-cz^{-1}},\qquad c=\frac{1}{1+\frac{2RT}{m}}\tag{3}$$

and, consequently,

$$y[n]=c\big(-x[n]+x[n-1]+y[n-1]\big)\tag{4}$$


Autocorrelation of a telegraph process/constant signal I am trying to calculate the autocorrelation function for the telegraph process, but I somehow don't get the right results.

I am calculating the correlation function of a signal X which can take values 0 and 1 by calculating the spectral density $S(\omega)=|FFT(X)|^2$. The autocorrelation function is then given by the inverse fft of the spectral density $iFFT[S(\omega)](\tau)$.

If I do this the autocorrelation functions never drops down to zero for $\tau \rightarrow\infty$.

Indeed if I take a constant signal $X = 1$, the autocorrelation should be $0$. But if I use the steps mentioned above the autocorrelation function will give $1$ instead.

I really don't understand this, could you explain why I get these results?

The random telegraph signal with parameter $\lambda$ is derived from a Poisson process
of arrival rate $\lambda$. In the Poisson process, the number of arrivals in $(t_1, t_2]$ is a Poisson random variable with parameter $\lambda(t_2-t_1)$. Consider counting the
Poisson arrivals in $(0,\infty)$ with a digital counter whose least significant bit (LSB)
is $0$ or $1$ with equal probability at time $t=0$.  The random telegraph signal is just
a model for the value of the LSB, or rather for the value of $(-1)^{\text{LSB}}$ which
has value $+1$ or $-1$ according as the LSB is $0$ or $1$.

The random telegraph signal is a random process (collection of random variables)
$\{X_t \colon 0 \leq t < \infty\}$ where each $X_t \in \{+1,-1\}$ and $X_0$ is equally likely to be $+1$ and $-1$. For any $t_2 > t_1 \geq 0$, the number of arrivals in
$(t_1,t_2]$ is Poisson with parameter $\lambda(t_2-t_1)$ and thus is an
even number with probability $\exp(-\lambda(t_2-t_1))\cosh(\lambda(t_2-t_1))
= \frac 12 + \frac 12 \exp(-2\lambda(t_2-t_1))$.
It follows that the LSB is in the same state at time $t_2$ as it was at
time $t_1$, that is, $X_{t_2} = X_{t_1}$
with probability $\frac 12 + \frac 12\exp(-2\lambda(t_2-t_1))$ while $X_{t_2} = -X_{t_1}$
with probability $\frac 12 - \frac 12\exp(-2\lambda(t_2-t_1))$. It is
readily verified that $P\{X_t = +1\} = P\{X_t = -1\} = \frac 12$ and so the
process has mean function $0$ while the autocorrelation
function is

$$\begin{align}R_X(t_1,t_2) &= E[X_{t_1}X_{t_2}]\\
&= (+1)\left[\frac 12 + \exp(-2\lambda(t_2-t_1))\right]
+ (-1)\left[\frac 12 - \exp(-2\lambda(t_2-t_1))\right]\\
&= \exp(-2\lambda(t_2-t_1))\end{align}$$

which is a function of the time difference only.  Thus

  The random telegraph signal is a zero-mean WSS process with
  autocorrelation function $R_X(\tau) =  E[X_tX_{t+\tau}]= \exp(-2\lambda \tau)$
  and its power spectral density is of the form of a Cauchy density
  function.

It is worth emphasizing that if one takes a finite segment of an
actual realization or sample path of a random telegraph signal,
one does not get the nice and smooth
autocorrelation function described above but rather a piecewise linear
function that decays away slowly.

So, can we get any of this from a DFT?  The answer is that we could
get an approximate answer, but we need to take care in setting up
the problem.

The random telegraph signal shifts back and forth between $+1$ and
$-1$ at random times.  Thus, when we take $N$ samples spaced $T$
seconds apart of a realization of this signal, each sample value is
$+1$ and $-1$,
but since this is not a band-limited signal, the samples cannot be
used to reconstruct the realization at all. Note also that unless
$T$ is carefully chosen to be quite small so that we can be virtually
certain (i.e. with very high probability)
that there is at most one transition from $+1$ to $-1$ (or vice versa),
we will have missed a lot of the nuances (a.k.a transitions)
in the signal.  But, if $T$ is very very small and $N$ is very very large,
our time-domain vector will have reasonably long runs of consecutive
$+1$s or consecutive $-1$s during which interval there are no
transitions at all.

  Now, the autocorrelation function that we need
  is the aperiodic autocorrelation function, and not the periodic
  autocorrelation function, and so we need to zero-pad our signal
  vector of length $N$ with $N-1$ or more $0$s (say $N$ for
  concreteness) and compute the DFT of this vector of length $2N$.
  An FFT algorithm could be used for this. Now  $|X[k]|^2$ will
  roughly resemble the Cauchy spectrum of the random telegraph
  signal.


Image processing for extremely small feature detection I'm new to the world of image processing, and I need some hints to analyze an image.

I have a sensor which is made up of a high resolution CMOS camera sensor and there's a mask with a micrometric hole on it.

A light shines through the mask and the pinhole creates a circular projection on the camera's CMOS. Since the projection affects only 10-20 pixels out of the 8MP of the sensor, what would be a smart way to process the image the fastest way possible? I'm now using Matlab's native code but I want to code it myself.

The sensor is self-powered, so CPU processing is an issue.

Someone suggested wavelet analysis, I was wondering if that's really my best shot.

What do you think?

I'm going to go ahead and assume currently using Matlab implies you're on a laptop or PC-alike device with a modern x86 CPU; optimal/efficient/power-saving DSP depends on the hardware you're on, but if we have the features of a fully fledged PC platform, it becomes pretty easy :)

I think this is the classical case where a linear search doesn't hurt.

Assume your pictures come in a in a row-major format, i.e. each pixel lies in memory next to its left and right neighbours, and the top and bottom neighbours are one row width of pixels away.

Since memory bandwidth is the limiting factor for simple operations on data, it's always important to minimize the times you fetch new memory – things that are in the CPU registers can directly be used for calculations, whereas things that still are in RAM easily take a couple hundred multiplications worth of CPU time to just get loaded, which is why there is an intermediate RAM cache in your CPU, which your processor can use efficiently – unless you jump wildly through memory as you process, because then it newer knows what it should keep in cache and what to throw out to make space for new stuff.

So, row-major greyscale.

You go to the first row in your picture, and search for the brightest pixel; you keep a variable storing both the value (necessary for searching) and the index of that pixel (column). At the same time, you keep a sum of pixels. At the end of the row, you store that sum and the index of the brightest pixel.

Move on to the next row, and so on.

At some point, you might want to terminate; either you've run out of rows, or your algorithm has produced enough high row sums, and isn't producing them anymore, so you've probably found your circular spot.

With the brightest pixel column indices and knowing in which rows that happened, you have your region of interest in linear time, without thrashing the CPU cache.

If this isn't fast enough for you, the nice thing about the operations you've been doing (sum, comparison with value) is that modern x86 CPUs have instructions that let you these operations with up to 16 operands at once (depending on the data type you're using underneath), so this might be a good time to look into that. I'd recommend not writing SIMD intrinsics from hand for a start, but looking into VOLK (Vector Optimized Library of Kernels) which should have the routines you need, at least for floating point, already.

Yet another thing: this can be highly scaled on multiple CPUs. If you've got four CPU cores, the first could do the upper fourth of the image, the second the next fourth and so on. Same goes for doing this with OpenCL on your graphics card – but do remember that the whole data marshalling in and out of graphic cards and GPU drivers might eat its performance children rather easily.

So, with your Region Of Interest detected, you can do multiple things. I'd just go ahead and copy the smaller subimage of which I'm sure contains the whole spot into a new, compact image, and probably store it for later verification (at least as long I'm still developing); then I'd go and let one of the classical circle-finding approaches loose on the image: Typically, the circular Hough transform will do just fine. Notice that you could have done that directly with the original image, but you'd have lost a lot of time spent on transforming things that definitely don't contain circles. The OpenCV has circular Hough built-in.

A bit of development advice: don't prematurely optimize! Get your raw ROI detection working with whatever programming language you prefer (since we're talking about most efficient implementations here, this is pretty likely to be C/C++, but needn't be! A functional programming language might make it a lot easier for the compiler and runtime to parallelize things, though it will be harder (if possible at all) to later on sneak SIMD intrinsics into the code.). You'll probably notice that pretty much everything beats the speed of the Matlab interpreter by orders of magnitude, whereas the natively implemented low-level algorithms within Matlab are highly optimized (you pretty much will have a hard time beating the speed of Matlab's FFT implementation – which is FFTw, by the way, or it's matrix/vector operations, which are BLAS/LAPACK, and so on). But practical algorithms not only do elementary operations but also have program flow defined in matlab M-code, and that, performance-wise and from the point of elegance, is stuck somewhere in the early 80's, at times. So you'd be surprised how much you often can win purely by practically writing the same program in Python, using numpy, scipy, pyOpenCV, without touching low-level coding at all.


2D DFT in image processing in python I try to compute 2D DFT in a greyscale image with this formula:

I write the code bellow with python

def DFT2D(image):
    data = np.asarray(image)
    M, N = image.size # (img x, img y)
    dft2d = np.zeros((M,N))
    for k in range(M):
        for l in range(N):
            sum_matrix = 0.0
            for m in range(M):
                for n in range(N):
                    e = cmath.exp(- 2j * np.pi * ((k * m) / M + (l * n) / N))
                    sum_matrix +=  data[m,n] * e
            dft2d[k,l] = sum_matrix
    return dft2d

I think I have some problems because I do not get the expected result. I get the image below:

I do not know where is the problem in my code!

I don't want to use a built-in function, I assume this should be work but maybe it has some mathematical issue. when I get a result matrix and convert it to image it does not as same as expected DFT image

So I think it's mostly missing the line

    dft2d = np.zeros((M,N),dtype=complex)

in your code. Otherwise your sum may not end up being complex.

And then you can't then plot the resulting image... because that needs to be real-valued.

I modified your code so that it sort of works (but is very slow). See below.  This results in something more like what you expect --- though because of the indexing, the peak is in the top-left corner rather than the middle.  There's probably something like Matlab's fftshift function that will center it for you.

Updated Code Below

import numpy as np
import matplotlib.pyplot as plt
import PIL
import cmath

def DFT2D(image):
    data = np.asarray(image)
    M, N = image.size # (img x, img y)
    dft2d = np.zeros((M,N),dtype=complex)
    for k in range(M):
        for l in range(N):
            sum_matrix = 0.0
            for m in range(M):
                for n in range(N):
                    e = cmath.exp(- 2j * np.pi * ((k * m) / M + (l * n) / N))
                    sum_matrix +=  data[m,n,1] * e
            dft2d[k,l] = sum_matrix
    return dft2d

img = PIL.Image.open("example.jpg")
img2 = img.resize((50,50))
plt.imshow(img2)
dft = DFT2D(img2)
plt.imshow(dft.real)


Is it customary to correct for the gain of a window? Consider how the Hanning window is defined:

0.5 - 0.5 * cos(n*2*Pi/(N-1))

By this definition, it has a gain of 0.5, which is simply the average value of the coefficients. By contrast, Flattop windows, as defined, have unity gain, presumably by design.

It would seem appropriate to scale the Hanning window by a factor of 2, but I have never seen this discussed anywhere. It would seem that all windows should be scaled for unity gain.

In practice, are windows typically corrected for their gain?  If not, why not?

EDIT:

Since nobody has given an answer, I'll elaborate a bit.

It is quite easy to find papers that report the gain of the more common windows. But nowhere have I seen anyone refer to correcting the gain before using it for spectral analysis. Maybe I have always missed that statement, or everyone assumes gain correction to be an obvious requirement.

It seems like common sense to set the gain of a window to unity so that signal's energy level is preserved. Furthermore, how can one compare the various windows for amplitude accuracy if one has 0 dB gain, as a flattop does, and the other has nearly 10 dB loss, as the Gauss does.

Windows are also widely used for FIR filter design. In this application, it should be clear that the signal to be windowed, a sinc pulse, has most of its energy in the center of the window. Consequently, the window does little to reduce the sinc pulse's total energy. Thus, when used for filter design, we don't want unity gain, but rather unity peak amplitude, as most windows have, except the flattops. Something other than unity peak amplitude would affect the gain of the resulting FIR filter.

Yes, it is customary to correct for the gain of a window, except for some cases I refer to later. (If you are interested only in the relative amplitude, of course you do not need to correct for the gain.)

Because the window reduces the gain of the original signal (time domain), the amplitude obtained through FFT need to be corrected. For example, if you use the Hanning window, you need to multiply all the amplitudes by 2 (the reciprocal of 0.5). As I understand it, most of the software packages for FFT automatically correct for the window used.

However, such correction is good only when all the frequencies of interest distribute throughout the time domain window. For example, suppose you have 1024 data with all signal levels of zero except for #512 point which has a value of 1 (impulse signal). Obviously, any windows do nothing to the data. So, if you correct the amplitudes for the window gain (multiply by 2), then you will end up with overestimation of the amplitude. If your 1024 data are all zero except for the very 1st point with a value of 1, then every point has a value of zero after windowing, and you lose the signal.

So, if you are dealing with random signals, with all the frequency components expected to lie almost evenly over the length of the signal, you need (or should) correct for the gain of the window you use.


Can I set a constraint on the first tap of an FIR filter such that its inverse is stable? Let's say I have the following FIR filter $h[n]$, so the output $y[n]$ for an input $x[n]$ is
$$
y[n] = \sum_{k=0}^{m-1}x[n-k]h[k]
$$

The inverse of this filter is given by the IIR difference equation

$$
y[n] = \frac{1}{h[0]}\bigg(x[n] - \sum_{k=1}^{m-1}y[n-k]h[k]\bigg)
$$

Is there some constraint I can put on the filter taps such that the inverse is stable?

I know that if I keep the zeros of the FIR filter inside the unit circle, then the poles of the inverse filter will also be inside the unit circle implying stability. However, are there more simple constraints I can make on $h[n]$ to ensure the stability of the inverse? I was thinking something like $h[0] > \sum_{k\neq0} h[k]$, but I can't prove stability of that.

The Invertible FIR Filter

A constraint based on the first coefficient alone is developed as follows:

From Cauchy's argument principle any FIR filter that meets the following constraint will be invertible (including marginal invertibility, change $\le$ to $<$ otherwise):

$$\max\left(\arg \left( H(e^{j\omega}) \right)\right)-\min\left(\arg \left( H(e^{j\omega}) \right)\right) \le \pi\tag{1}\label{1}, \space\space \omega \in [\omega_o, \omega_o+\pi) $$

Where:

$H(e^{j\omega}) = \sum_{n=0}^{N-1}h_ne^{-j\omega n}$: Frequency response of FIR filter

$\arg()$ : unwrapped phase of $H(\omega)$

Note: In an earlier version of this post, I had a simpler subset of the above that constrained the absolute value of the phase to not exceed $\pi$, but consider given any such solution that we can rotate the filter response by a fixed angle without effecting invertibility. Thus it is more generally constrained to be that for a FIR with real coefficients, the difference of the unwrapped phase for $\omega \in [0, \pi)$ cannot exceed $\pi$. Extending this to an FIR with complex coefficients, and given that we rotate the filter response (which shifts the frequency response by a fixed amount), results in the generalized constraint above applicable to any interval over $\pi$ in frequency.

Thus to constrain the first coefficient $h_o$ we can derive from $\ref{1}$:

$$\arg( h_0) + \max\left(\arg\left(\sum_{n=1}^{N-1}h_ne^{-j\omega n} \right)\right)- \min\left(\arg\left(\sum_{n=1}^{N-1}h_ne^{-j\omega n} \right)\right) \le \pi\tag{2}\label{2}$$

Which shows the complexity of a constraint based on the first coefficient alone but that it can exist. It is equivalent from this and simpler to state that if all the zeros are inside the unit circle (an invertible FIR filter), then the plot of the frequency response on a complex plane as we sweep $\omega$ from $0$ to $2\pi$) cannot encircle the origin.

Further details below:

An invertible FIR filter is a minimum phase filter, since all zeros must be inside the unit circle (or on the unit circle for marginal stability). However a subset of all possible minimum phase filters given a first tap would be the set of decreasing coefficients as the OP has hypothesized as the would meet the constraint given below under “by inspection”.  However other minimum phase filters also exist where subsequent taps are larger than the first.

The easiest constraint I can think of that would be on all coefficients specifically beyond solving for the roots is given by Cauchy's argument principle: The frequency response as given by the coefficients cannot encircle the origin for an invertible FIR filter. With the frequency response given as:

$$H(z) \bigg|_{z=j\omega} = \sum_{n=0}^{N-1}h_n e^{-j\omega n}$$

Further details below:

The invertible filter must have all zeros inside the unit circle, since zeros become poles once the filter is inverted, and any poles outside the unit circle means instability for causal systems. (For marginal stability consideration; meaning a system that neither grows nor decays, the zero can be on the unit circle.)

An FIR filter with all zeros inside the unit circle is a minimum phase filter Any other FIR filters would not be invertible. This includes maximum phase filters, which have all zero's outside the unit circle, and mixed phased filters with both minimum and maximum phase components (linear phase filters are mixed phase).

So the constraint is the filter must be a minimum phase filter to be invertible. Below I list 4 tests for detecting a minimum phase filter, with the Cauchy's argument principle the closest to providing a simple coefficient constraint rule.

By Inspection: A tell-tale sign of a minimum phase filter is the concentration of coefficients toward the start of the filter. Given this, you can rule out many filters by simple inspection if the coefficients are either concentrated in the center or end of the filter. Specifically when considering all filters with the same magnitude response, the coefficients for the minimum phase filter (which is the filter's impulse response), will decay the fastest in time. A detailed proof of this property specific to minimum phase polynomials is given in Oppenheim and Shafer's Digital Signal Processing book and is summarized as:

$$\sum_{n=0}^N|h[n]|^2\ge \sum_{n=0}^N|g[n]|^2$$

Where $h[n]$ is a minimum phase filter and $g[n]$ is any other filter with the same magnitude response and $N$ can be any positive integer. This does not mean all coefficients are in decreasing order, for example [5 6 3 2 1] is linear phase, while [5 8 3 2 1] is not, so this is not necessarily a simple constraint that can be applied but can certainly identify obvious non-linear phase solutions.

Cauchy's Argument Principle: A very simple approach to test for this condition applicable to FIR filters is to use Cauchy's argument principle (see Nyquist's Stability Criterion and Cauchy's Argument Principle) by plotting the frequency response on a complex plane. For a causal FIR filter, the number of clock-wise encirclements of the origin will be equal to the number of zeros outside the unit circle. If all the zeros are inside the unit circle, there will be no encirclements of the origin (I show an example below). For an FIR filter, there cannot be any counter-clockwise encirclements (as all poles are at the origin), so if any encirclements do happen, they will be clock-wise only.

Solving for the Roots: Confirm the magnitude of all roots of the polynomial given by the coefficients of the filter are all $|z|\le 1$

Hilbert Transform: Compare the magnitude and phase response for the filter. Since there is a unique relationship between the magnitude response and phase response for minimum phase filters, the two can be compared to determine if the filter in question is indeed the minimum phase solution for that magnitude response. This is detailed further by PeterK in this post: Derive minimum phase from magnitude with relationship copied below:

$$
\theta(\omega) = - {\scr H}\left[ \ln(G(\omega)) \right]
$$
where ${\scr H}$ is the Hilbert transform.

Every magnitude response has a minimum phase solution, therefore every FIR filter can be decomposed into a minimum phase filter (invertible) cascaded with a all-pass filter (constant magnitude response with phase change only over frequency and not invertible).

These concepts are demonstrated with a simple example of a 2 tap FIR filter with coefficients [1 0.5] and it's reverse [0.5 1]. In the first case the filter is minimum phase with a transfer function $1+0.5z^{-1}$ and the second filter is the reverse, which is a maximum phase filter with the transfer function $0.5+z^{-1}$. The magnitude response of both filters is identical but the phase response is very different, as given by the vector diagram showing phase versus frequency for both filters. (This diagram is created by simply replacing $z^{-1}$ with the phasor $e^{-j\omega}$ which is exactly how to get the frequency response from the z-transform).

It is a bit of an optical illusion but the resulting phasor given by the sum $\sum_n h[n]e^{-jn \omega}$ does have the exact same magnitude for all frequencies as we sweep the frequency $\omega$ from $0$ to $2\pi$. However notice how constrained the angle will be in the diagram on the left, compared to the one on the right! Minimum phase versus Maximum phase.  This plot also demonstrates Cauchy's argument principle showing how the frequency response of a minimum phase filter when plotted on the complex plane will not encircle the origin.

Below is a plot of the magnitude and phase response for the two example filters above (with the same magnitude response for either). Since group delay is $d\phi/d\omega$ the minimum phase filter will have the lowest delay while the maximum phase filter will have the largest delay, which makes sense when you consider the placement of the largest tap in the FIR filter which is a series of summed and weighted delays (the energy will emerge from the minimum phase filter sooner).

Another example demonstrating the Cauchy argument principle is below, with a plot of the frequency response on the complex plane for filters [5 6 3 2 1] and [5 8 3 2 1]. The filter [5 8 3 2 1] is proven to not be minimum phase since the frequency response encircles the origin:

And here is another example for a maximum phase filter, which is also confirmed using Cauchy's Argument Principle. This is for the fourth order filter with coefficients [1 -3 -3 2 5] where we see all four zeros are outside the unit circle as we have four encirclements of the origin. (The easy way to count encirclements is to note a direction on the frequency response with a forward direction consistent with increasing $\omega$, and then draw a vector from the origin out toward infinity at any angle and count how many crossings of the frequency response take place: if the cross is of a forward direction the count increases, and if of a negative direction the count decreases).

And here is another simple example with the FIR filter given by coefficients [1, 1,5, 0.6] showing how the frequency response of an FIR with a real positive first coefficient can enter the LHP and still be a minimum phase filter. Specifically, it is the origin that is not encircled consistent with Cauchy's Argument Principle. With that constraint the phase response can never exceed $\pm \pi$. Below the plot is the standard magnitude and phase frequency response as two separate plots.


Normalisation of different length signals I have a collection of audio files of subjects sleeping, which have different lengths,  ranging from 8 hours to 10 hours.

I want to detect when a snore occurs in them.  First of all I want to normalise them.  What is the best method to normalise them so they all have similar values.

Currently I'm am using the following code which basically sets the sum of all the values to 1.

recording = recording / norm(recording, 1);

However, for example, the energy value of a snore differs greatly from one recording to another.

Is there any better methods of normalization which will cope with a range of different length recordings?

In audio, "normalization" means making all values be between 0 and 1. So

recording = recording / max( abs(min(recording)), abs(max(recording)) );


Kalman Filter : How measurement noise covariance matrix and process noise helps in working of kalman filter , explain intuitively please? How noise covariance matrix and process covariance matrix helps in improving the state estimate, can some one explain intuitively without mathematics ?

Considering a linear dynamic system framework, if there were no noises in the process equation and measurement equation, then one could use deterministic means to compute the state exactly, without any need for an estimator.

However, the unavoidable presence of noise in the state transition and the measurements make it impossible to compute the state exactly but only approximately. As the measurements are not true but in error.

Kalman filter (similar to other statistical estimators) improve the state computation by referring to the statistical character of the noise present in the measurements and state transtitions. Kalman filter also benefits from the dynamic equation of the system in predicting the state update.

So those noise covariance matrices are used to specify the Kalman gain K, which tries to correct the error in the deterministic computation of the state update.


How to Design Nonlinear Phase FIR Filter? I've recently encountered a rather clever low-pass filter that appears to have been designed to sacrifice uniform group delay in the pass-band to get better frequency response relative to the filter length. What methods are available to design FIR filters that maximize frequency performance relative to length while allowing phase to become nonlinear (but still continuous) in the pass-band?

Although not a comprehensive list, the following are filter algorithms available in Matlab for designing non-linear phase FIR filters:

firlpnorm: Nonlinear phase equiripple or least squares, uses Lp norm: for equiripple filters p = infinity, for least squares p =2.

cfirmpm: complex and nonlinear phase equiripple filters.

See the Matlab documentation for full details on how to design non-linear phase using these algorithms.