FFT of square wave: maths and actual numbers do not coincide this fact is driving me insane... Theory states that the Fourier transform of a square wave is a sinc function, but if I compute the fft of a synthetic (perfect) square wave of amplitude 1 and frequency 1 and I check the numbers, I see them drift from a sinc function. Phases seem to change linearly thru the bins, instead of being all constant as expected, and real values (the odd ones of course) are all strangely equal to 2/framesize instead of being zero as I would expect. In fact, if I synthesize the same square wave in frequency domain by simply using a sinc function and I do the IFFT and I check the output, I notice that the square wave so obtained is almost perfect if it were not for some tiny ringings at the edges (Jibbs ?)
How is it ? What is then the "correct" function representing it in frequency domain if not a complex sinc as it should be ??? Or perhaps what I noticed is simply due to the fact that, since we are dealing with discrete quantities (fft =dft) things work slightly different than in the abstract case (fourier transform in pure math sense), and therefore some little modifications (scaling ?)are required ?  Thanks in advance

The continuous-time Fourier transform of a single rectangular pulse $p(t)$ of duration $[-d,d]$ is :

$$ P(\Omega) = \frac{ 2 \sin(\Omega d) }{\Omega } = 2 d \cdot \text{sinc}( \frac{\Omega d}{\pi} ) \tag{1} $$

where $\Omega$ is the frequency in radians per second.

You can not represent $p(t)$ or $P(\Omega)$ using a sampled-data computer system because $p(t)$ is not bandlimited. However, there is an analogous definition of a rectangular pulse in discrete-time, defined as a sequence $x[n]$ of duration $[-d,d]$, with $d$ being integer, whose discrete-time Fourier transform is :

$$X(e^{j\omega}) = \frac{ \sin( \omega (2d+1) /2 ) }{ \sin( \omega /2 )} \tag{2}  $$

which is called a periodic-sinc or a Dirichlet kernel. Note that this is not a simulation of the continuous case.

Using the computer with N-point DFT/FFT, implementation of discrete-time case will be

$$X[k] = X(e^{j\omega})|_{w_k = \frac{2\pi}{N}k} = \frac{ \sin( \frac{\pi k}{N} (2d+1)) }{ \sin( \frac{\pi k}{N} )} \tag{3}  $$

A Matlab/Octave implementation is as follows

d = 3;       % rectangular pulse in [-d,d]
N = 32;      % DFT/FFT length

k = 0:N-1;   % DFT index

X = sin( (pi/N)*(2*d+1)*k)./sin( (pi/N)*k));
X(1) = 2*d+1;

x = real(ifft(X,N));   % compute inverse DFT on Xk

figure,stem(k,xp);
title('rect-pulse from inverse DFT of periodic-sinc');

Note that the circularity of DFT implies that negative indices of the pulse $x[n]$ are at the right edge.


How to get a real signal back to a non-complex form after processing it in frequency domain (using fft) I have an audio signal that is composed of 2 channels. The signal is all real numbers.

I transform this signal into frequency domain using fft, multiply it with a filter frequency response, and then transform it back into time domain.

My problem is that the result is a complex signal.

How can I convert this complex signal back into a real time signal.

The filter is not a real filter.

First I would advise against your filtering approach, which is the "Frequency Sampling" method of filtering (if you are using multiplying your FFT by a target response directly) which has poor performance and efficiency. (The frequency sampling approach will provide your exact solution at your FFT bin centers only, and then a lot more ripple in between vs the algorithms I suggest below.)  However if you need to be in the frequency domain anyway, and are multiplying by the FFT of your desired filter coefficients, this would be ok as long as you are properly dealing with the circular convolution involved.

Also I would use a real filter unless you are intentionally trying to get an asymmetric spectrum response (meaning positive and negative spectrums are differrent, which implies a complex signal-- given that is what you are trying to avoid then there is no reason to use a complex filter).

To implement a real filter, with a real signal and a real output, consider using either the Parks-McClellan or Least Squares algorithms, with design tools readily available in Matlab, Octave and Python, and then implement your filter as an FIR filter in the time domain.

If you need to implement filtering in the frequency domain, you can force it to be a real filter by ensuring that the filter is conjugate symmetric; for an FFT the 0 bin is the center (DC value), and then the samples from 1 to N/2-1 should be conjugate symmetric to the samples from N/2 to N-1; where sample 1 is the conjugate of sample N-1, sample 2 is conjugate of sample N-2 etc...


How can I get the continuous-time transfer function coefficients (or poles and zeros) from the corresponding discrete-time TF and vice versa? Let's say I have a continuous time transfer function which I know its numerator coefficients $(B^c = [b^c_m, ..., b^c_1, b^c_0])$ and denominator coefficients $(A^c = [a^c_m, ..., a^c_1, a^c_0])$. Also I know its zeros $(Z^c = [z^c_m, ..., z^c_1])$, and poles $(P^c = [p^c_m, ..., p^c_1])$ and gain $K^c$.

How can I get the numerator $(B^d = [b^d_m, ..., b^d_1, b^d_0])$, denominator $(A^d = [a^d_m, ..., a^d_1, a^d_0])$, zeros $(Z^d = [z^d_m, ..., z^d_1])$, poles $(P^d = [p^d_m, ..., p^d_1])$, and gain $K^d$ of an equivalent discrete-time transfer function (considering sampling time of $T_s$)

I know about the transformation methods (like Tustin or bilinear approximation), but those transformation methods convert the transfer function from continuous to discrete and it is not obvious from that how the numerator and denominator coefficients change. In other words I need the implementation of two functions like below:

$$[A^d, B^d] = \text{function_1} \,(A^c, B^c)$$
$$[Z^d, P^d, K^d] = \text{function_2} \,(Z^c, P^c, K^c)$$

I know for function_1 there is a Matlab function called c2d() which does that, but that is not compatible with Matlab coder, and I need to write my own function. So that is why I am interested to know the explicit mathematical formulas that does the conversion between coefficients (or zeros and poles) and not the transfer function.

For the function_2 maybe I have the answer but I am not sure it can handle all cases (including single/multiple poles or zeros):

$Z^d = e^{Z^c \times T_s}$

$P^d = e^{P^c \times T_s}$

$K^d = 1-e^{-K^c \times T_s}$    ---> (is the minus sign right?)

This is the formula that you need:
$$
H_d(z) = H_a \left( j\frac{2}{T} \tan(\omega_d T/2) \right),
$$
where $H_a(s)$ is the transfer function in the Laplace domain, and $\omega_d$ is the "warping frequency", i.e. the frequency where you want the transfer functions to line up exactly (usually chosen as the cutoff frequency or center frequency).


What information does the abs(normalized_image - normalized_quantised_image) carry? Given a grayscale image X2

And its (under)quantised version with 3 bits depth X3

What information does the abs(X2_normalized - X3_normalized) carry?As an example the following image X5 generated with the code below

X2norm= X2-min(X2(:)); X2norm = X2norm / max(X2norm(:));
X3norm= X3-min(X3(:)); X3norm = X3norm / max(X3norm(:));
X5 = abs(X2norm - X3norm);
image5=figure;imagesc(X5);title('image 5');
saveas(image5,'image5.jpg');

Assuming your image is a $K$ bits depth image, and the quantized version is a $3$ bits version, pixels of the difference image, represent $(K-3)$ least significant bits of pixels of the original image. It means it completely is similar to original image in darker parts and completely masks the brighter parts.


Low pass filter with DFT This is my first question at this site, please excuse me if I'm not asking the question properly.

I have an acoustic signal of a test and I want to remove the "thin" peaks, because they don't really affect the structure that I'm studying. For that I've thought of a Low pass filter in order to retain only frequencies below 300 Hz.

My aim is to count the peaks of the filtered signal above 3 sigma (Standard Deviation).
I'm programming in C# and I don't know if there is a direct (and simpler way) of doing this but what I'm doing now is

Discrete Fourier Transform to my signal (SPL/s -> SPL/Hz)
rescale what I obtain (bins = sampling frequency)
Cut frequencies above 300 Hz
Inverse DFT

It seems logical to me but I'm not sure if I'm doing it correctly.
I have this little problem while doing the FT, is that I don't know how to scale properly the magnitude axis (y): if I transform and then I do the inverse transform I don't get the original signal but an attenuated one. What's happening here?

And my other question is about the symmetric nature of a FT. If I do a Fourier Transform to my signal I get a mirrored spectrum that goes till 8000 Hz, can I simply cut at my desired frequency (is what I'm doing now)? Thanks for your time.

I have an acoustic signal of a test and I want to remove the "thin"
  peaks, because they don't really affect the structure that I'm
  studying. For that I've thought of a Low pass filter in order to
  retain only frequencies below 300 Hz.

If this is accurate, that is, you have an acoustic signal that exists from baseband to 300 Hz, then yes, what you want to do is a digital low-pass filtering operation to retain all energy below 300 Hz, and discard energies above that.

  My aim is to count the peaks of the filtered signal above 3 sigma
  (Standard Deviation). I'm programming in C# and I don't know if there
  is a direct (and simpler way) of doing this but what I'm doing now is

Your intuition is correct - you want to remove frequencies below 300 Hz, so naturally you think to transform your signal into the Fourier domain, clobber unwanted frequencies, and then transform back. There are some issues with this method:

Firstly, this method is taboo in the DSP community, (and for decent reasons), although there are cases where you can get away with it, depending on your application. Some background: Nulling out unwanted co-efficients in a transform domain and inverse transforming back into the original domain is used in denoising especially in wavelets, and it also makes intuitive sense, (and correctly so). The problem with doing this in the Fourier domain is that the basis function of the DFT are global, and not local like they are in wavelets. In other words, you can clobber our wavelet co-efficients and affect only local areas in your signal (ie, small parts of your signal) - but doing so in the Fourier domain affects all your signal.
Expanding on the above background, when you "boxcar" your frequency domain, (ie, mask it with 1s and 0s that either allow or disallow certain frequencies), this is the same as convolving your time-domain signal with a sinc function, that carries on forever. This is why simply culling out unwanted frequencies - while acceptable in some applications - will introduce 'ringing/smearing' in the time domain.
Instead what is usually done, is that people create filters (masks in the frequency domain that are not sharp 1s and 0s, but more graceful), such that their equivalent transform back into the time-domain does not have the effective infinite extent of sinc function, (that is, the decay of the sidelobes is much faster than that of a sinc). This then provides a nice trade-off between removing unwanted frequencies, and not having to pay the price of extreme ringing/smearing in the time domain.

Armed with this task at hand, what you need to do is the following:

Create a digital low pass filter in some other software such as MATLAB or Python. You can create a simple FIR filter for your purposes. (Do not use an IIR because you are working with an audio application and default IIR filters will warp your phase, but thats another story). When you make your filter, you will simply have a vector of numbers, of a length you chose. The larger the length of the filter, the sharper the transitions in the frequency domain, but also the larger the delay, so this is a trade-off you can experiment with. At the end of the say though, you will $1$ by $N$ vector of numbers. This is your filter. Let us call it $h[n]$
You have your audio signal, let us call that $s[n]$. The planets are not aligned, so most likely your audio signal is of a different length than your filter, and let us say it is of length $M$. So your signal $s[n]$ is simply a $1$ by $M$ vector of numbers.
Your objective is to filter $s[n]$ by $h[n]$. There are a number of ways to do this. Since you already seem to have access to the FFT in C#, we will use the FFT method. In the FFT method, you want to multiply the Fourier transform of your filter, (call it $H(f)$), with the Fourier Transform of your signal, (call it $S(f)$).
We doing a linear convolution, (as opposed to a circular one), but in the frequency domain, so we have to take extra care when it comes to FFT sizes. If you did the linear convolution in the time domain, the resulting filtered signal would be of length $N + M -1$. This is then the FFT size you need to use.
You are almost done. All you need to do now is take the $N+M-1$ length FFT of your filter, $h[n]$. Take the $N+M-1$ length FFT of your signal, $s[n]$. Multiply them together. Then take the $N+M-1$ length IFFT of this element-by-element product, and viola, you have filtered your signal.

Hope this helped.


Help finding the output for input with conditions I'm trying to find the output where $h[n]$ is the impulse response of a real stable causal LTI system, and $s_0 \epsilon \mathbb{R}$, and where the input is $x[n] = e^{2\pi i n s_0}$ for $n \epsilon \mathbb{N}$.  Supposedly my answer should be a compact formula with no infinite sums...

So since since my system is causal I obtain, $y[k] = e^{2\pi i ks_0}\sum_{n=0}^{\infty}h[n]e^{-2\pi i ns_0}$.  From here I can conclude that since the system is stable this series converges, but without knowing more information about $h[n]$ I can't see any way to further simplify the expression and rid myself of this infinite sum.  Of course I also know by assumption that $y[k],s_0 \epsilon \mathbb{R}$, but I can't see anyway to use this fact without, again, knowing more about $h[n]$.

I suppose my sum is looking a lot like the standard Fourier transform but I don't see how I can use that with out, once again, knowing more about $h[n]$.  Hopefully someone can help me here, thanks.

Edit: Could it be that a real system implies that $h[n]$ is real?  Although that's still not enough info to allow me to find a closed form for this series.

Edit 2: Or could "real" in this context be referring to a system which is actually implemented?

So since since my system is causal I obtain, $y[k] = e^{2\pi i ks_0}\sum_{n=0}^{\infty}h[n]e^{-2\pi i ns_0}$.  From here I can conclude that since the system is stable this series converges, but without knowing more information about $h[n]$ I can't see any way to further simplify the expression and rid myself of this infinite sum.

Hint:  The function (of $z$) defined by $\sum_n h[n]z^{-n}$ is given a special
name and is of some importance in signal-processing circles. Do you know about
this function?
What is the value of this function at $z = e^{2 \pi i s_0}$? which by the way
lies on the unit circle in the complex plane, seeing as how we are talking of circles.


Fast/efficient way to turn real passband samples to complex baseband I have real-valued samples of a bpsk signal i am receiving (see frequency spectrum in figure A), that I want to turn into complex samples for further processing.

However when I just use the received signal as the real part of my complex symbols, I get a mirrored version of the original (like in figure B) with two times the bandwidth of A. By shifting the signal to baseband, and then downsampling it to half the samplerate I get something like in figure C.

Question: Is there a more efficient way to get directly, i.e. in one step from A to C?

In the particular case of BPSK transmitted over wireline, there is no quadrature component, so the equivalent low-pass representation (the complex envelope) is real. You can get the complex envelope simply by downcoverting the signal from the carrier frequency down to DC.

For general baseband signals, or over the wireless channel, the process can be outlined as follows:

Find the analytic signal. Say the passband signal is $r(t)$. Then, $r_+(t)=r(t)+j\hat{r}(t)$ has no negative frequencies ($\hat{r}(t)$ is the Hilbert transform of $r(t)$).
Downcovert $r_+(t)$. The signal $x(t)=r_+(t)e^{-j2\pi f_c t}$ is the complex envelope.

Note that implementing this algorithm is not necessarily trivial. In practice, you almost always want to downconvert $r(t)$ to some intermediate frequency before attempting the Hilbert transform.

An alternative procedure is as follows:

Center the positive part of the spectrum around DC, by calculating $y(t)=r(t)e^{-j2\pi f_c t}$. Note that the resulting signal is complex even in the BPSK case.
Get rid of the negative part of the spectrum by low-pass filtering $y(t)$.

The downside of the second procedure is that it doubles the maximum frequency of the signal. Mathematically both procedures are identical (up to a constant amplitude factor).


why companding is used? Which one is the best answer?

Companding is used []

a. to overcome quantizing noise in PCM
b. in PCM transmitters, to allow amplitude limiting in the receivers
c. to protect small signals in PCM for quantizing distortion.
d. in PCM receivers, to overcome impulse noise.

I think c is the best answer. am i right?

Some thoughts about the other possible answers:

a. in general, it can mitigate quantization noise, but cannot completely overcome it
b. the limited range is over the whole channel, not only the receivers.
d. impulse noise is not suppressed

Answer c. is apparently the closest match.
Did you check the definition from the  companding tag at SE.DSP? Could you propose some improvement?


Is the Auditory Spectrogram the same as the Magnitude Spectrogram? This is in the context of speech signals. The magnitude spectrogram is computed by framing the signal into short time windows, applying a hamming (or similar) window, computing the fft over each window, and taking the absolute value of the complex numbers obtained. I came across the Auditory Spectrogram in a paper and several toolkits (http://amtoolbox.sourceforge.net/amt-0.9.5/doc/general/audspecgram.php, http://www.isr.umd.edu/Labs/NSL/Software.htm). Do they refer to the same thing? Or else, how is the auditory spectrogram computed and what are the differences between the two representations?

Edit: As Jazzmaniac pointed out, the below answer I've given is incorrect. See Jazzmaniac's comment for explanation!

The first source you provide contains this image as an example output of the function audspecgram():

This very much appears to be a standard spectrogram of an acoustic signal. To describe it as an 'auditory spectrogram' might just hint to that fact—that it's derived from an acoustic signal. In other words: you are right in suspecting that what is described as an 'auditory spectrogram' in your source is just a spectrogram that is a result of exactly the procedure you describe to obtain it.


Estimating Filtered Noise Variance Assuming I have samples of a colored noise.
I know the empirical Variance of those samples.
I also know the bandwidth of the filter which colorized the noise (Only the bandwidth, not the actual filter).

Now I want to estimate the variance of the noise if I filter those samples by a (Different) given filter.

Is there a way to estimate or bound from above and below the variance of the noise?

Thank You.

I assume that you have a wide-sense stationary discrete-time input noise signal $x(n)$, a linear time-invariant filter with impulse response $h(n)$, and an output noise process $y(n)$. If we assume that $x(n)$ is zero mean (which is not necessary, but makes it easier) and if we model $x(n)$ as perfectly band-limited, we get for its variance

$$\sigma_x^2=\frac{1}{2\pi}\int_{-\pi}^{\pi}S_x(e^{j\theta})d\theta=
\frac{1}{\pi}\int_{\theta_1}^{\theta_2}S_x(e^{j\theta})d\theta\tag{1}$$

where $S_x(e^{j\theta})$ is the power spectral density of $x(n)$, $\theta$ is the normalized frequency $\theta=2\pi f/f_s$, and $\theta_1$ and $\theta_2$ are the lower and upper band-edges, respectively. With the frequency response of the filter

$$H(e^{j\theta})=\sum_{n=-\infty}^{\infty}h(n)e^{-jn\theta}$$

we can write the power spectral density of the output noise as

$$S_y(e^{j\theta})=S_x(e^{j\theta})\left|H(e^{j\theta})\right|^2$$

and its variance is

$$\sigma_y^2=\frac{1}{2\pi}\int_{-\pi}^{\pi}S_y(e^{j\theta})d\theta=
\frac{1}{2\pi}\int_{-\pi}^{\pi}S_x(e^{j\theta})\left|H(e^{j\theta})\right|^2d\theta=
\frac{1}{\pi}\int_{\theta_1}^{\theta_2}S_x(e^{j\theta})\left|H(e^{j\theta})\right|^2d\theta$$

From this equation we can get a lower and upper bound for $\sigma_y^2$ as follows

$$\min_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|^2
\frac{1}{\pi}\int_{\theta_1}^{\theta_2}S_x(e^{j\theta})d\theta\le\sigma_y^2\le
\max_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|^2
\frac{1}{\pi}\int_{\theta_1}^{\theta_2}S_x(e^{j\theta})d\theta$$

and from (1) we finally get

$$\min_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|^2\sigma_x^2\le
\sigma_y^2\le\max_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|^2
\sigma_x^2\tag{2}$$

These bounds can be useful (i.e. relatively tight) if the noise bandwidth is small compared to the bandwidth of the filter $h(n)$. Unfortunately, they are almost useless if the filter bandwidth is small compared to the bandwidth of the input noise process $x(n)$.

In the latter case, i.e. if the filter bandwidth is much smaller than the noise bandwidth, you can come up with similar bounds as in (2). You just need to estimate the minimum and maximum of the power spectral density $S_x(e^{j\theta})$ within the filter bandwidth.


Special case peak/notch filter for exactly 1/32 sampling rate? I need a filter at 1/32 my sampling rate. Because it's 1/2n, I was thinking there might be a special case filter for this case. Anything like this?

The period, half-period, and quarter-period of the peak or notch frequency are integer multiples (32, 16, and 8) of the sampling period. If you want a linear-phase filter then compared to an arbitrary frequency this enables savings.

When using a windowed finite impulse response (FIR) filter, every 16th coefficient will be zero. Small savings.

When using a truncated infinite impulse response (truncated IIR, TIIR) and Goertzel algorithm based filter, if the phase of the frequency is required to be zero at the center of the impulse response then the length of the impulse response can be chosen such that the phase of the frequency will be an integer multiple of $\frac{\pi}{2}$ at the ends of the impulse response. Such phase shifts of a TIIR marginally stable resonator can be implemented without multiplication that would be needed in case of an arbitrary frequency requiring an arbitrary phase shift to obtain a symmetrical impulse response.

The following band-pass composite TIIR filter has an impulse response that is a raised cosine window function of integer length $N+1$ (counting also the zero-valued endpoints) normalized to unit area and multiplied by a cosine at $1/M$ sampling frequency. $N$ must be even and an integer multiple of $M$ (which needs not be integer). The impulse response is symmetrical around time $N/2$.

$$\text{mid}_0[t] = \text{in}[t] - \text{in}[t-2] - \text{in}[t-N] + \text{in}[t-N-2]\\
\text{mid}_1[t] = \text{mid}_0[t] + 2\cos(\frac{2\pi}{M})\text{mid}_1[t-1]-\text{mid}_1[t-2]\\
\text{mid}_2[t] = \text{mid}_0[t] + 2\cos(\frac{2\pi}{M} + \frac{2\pi}{N})\text{mid}_2[t-1]-\text{mid}_2[t-2]\\
\text{mid}_3[t] = \text{mid}_0[t] + 2\cos(\frac{2\pi}{M} - \frac{2\pi}{N})\text{mid}_3[t-1]-\text{mid}_3[t-2]\\
\text{out}_\text{band}[t] = \frac{-1^\frac{N}{M}}{N}\left(\text{mid}_1[t] - \frac{1}{2}(\text{mid}_2[t] + \text{mid}_3[t])\right).$$

The filter takes advantage of that multiplication by a window function is frequency domain convolution and a raised cosine has a 3-term frequency domain representation. The filter can be easily modified to implement other generalized cosine windows. The three middle filter stages are marginally stable and require periodical refreshing for numerical stability. For example for $M = 32$ and $N = 256$:

Figure 1. Impulse response (horizontal axis: time in samples)

Figure 2. Magnitude frequency response (dB, horizontal axis: frequency in radians)

Because the filter is linear-phase, it is easy to convert it to a notch filter or a peaking equalizer filter:
$$\text{out}_\text{notch}[t] = \text{in}[t-\frac{N}{2}] - \text{out}_\text{band}[t]\\
\text{out}_\text{peaking}[t] = \text{in}[t-\frac{N}{2}] + a\,\text{out}_\text{band}[t],\ a \ge -1.$$

The following slightly faster reacting version of the filter has no zero as the first or last sample of the impulse response of length $N-1$ and requires that $N$ is an even integer that is an integer-plus-half multiple of $M$. The impulse response is now symmetrical around time $\frac{N}{2}-1$, so conversions to $\text{out}_\text{notch}$ and $\text{out}_\text{peaking}$ must be modified accordingly.

$$\text{mid}_0[t] = \text{in}[t] + \text{in}[t-N]\\
\text{mid}_1[t] = \text{mid}_0[t] + 2\cos(\frac{2\pi}{M})\text{mid}_1[t-1]-\text{mid}_1[t-2]\\
\text{mid}_2[t] = \text{mid}_0[t] + 2\cos(\frac{2\pi}{M} + \frac{2\pi}{N})\text{mid}_2[t-1]-\text{mid}_2[t-2]\\
\text{mid}_3[t] = \text{mid}_0[t] + 2\cos(\frac{2\pi}{M} - \frac{2\pi}{N})\text{mid}_3[t-1]-\text{mid}_3[t-2]\\
\text{out}_\text{band}[t] = -1^{\frac{N}{M}-\frac{1}{2}}\left(\frac{2\sin(\frac{2\pi}{M})}{N}\text{mid}_1[t] - \frac{\sin(\frac{2\pi}{M} + \frac{2\pi}{N})}{N}\text{mid}_2[t] - \frac{\sin(\frac{2\pi}{M} - \frac{2\pi}{N})}{N}\text{mid}_3[t]\right).$$

For example with $M = 32$, $N = 272$, and $\frac{N}{M} = 8.5$:

Figure 3. Impulse response (horizontal axis: time in samples) with a time span from 0 to 270 inclusive.

The frequency response looks pretty much the same as that of the previous filter, except for that the first zero is not at frequency 0. So I did not bother to plot it.

Very similar filters are presented in:

Eric Jacobsen & Richard Lyons. DSP Tips & Tricks: The Sliding DFT. IEEE Signal Processing Magazine, March 2003, pp 74–80.

Eric Jacobsen & Richard Lyons. DSP Tips & Tricks: An Update to the Sliding DFT. IEEE Signal Processing Magazine, January 2004, pp 110–111.


Combining Kalman Filters Say I want to use Kalman filters for predicting the price of items at a supermarket. I have a Kalman filter for each item (apple/beef/brooms/etc). I notice that some items are sort of related, like apple with orange, beef with pork.

My question is how can I make these Kalman filters interact with each other? For example let's say I get observations about the price of beef and they show that they are going up. I want to tell the Kalman filter for pork somehow that it should also update. How can I do this?

Any literature or examples of this usage will be greatly appreciated, thanks

EDIT:

I am actually trying to ask this question in the general sense, but I wanted to frame the problem in a way that would be simpler for others to understand. So I understand that the Kalman filter uses observations to update its state. Suppose I want to incorporate observations for some Kalman filter into another, is there an accepted way of 'transforming' these observations? Or is it more of an art?

I know of no "standard" way of incorporating the output of a set of Kalman filters into another. The reason is that each application will have a different "signal model" (way of combining the multiple outputs into different measurements).

For example, suppose you have two Kalman filters that give you state estimates:
$$
\mathbf{x}_k^1 = \mbox{state of first Kalman filter}\\
\mathbf{x}_k^2 = \mbox{state of second Kalman filter}\\
$$
then you'll need to decide the parameters in your new model:
$$
\mathbf{x}^{\tt tot}_k = [ {\mathbf{x}^1_k}^T {\mathbf{x}^2_k}^T ]^T\\
\mathbf{x}^{\tt tot}_{k+1} = \mathbf{A} \mathbf{x}^{\tt tot}_k + v_k\\
\mathbf{z}^{\tt tot}_k = \mathbf{C} \mathbf{x}^{\tt tot}_k + w_k
$$
where $\mathbf{A}$ is the state update matrix, and $\mathbf{C}$ is the output matrix, $v_k$ is the process noise driving the interactions between the different states and $w_k$ is the measurement noise.


Basic Questions on Wiener Filtering I read a lot about Wiener Filters (focusing on discrete time case). I understand the math, but I am quite disconnected from the real life assumptions behind using such a filter. Unfortunately textbooks are just laying out the formulas without putting forward much intuition in it.

The task is to extract x(t) given y(t), where y(t) = x(t) + n(t). Great but first we need to build an appropriate filter. At this point:

Are we assuming that we have access to some training samples of x(t) and/or n(t) in addition to y(t)?
Or are we just assuming that we have access to statistical properties(auto correlation, cross correlation) of x(t) and/or n(t) in addition to y(t)? If so since x(t),n(t) are not explicitly available, are these statistics computed from other signals available to us with similar statistical properties?
Or is it a combination of both, such as having access to some training samples of x(t), but having only access to second order statistics of n(t) (or vice versa)?
Is the noise n(t) always assumed to be uncorrelated to x(t) in Wiener filtering, or are the textbooks going for the simplified version? What is the default scenario?
Why are all the examples coming up with all-zero filters but not using pole-zero filters? Is it just a matter of simplified textbook treatment or is the problem definition not amenable to using transfer function modelling approach? Yes the first approach is computationally easier, but these days we have significant computing power and I don't think it is a good justification anymore, if that's the only reason.

Looking forward to hear from practitioners.

The task is to filter x(t) when given y(t), where y(t) = x(t) + n(t).
  Great but first we need to build an appropriate filter. At this point:

No. The task is to filter $y[n]$ to achieve $x[n]$. Your filter acts on the available signal $y[n]$ to get a best estimate of $x[n]$ from a noisy observation of it. BTW. discrete-time uses $y[n]$ or at least $y(n)$ notation, and not $y(t)$. In DSP, $t$ is reserved for continuous-time.

  1) Are we assuming that we have access to some training samples of
  x(t) and/or n(t) in addition to y(t)?

trainig samples is a term related with adaptivity, pattern recognition, neural networks etc. In the classical Wiener context, you do not need training samples but you need data statistics. See answer 2 for how to obtain the statistics.

  2) Or are we just assuming that we have access to statistical
  properties(auto correlation, cross correlation) of x(t) and/or n(t) in
  addition to y(t)? If so since x(t),n(t) are not explicitly available,
  are these statistics computed from other signals available to us with
  similar statistical properties?

Yes. In the Wiener filter context, you either assume (for theoretical development) or estimate (for practical applcations) the correlations (or Power Spectral Densities) between all those signals. In practice you estimate them from avaiable data by several means, depending on the application. You need data for this; and that's probably what you call training samples.

  3) Or is it a combination of both, such as having access to samples of
  x(t), but having only access to second order statistics of n(t) (or
  vice versa)?

You do not have access to relevant noise-free $x[n]$ by any means; otherwise why filter $y[n]$ to get $x[n]$? But, in some cases, one can have noise free irrelevant blocks of $x[n]$ samples, in addition to the relevant noisy blocks. So you can use those noise free blocks to estimate certain parameters of the process $x[n]$, so as to filter out the noisy blocks.

  4) Is the noise n(t) always assumed to be uncorrelated to x(t) in
  Wiener filtering, or are the textbooks going for the simplified
  version? What is the default scenario?

Yes. The noise $n(t)$ (or better $v[n]$ in discrete-time context) is almost always modeled as additive, white, zero-mean Gaussian noise uncorrelated with $x[n]$. This simplifies all the mathematical development and in typical cases is a good model. But you can relax some of those assumption, if that would yield a better (but more complex to derive) filter...

  5) Why are all the examples coming up with all-zero filters but not
  using pole-zero filters? Is it just a matter of simplified textbook
  treatment or is the problem definition not amenable to using transfer
  function modelling approach? Yes the first approach is computationally
  easier, but these days we have significant computing power and I don't
  think it is a good justification anymore, if that's the only reason.

IIR filters are good for their computational efficiency (and in Wiener context, they also yield the minimum MSE compared to FIR) but they may suffer from unstability and causality problems, and their design and analysis is more complex in adaptive filtering context. For this reasons, it's simpler and safer to use FIR based filters for adaptive Wiener applications.

To learn more about the Wiener filters and adaptive Wiener filtering:

Simon Haykin - Adaptive Filter Theory.
Monson Hayes - Statistical Digital Signal Processing and Modeling.


FFT Result - which bins are useful Which bins of FFT result are useful? For example I made FFT on 1024 real samples, and it produced result with 1024 complex samples.
I read somewhere, the only half of that result is useful if input samples are of type real.

Is bin number 513 (nyquist frequency) useful? Should I cut this bin or it should stay?

I believe you should use the Nyquist bin every time. By having it, you can always recreate the full DFT of your real-valued signal from the first half of spectrum. Whereas if you discard this value there is no way back and you can recreate your time domain signal.

As in answer you pointed to, usually the Nyquist bin shouldn't contain any data when while using anti-aliasing filter. Nevertheless sometimes you might want to estimate the energy of your signal from it's DFT. You might therefore get inaccurate value. Question is, do you really care about such accuracy?

Just to wrap up. I always use Nyquist bin and I don't need to ask myself these questions.


Inverse Fourier of Two-Pole Transfer Function I would appreciate if someone could walk me through this derivation.

I have a transfer function in the frequency domain, which has two poles

$$\tilde{H}(\omega) = \Big(\frac{1}{1 + i \omega \tau_1}\Big)\Big(\frac{1}{1 + i \omega \tau_2}\Big)$$

For a single pole, the inverse Fourier transform is:
$$\mathcal{F}^{-1}\Big\{\Big(\frac{\tau}{1 + i \omega \tau}\Big)\Big\} = \Theta(t) ~ e^{\displaystyle -\frac{t}{\tau}} $$

where $\Theta(t)$ is the Heaviside function since this is a causal filter. (See this answer.)  I also plugged the above $\tilde{H}$ into Wolfram Alpha and got the expected result:

$$\mathcal{F}^{-1}\Big\{\tilde{H}(\omega)\Big\} = \frac{1}{\tau_2 - \tau_1} \Bigg( ~ e^{\displaystyle -\frac{t - t_d}{\tau_2}} - e^{\displaystyle -\frac{t - t_d}{\tau_1}} \Bigg)\Theta(t - t_d)$$

where $t_d$ is the time after which the causal filter is effective.  I'm not sure if it's a reversal of integration limits, integration by parts, or simply some of the terms cancelling after integration and plugging in the limits that might be responsible for the subtraction.

And I'm also not 100% sure how to work in the $t_d$ term.

The comments to the question already hint at the solution, but I'll write it out.

$\tilde{H}(\omega)$ can converted from a product into a sum using Partial Fraction Expansion:

$$\tilde{H}(\omega) = \frac{1}{1 + i \omega \tau_1} \frac{1}{1 + i \omega \tau_2} = \frac{A}{1 + i \omega \tau_1} +\frac{B}{1 + i \omega \tau_2}$$

$$ 1 \omega ^0 + 0 \omega ^1 = A(1 + i \omega \tau_2) + B(1 + i \omega \tau_1)$$

Which gives
$$\implies A + B = 1$$
and
\begin{equation}
\begin{split}
A i \omega \tau_2 + B + i \omega \tau_1 &= 0\\
(1 - B) \tau_2 + B \tau_1 &= 0\\
B(\tau_1 - \tau_2) &= -\tau_2\\
B = \frac{\tau_2}{\tau_2 - \tau_1}
\end{split}
\end{equation}

$$A = 1 - B = \frac{- \tau_1}{\tau_2 - \tau_1}$$

Leading to:

$$\tilde{H}(\omega) = \Big(\frac{1}{\tau_2 - \tau_1}\Big) \Big[ \frac{\tau_2}{1 + i \omega \tau_2} -\frac{\tau_1}{1 + i \omega \tau_1}\Big]$$

$\tilde{H}(\omega)$ also contains a term $e^{-i\omega t_d}$, so now for each of the terms in the sum above the inverse Fourier transform is:

$$\mathcal{F}^{-1}\Big\{\Big(\frac{\tau}{1 + i \omega \tau}\Big) e^{-i\omega t_d}\Big\} = \Theta(t) ~ e^{\displaystyle -\frac{t - t_d}{\tau}} $$

and their sum has the desired form in the question.


Where can I find an authoritative (peer-reviewed or textbook) reference to sampling-induced beating? I presume we are all here well aware about foldback aliasing when sampling signals above the Nyquist frequency; i.e. half the sampling rate.

By contrast, the phenomenon of beating occurs when sampling a signal with a frequency slightly below the Nyquist frequency. You can see and hear some beats on this educational site.

Here is an example of a 60Hz signal sampled at 128Hz resulting in a 4Hz beat tone with a 8Hz beat. The beat tone frequency is calculated as follows: $f_{t}=\frac{f_{s}}{2}-f_{m}=\frac{128}{2}-60=4\,Hz$. The beat has a frequency twice this tone frequency.

In preparation of a paper, I am looking for a peer-reviewed or textbook reference about beating. I browsed through almost 20 DSP textbooks and found nothing.

On the web, I only found a Proceedings of the 2003 American Society for Engineering Education Annual Conference & Exposition session paper by Kostic which included the nice explanatory figure shown below. I am looking for something more authoritative though; peer-reviewed or textbook.

I think you are mixing two things that are actually not related. "Beating" happens if you add two sine waves that are close in frequency. What you describe is sampling sine wave close to the Nyquist Frequeny. If you plot the samples, it looks like there is beating going on, but that's not actually the case. All information is properly preserved and if you were to properly reconstruct the analog signal you would see a perfectly good sine wave with no beats at all.

The somewhat tricky issue here is "proper" reconstruction. Since your signal is very close to the Nyquist frequency you need a very steep anti-aliasing filter.

Here is a code example of how to reconstruct the sine wave with no beats. Any residual beating is simply an artifact of an insufficient anti-aliasing or interpolation filter.

%% create a 60 Hz sine wave sampled at 128 Hz
fs = 128;
n = 256;
x = sin(2*pi*(0:n-1)'*60/fs);
figure(1); clf;
plot(x);
% this looks like beating but it isn't

% let's upsampe to 32kHz so we can actually listen to it. We need a very
% steep anti-aliasing filter, so we go with a 64 tab Kaiser Window and a
% beta of 15
y = resample(x,250,1,64,15);
% plot a few periods from the steady state portion
figure(2); clf;
plot(y(10000:20000,:));
% looks like a perfectly good sine wave
sound(.9*y, 32000);
% sounds like one too.


Discrete inverse Fourier transform I have a question regarding discrete inverse Fourier transform, and no answer I found on the internet seem to be satisfying. This might be because I do not fully get some of them, so please excuse my ignorance.

Suppose you have full knowledge of a frequency-domain function $\hat{f}(\omega)$. It is a continuous, well-behaved bump function. You want to compute its inverse Fourier transform $f(t) = \Re \int_{-\infty}^{+\infty} \hat{f}(\omega) e^{j\omega t} d\omega$, but this cannot be done analytically, so you compute a discrete inverse Fourier transform. You know the interval $ \left[ \omega_{\min}, \omega_{\max} \right]$ on which it is reasonable to numerically compute the IDFT, so you now have to sample $\hat{f}(\omega)$ on this interval.

My question is the following: is there an equivalent of the sampling theorem for the inverse Fourier transform ? Is there a generic rule to sample my frequency-domain signal, as there is for time-domain ? I noticed that I obtain something similar to aliasing when I undersample $\hat{f}(\omega)$.

Side note : all of this is used in a finite element code to compute a complex source. As such the IDFT will be performed at least hundreds of thousand of times, so the question to limit its cost.

Thanks in advance !

The Fourier Transform is (mostly) symmetric in both directions so the same sampling rules apply.

The minimum frequency resolution is determined by the length in the time domain, the same way that the minimum time resolution is determined by the "length" in the frequency domain.

That creates a bit of theoretical problem: in order for the required frequency resolution to be finite, the length in the time domain has to be finite as well. However, a signal that has finite length in one domain has infinite length in the other. So you can't actually sample a signal in both domains at the same time.

In practice, you need to work around this by replacing "must be zero outside the supported interval" by "must be small enough outside the supported interval". So you accept some residual aliasing in either domain, provided you can make it small enough for your specific application


How to find a maximal possible value of a cross correlation Is there any way of calculating a maximal possible value of a cross correlation of two signals whose maximal amplitude and sample rate is known?

To make things clearer, , signals in question are mostly the same signal but one is shifted in phase.
So, basically I want to know what will be the output of a cross correlation if they are not shifted, but without actually shifting one back and doing the CC to find out.

Yes, but only if you know the mean and variance (or standard deviation) of the signal. The maximum correlation will be the (variance of the signal + the mean-squared) times the total number of samples.   If you do not already know the variance, calculating it would be similar to calculating the auto-correlation in terms of processing so nothing is really saved in that case. But in some cases you may have prior knowledge or ability to more easily estimate the variance (from a smaller sample set for example if the process is white and stationary).

To check my logic, consider what is involved to calculate correlation for a sampled signal when it is aligned with itself (for the maximum autocorrelation) compared to calculating the variance for the signal:

To correlate two signals: You multiply the two signals sample by sample, and then sum the result.

The equation to calculate the variance is similar, but you must subtract the mean from each sample before doing the product, and then you must divide by the total number of samples. (Important to note, this is not the unbiased estimate of the variance which scales by N-1, but the calculation for the 2nd Moment about the mean, which scales by N). The variance therefore is just a normalized autocorrelation at lag=0 with the mean removed.

To put very simply, CorrMax is just 1/N * the second moment. (The second moment is the average of the squares or mean-squared value). Variance is the 2nd moment with the mean removed (about the mean), so we have to add it back in as we've done above if we were starting with that particular factor known. However consider if we already know the mean-squared value, (or the square of the rms value since rms is commonly used), then the relationship is the mean-squared value times the number of samples.

Therefore in summary:
To determine the maximum correlation for a sequence X that is N samples long,  assuming you already have either the mean and variance (2nd moment about the mean), use:

$$ CorrMax(X) = N(\sigma_X+\mu_X^2)$$

Where:

X is a vector of length N

$ \sigma_X $ is the second moment of X about the mean

$ \mu_X $ is the mean of X

Or alternatively calculated from the mean-squared value for the sequence,

$$ CorrMax(X) = N(X_{rms}^2)$$

Where:

X is a vector of length N

$ X_{rms}^2 $ is the second moment of X (mean-square of X)


Shift property of the Short-term Fourier transform (STFT) Properties of STFT (Short-term Fourier transform) say that it preserves time shift up to modulation. Does it mean it is sensitive to time-shift in signal?

Of course the STFT changes when the signal is shifted. According to the definition we have

$$\text{STFT}\{x(t)\}=X(\tau,\omega)=\int_{-\infty}^{\infty}x(t)w(t-\tau)e^{-j\omega t}dt\tag{1}$$

If we define a shifted signal $y(t)=x(t-T)$, the corresponding STFT becomes

$$\begin{align}\text{STFT}\{y(t)\}=Y(\tau,\omega)&=\int_{-\infty}^{\infty}x(t-T)w(t-\tau)e^{-j\omega t}dt\\&=\int_{-\infty}^{\infty}x(u)w(t-(\tau-T))e^{-j\omega u}e^{-j\omega T}du\\&=e^{-j\omega T}X(\tau - T,\omega)\tag{2}\end{align}$$

From $(2)$ you see that time shifting results in a phase factor and a shift in the STFT. The magnitudes of the STFTs are related by just a time shift, which is of course identical to the time shift of the signal:

$$|\text{STFT}\{x(t-T)\}|=|Y(\tau,\omega)|=|X(\tau-T,\omega)|\tag{3}$$


Calculating spectral coherence coefficient (scc) - implementation help I am trying to calculate the SCC of a cyclostationary signal. I have estimated the Spectral Correlation Function (scf) of the signal using FFT Accumulation Method (FAM). I implemented the following pseudocode in matlab:

function [Sx, alphao, f] = scf(x, fs, d_f, d_alpha)
    Compute N' = pow2(nextpow2(fs/df))
            L = N'/4
            P=pow2(nextpow2(fs/dalpha/L));
    N = P*L;    %total no of datat points to process
    x = x(1:N)
    Divide x into chunks of length N' with L overlapping points
    Apply Hamming window across each chunk
    Perform FFT
    Perform the downshift to baseband
    Multiply the baseband with it's complex conjugate
    Execute smoothing by means of P-points FFT.
end

where fs is sampling frequency, d_f is frequency resolution and d_alpha is cyclic frequency resolution

Now I have three matrices : $S_x^\alpha$, which is the SCF, $alphao$, which is a 1D vector the cyclic frequencies and $f$, which is a 1D vector the spectral frequencies. I need to calculate the Spectral Coherence Coefficient which is given by the following equation:
$C_x^\alpha(f) = \frac{S_x^\alpha(f)}{\sqrt{S_x^0(f-\alpha/2)S_x^0(f+\alpha/2)}} \forall f$

My questions are:

If I put $S_x^0(f-\alpha/2)$ in the denominator then won't the denominator essentially becomes $S_x^0(f)$ since $\alpha = 0$?
Will $C_x^\alpha$ be of the same size as $S^\alpha_x$?

I am just working through this information myself, but I believe I understand your questions.

  If I put $S_x^0(f-\alpha/2)$ in the denominator then won't the denominator essentially become $S_x^0(f)$, since $\alpha = 0$?

It will become $|S_x^0(f)|$ at $\alpha = 0$, but at other values of $\alpha$, you would need the conjugate multiply at offsets of $f$.

  Will $C_x^\alpha$ be of the same size as $S^\alpha_x$?

Yes, one would think. The coherence SOF function, is just supposed to be some normalized version of the SCF.

I am struggling with the FAM, it appears to output data on a tilted spectral plane, where $\alpha=0$ lies on the matrix diagonal.


FFT on non-rectangular part of image I need to do a FFT on an image for noise reduction, but the problem is that I do not need the complete image, but only a circle in the middle. The borders are a fixed rig, thus I am not interested in what it displays, but it has an impact on the result of the FFT.

Is there any way to just cut out a circular part and use that for the FFT? Note that if I use black background, the edge between background, and image data will have quite an impact.

Instead of having a hard edge between the image data of interest and the background, you could use a two-dimensional tapered window function, as is often done in spectral analysis. You might start by trying a Gaussian window, which for a two-dimensional case would look something like:

$$
w[x,y] = e^{-\frac{\left(x-\frac{N_x-1}{2}\right)^2}{2\left(\sigma_x \frac{N_x-1}{2}\right)^2}} e^{-\frac{\left(y-\frac{N_y-1}{2}\right)^2}{2\left(\sigma_y \frac{N_y-1}{2}\right)^2}}
$$

$N_x$ and $N_y$ are the dimensions of the desired transform in the $x$ and $y$ directions, respectively, and $\sigma_x$ and $\sigma_y$ are parameters that allow you to control the shape of the window; for small $\sigma$ values, most of the energy in the window function will be concentrated toward the center, with that effect decreasing as you increase $\sigma$.


Correlation Performed by Convolution Background: The question here is related to images in particular and not signal/waveforms.

I have been reading a lot of answers about the difference between convolution and correlation but I am stuck on my doubt as to how are they related to each other while performing operations. Most of the answers are individual answers for each of these terms.

For example:

A statement I am trying to understand says

  Here correlation operations are performed by replacing the exhausted
  convolutions with element-wise multiplications using Discrete Fourier
  Transform (DFT).

I am not able to relate the individual definitions here. Can anyone please explain this?

What that statement is saying that

a. Correlation is performed the same way as one would perform convolution (you must implicitly know that one of the sequences is conjugated and time reversed to express a correlation as a convolution, as it was not stated there),

and

b. Mathematically the convolution is performed using this relationship

$$x[n]*h[n] = \mathscr{F}^{-1}\left\{\mathscr{F}\left\{x[n]\right\}\cdot\mathscr{F}\left\{h[n]\right\}\right\}$$

Where $*$ denotes convolution, $\mathscr{F}\left\{f[n]\right\}$ denotes the Discrete Fourier Transform of $f[n]$, and $\mathscr{F}^{-1}\left\{F[k]\right\}$ denotes the Inverse Discrete Fourier Transform of $F[k]$.

Note that these operations with Discrete Fourier Transforms actually yield the "circular" convolution of $x[n]$ and $h[n]$, which is equal to the linear convolution of $x[n]$ and $h[n]$ only under a certain condition: there must be enough zero padding to avoid aliasing.


Is it fast time or slow time for the delay in the ambiguity function domain for radar signal? I am confused about the delay concept in the ambiguity function. The FFT of slow time is the doppler frequency. However, people usually relate the range-doppler domain with the ambiguity function domain. However, the range corresponds to the fast time.

Your question is a little bit unclear, but I'm going to try to answer according to my understanding of your question:

So it sounds like you're asking about pulse-Doppler radar processing, where we have two time dimensions: fast-time and slow-time.  Let us review those terms really quick.

In pulse-Doppler processing, multiple pulses are sent out by the radar at a steady rate/frequency (PRI/PRF, pulse repetition interval/frequency).  Each time a new pulse is sent out, the radar receives signals for a period of time, and then it switches back to a transmit operation and starts the cycle over again.  The data that is received is stored in a 2D matrix.  Each time the radar receives a new pulse, a new row is added to the matrix, and the data that is recorded at each receive interval is placed in this row.

This means that, like you said before, we have two time dimensions: fast-time and slow-time.  Fast-time corresponds to a small-scale time measurement of radar data at each pulse; slow-time is a (typically) much larger time scale, and corresponds to the dimension where each pulse was sent out.

Note that in this explanation so far, we have not taken the DFT of anything, so our 2D matrix represents time and time.  Like you alluded to before, if you take the DFT of the slow time dimension's IQ data, what you will get is a transformation to the Doppler domain. This happens because over the slow-time dimension, your target signal's phase will be changing if it is moving.  Assuming the target is moving at a constant rate, this means that the phase of the target signal will linearly change, and a liner change in phase is equivalent to oscillation at some frequency; this is  why we use the DFT for this: it allows for an easy way/convenient way to measure phase.

So now that we've reviewed that, lets talk about ambiguity functions:

Think of how range-Doppler data recording is structured.  There are multiple transmit/receive options, and in order to keep the target's phase stable, every transmit pulse is (typically) exactly the same.  Therefore, how can I tell whether a target signal appeared in range during my reception of the first pulse, or if the target appeared at any other range?  This is the concept of range ambiguity.  When you look at an ambiguity function, typically the x-axis is labeled as Tau, representing a time delay: this is equivalent to the range dimension.

Similarly, let's think about how Doppler is being measured: with pulse-Doppler processing, we're counting on being able to adequately measure the rate of change of the phase in question using a set sampling interval (our PRI/PRF from earlier).  It turns out that our PRF/2 is our effective Nyquist frequency for our slow-time sampling operation  Well, just like in sampling say, digital audio signals, we have to worry about aliasing any signals we record.  For audio signals, if my Nyquist frequency is, for example, 1kHz, any signals with frequency content above 1kHz will be aliased to look as if they are at another frequency.  This is exactly what can happen in the Doppler domain when using pulse-Doppler processing: if I have target velocity that is greater than my PRF, my DFT cannot accurately represent it, and that target signal "appears" somewhere else in the Doppler domain ambiguously.  This is called Doppler ambiguity.

So how does the range-Doppler map equate to the ambiguity function? Well, they've both effectively measuring the same thing!  The ambiguity function is a measurement of Doppler speed and range/time.  The range-Doppler map (RDM) from pulse-Doppler processing is a measurement of... Doppler speed and range/time!

There are two types of ambiguity functions: single pulse, and multi-pulse.  A multi-pulse ambiguity function simple shows you the "full" picture of both range/Doppler that you cannot see on the RDM due to measurement ambiguities.

Now a single pulse ambiguity function is also a 2D image, but we have no slow-time dimension.  Instead, our Doppler dimension is derived by using multiple hypotheses of the target speed and baking our guesses in to the pulse-compression operation that gets us our range dimension.

So hopefully all of that helps answer your question a bit.  If you'd like to learn more and have access to MATLAB, try walking through the examples on this MATLAB page


FIR estimator for IIR system Suppose that we have a dynamical system of which the impulse responses are infinite (IIR). Now I found methods on papers (http://dx.doi.org/10.1109/9.839942) estimating states or outputs of such a system with a FIR estimator. So this FIR estimator gives me approximately the same (state) output as my original system.

I find it hard to conceive how an estimator with only poles in $z_i=0$ is able to duplicate a system with poles $|z_i|<1$ (but not necessarily $z_i=0$). I think it has something to do with interpolation conditions, but I would really appreciate it if someone could make this clear to me. Thanks :)

I understand that for a stable dynamical system IIR will decay to zero and thus can be truncated to FIR. But my question remains how a FIR could 'mimic' the poles of an IIR?

Simply truncating the impulse response of an IIR filter will give the optimal $\ell_2$ FIR approximation to the IIR system.

Whether this form of approximation is good enough to "mimic" the IIR system depends on what you mean by that.


Fitting a gaussian image using opencv I'm reading through the opencv documentation and some questions in SO but it doesn't seem to provide this information.

I've an image $I(x,y)$ and I want to find a gaussian function $f_{\mu,\Sigma}(x,y)$ such that
$$
f_{\mu,\Sigma}(x,y) \approx I(x,y)
$$
using LSE or MLE estimation. $\mu$ is the mean vector and $\Sigma$ is the covariance matrix. My question is if there's any opencv function that perform such task.

People suggest to implement the estimator from scratch, though some suggested to use LM algorithm for such purpose, and in general doesn't seem very difficult I wonder if there's any function for such purpose.

Is there any?

Thank you


Orthogonal basis of signal space and the projection of white noise I am trying to make the relation between linear modulation and nonlinear one by using orthonormal expansion. The purpose is to understand what is orthonormal set in each case and to understand the operation of projection.

Linear modulation

In digital linear modulation, the baseband signal can be written as

$$x(t) =  \sum_n a_n p(t-nT)$$
where $a_n$ is data, $\lbrace p_n(t) = p(t-nT)\textrm{, }n\in \mathbb{Z}\rbrace$ is orthonormal set, i.e. $<p_n(t),p_m(t)> = \delta_{m-n}$ with $<x(t),y(t)> = \int x(t)y^*(t)dt$ is inner product defined for the signal space we are working on. To get $a_n$ back, we need project $x(t)$ to $p_n(t)$.

Use matched filter $q(t) = p^*(-t)$ thus $q(t-\tau) = p^*(\tau - t)$.

by defining $g(t) = p(t) \star p^*(-t)$ where $\star$ is convolution. $g(t=0) = 1$ and $g(t=nT) = 0$ because $\delta_{m-n} = <p_n(t),p_m(t)> = g((m-n)T)$.

\begin{align}
y(t) &= x(t) \star p^*(-t)\\
&= \sum_n a_n p(t-nT)\star p^*(-t)\\
&= \sum_n a_ng(t-nT)\\
\implies y(t=mT) &= \sum_n a_ng((m-n)T) = a_m
\end{align}

If we have white noise $n(t)$ : $y(t) = x(t) + n(t)$. The noise is processed as

\begin{align}
z(t) &= n(t) \star q(t)\\
&= \int n(\tau)q(t-\tau) d\tau\\
&= \int n(\tau) p^*(\tau -t) d\tau\\
\implies z(t=mT) &= \int n(\tau) p^*(\tau - mT) d\tau\\
&= <n(t), p(t-mT)>
\end{align}

Then we can say that for a given $y(t)$, each processing $m$ at the receiver which ends with sampling at $t=mT$, the noise process is projected to the vector $p_m(t)$ creating $z_m$ and the data part is projected to the same vector creating $a_m$. $\lbrace p_m(t) = p(t-mT)\textrm{, }m\in \mathbb{Z}\rbrace$ is orthonormal set, thus $z_m$ are uncorrelated, with Gaussian assumption, $z_m$ are independent.

This is famous process in classical technical books. Two conclusions I get

The set of orthonormal vectors $\lbrace p_m(t) = p(t-mT)\textrm{,
    }m\in \mathbb{Z}\rbrace$ has infinite number of elements as $m \in
    \mathbb{Z}$.
The orthonormal vectors $p_m(t)$ can be no-time-limited and the orthogonal property holds on the entier time support.

Non-linear modulation

Let's take a binary FSK modulation as example.
$M$ basis functions, $k \in \mathbb{S} = \lbrace 1, ..., M \rbrace$:

$$\phi_k(t) = \begin{cases}\sqrt{\frac 2T}\sin\left(\frac{k\pi t}{T}\right)& \text{for}\quad 0 \leq t \leq T\\0 &\text{otherwise}\end{cases}$$

I write the waveform as $x(t) = \sum_n \gamma_n(t-nT)$ where $\gamma_n(t) = \phi_k(t) \textrm{ if } a_n = k$.

Questions

Question 1 : At the receiver, for the $n^\rm{th}$ symbol, we try projecting to $M$ functions $\phi_k(t-nT)$. If $k=k_0$ gives the largest power, we decode $a_n = k_0$. Is $\lbrace \phi_k(t-nT) \textrm{, } n \in \mathbb{Z}\textrm{, }  k \in \mathbb{S}\rbrace$ the orthonormal set as $<\phi_k(t-nT),\phi_l(t-mT)> = \delta_{(m-n)\times(k-l)}$ ?

I mean I expect something equivalent to the two properties (infinite number of element and basis functions are not time-limited) of orthonormal set in the linear modulation case.
Question 2 : is it correct if I say the noise samples are uncorrelated because the projection are $<n(t),\phi_k(t-nT)>$ and $<n(t),\phi_l(t-mT)>$ and $<\phi_k(t-nT),\phi_l(t-mT)> = \delta_{(m-n)\times(k-l)}$ ?
Question 3 : The $M$ basis functions are interpreted as constellation of $M$ points. Is there any analogy with the linear modulation case ? If we consider M-QAM, we have also a constellation of $M$ complex points but if I write $a_n = a_{rn} + ja_{in}$, the two functions $p(t-nT)$ and $jp(t-nT) = e^{j\pi/2}p(t-nT)$ are not orthogonal.

Leave time dimension aside for a moment and deal with just one symbol at a time. In a PAM modulation the generated functions are not ortho and not normal. Actually, in a 4-PAM modulation you have four basis pulses: $p_0(t)=a_0 \cdot p(t)$, $p_1(t)=a_1 \cdot p(t)$, $p_2(t)=a_2 \cdot p(t)$ and $p_3(t)=a_3 \cdot p(t)$. The pulse $p(t)$ is explicitly defined just because it makes the math easier and coherent: it is a function that generates a 1-dimensional Hilbert subspace.

The task is about mapping data symbols to signals in the Hilbert space, in such a way that at the receiver side the mapping can be reverted, thus, such a mapping must be bijective. On the receiver side, we must compute some kind of metric over the symbol period to come back to the discrete domain. For a Hilbert space, such metric is the inner product. So, in both cases (FSK and PAM) you end up performing, one way or another, the same computation:

\begin{equation}
a = \int_T f(t) g(t) dt
\end{equation}

And this inner product is performed against each orthogonal basis function (one or more). By choosing more than one orthogonal basis functions, you get more parallel channels (more dimensions) to transmit onto if you wish. Then, choose how to encode the data:

BFSK encodes data by signaling one of those channels or the other. So for one of them, $a \neq 0$ (gets mapped to symbol $s_1$) and for the other (assuming orthogonality) $a = 0$ (mapped to symbol $s_0$). If both channels get signaled, you cannot reverse the mapping.
2-PAM encodes data by signaling, within the same channel, different levels. So that if $a \gt threshold$, symbol $s_1$ is decided. If two levels are sent, you cannot reverse the mapping.
M-QAM uses both methods: different levels and different channels within each symbol period.

In short, data can be encoded in the value of $a$ within the same channel (same dimension, like a scaled basis function) so you further have to define decision boundaries for mapping $a$ to discrete symbols, or can be encoded by which channel (or dimension the orthogonal basis function generates) you signal. This resembles vector-space math because Hilbert space is the space of continuous vectors.

So, for the questions:

Q1. Yes. Harmonic sines within the fundamental period are orthogonal under the inner product. Orthonormal just means the norm for each basis function equals 1.

Q2. No. When it is said that noise is uncorrelated it refers to the fact that AWGN has no memory (time dimension), the noise is already uncorrelated before projection onto any basis. Within a symbol period the noise is correlated, its value is its power over that period.

Q3. When the constellation is said to have M points, it means M points get mapped to M different functions in the chosen Hilbert sub-space. The way you choose those functions to be is not defined. You can choose a M-dimensional space (so you only need to find the orthogonal basis function which fits better) or a N-dimensional space with $M \gt N$ and define different regions in that space for symbol mapping. Also, in QAM, $p_r(t) = cos(\omega t) \ne sin(\omega t) = p_i(t)$, and the inner product is zero within one period.

I avoided to talk about the time dimension (when possible) because IMHO it would not help. Unless you are trying to exploit temporal correlation (viterbi, equalizers, noise shaping, etc...) it does not help and only introduces noise in the explanation. Usually, this kind of analysis is performed on a symbol based period. However, according to Fourier analysis, every function in the Hilbert space can be expressed in terms of amplitude, frequency and phase, and a shift in time is a shift in phase, thus the same logic applies. In real implementations you do not find neither infinite signal levels ($a$'s), since power is constrained, nor frequencies (basis function's shape) which are constrained by bandwidth. Time shift (phase) is not infinite neither since you need to start and stop sometime.


Moving average filter magnitude attenuation I am dealing with a moving average filter which filters power system signals to rid them of noise or harmonics etc. The fundamental or useful component of these signals is centered around 50 Hz or 60 Hz. The magnitude response of the filter at an angular frequency say  as per some standard papers is given as
$|H(e^{(j\omega_1)})|=|\frac{\sin(L_m \omega_1 /2)}{L_m \sin(\omega_1/2)}|$ ,
 where $L_m$ is the window length of the filter and the value of $\omega_1=\frac{2\pi f_f}{F_s L_m}$ . Sampling frequency is given by $F_s$
whereas the input signal frequency to the filter is given by $f_f$ .  The idea is when the filter input signal frequency is equal to the fundamental (no noise or harmonics etc)  or even close to it (say 59.5 Hz or 59.7 Hz considering a 60 Hz system) the magnitude response (gain) of the above mentioned filter should be equal to one  or close to it for any sensibly chosen value of the window length and the sampling frequency. For example if i choose a value of $L_m$=31 and and consider the input signal frequency $f_f$ to be 59 Hz, with a sampling frequency $F_s$  of say 1440 Hz, i end of up with a gain of 0.9972, then if i use $f_f$=58, i end up with a gain of 0.9973, similarly along the lines, for $f_f$  equal to 50 Hz i get a gain of 0.9980, showing that gain increases as frequency deviates from the nominal (which does not make any sense). I am having a hard time understanding this, could you clear this confusion for me. Could you show me how the gain is close to unity for small frequency deviations from the nominal

Question : How is the gain close to unity in the vicinity of 60Hz?
Answer   : It's not. Your $\omega_1$ mapping is not correct. Digital frequency $\omega \in [-\pi, \pi]$ in a filter response $H(e^{j\omega})$ is mapped to continuous frequency as follows :$$ f = \frac{\omega}{\pi} .\frac{f_s}{2}$$ $$\omega \in [-\pi, \pi] \rightarrow [-\frac{f_s}{2}, \frac{f_s}{2}]$$.

A Moving Average filter is a causal Low-Pass filter of which the cut-off frequency($\omega$) depends on $L_m$.

Intuition : It averages $L_m$ number of input samples, i.e. average of 1 current sample $x[n]$ and $L_m - 1$ past input samples. Averaging smoothens out variations from the input samples because the consecutive samples of output are changing very slowly. Why? Because, the variation in output sample can only be caused by the current input sample $x[n]$ which is small compared to sum of $L_m - 1$ past input samples. so, a large part of consecutive output samples remain approximately same. Hence, Moving average is a LPF and that is why useful to filter out noise with low computational complexity.

Now, since you have $|H(e^{j\omega})|$, so you can plot it against frequency and see where your desired band of 55Hz to 65Hz lies inside the Moving Average Filter response.
Figure 1 : Frequency response of the filter you have mentioned $L_m=30$, $F_s = 1440Hz$  and grey patch shows where the band 55Hz to 65Hz lies. It actually lies in the side lobe of MA response.

Figure 2 : Zoomed in Frequency response of the filter you have mentioned $L_m=30$, $F_s = 1440Hz$  and the patch shows where the band 55Hz to 65Hz lies.

You need to bring your interested band containing frequency from 55Hz to 65Hz inside the main Lobe and close to DC to get a less distorted 60Hz sinusoidal signal.
You can do 2 things to achieve your goal :

Decrease $L_m$ so that cut-off frequency of Moving Average(MA) filter moves away from DC and 60Hz sinusoid falls inside main lobe.

Increase sampling frequency $f_s$, so that input signal spectrum shrinks and as a result 60Hz comes inside the main and close to DC.

Option 1 : Decrease Filter length, $L_m = 10$, and $F_s = 1440Hz$

Option 2: Increase Sampling rate, $L_m = 30$, and $F_s = 30*1440Hz$

As you can see that the effect both methods is that 60Hz now falls inside Main Lobe and hence you will retain you r desired signal of 60Hz. But the 2 options are not same and $2^{nd}$ option of increasing sampling rate while keeping $L_m = 31$ will perform way better than decreasing $L_m$ to 10. Why?
Two main reasons:

Decreasing $L_m$ to a lower value will not have as good denoising effect as a higher $L_m$, since the averaging is done over a shorter length.

Increasing sampling rate and then low-pass filtering your spectrum will inherently increase SNR. You can think about how!


Which order to perform downsampling and filtering? I have EEG data recorded with 128Hz sampling rate.
As my goal is to reduce the amount of data (and maybe noise), I want to downsample the data to 64Hz (I am only interested in range 0.5 - 30Hz).

I would perform downsampling (64Hz) and use a bandpassfilter (0.5 - 30Hz).

Now my question is, what is the best order to perform those two steps and why?
Which impacts will these methods have to the quality of my data?

You need to filter first and then downsample. Otherwise, you will run into aliasing problems. I.e. frequencies that are above 30 Hz will create images within your frequencies of interest. You can consider the little script below to compare both methods:

Fs = 128.0
t = np.arange(0, 10, 1/Fs)
signal = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*50*t)
sigma2 = 0.5

rx = signal + np.sqrt(sigma2) * np.random.randn(len(signal))

downsampled = rx[::2]
(b1, a1) = scipy.signal.butter(6, 30/(Fs/4))
down_and_filtered = scipy.signal.lfilter(b1, a1, downsampled)

(b2, a2) = scipy.signal.butter(6, 30/(Fs/2))
filtered_and_down = scipy.signal.lfilter(b2, a2, rx)[::2]
t2 = t[::2]

plt.figure(figsize=(20,20))
plt.subplot(221)
plt.plot(t, rx)
plt.plot(t, signal)
plt.xlim((1, 1.3))
plt.title('Received signal')

plt.subplot(222)
plt.plot(t2, down_and_filtered)
plt.plot(t2, filtered_and_down)
plt.title('Compare of both methods in time domain')
plt.xlim((1,1.3))

f = np.linspace(-Fs/2, Fs/2, 4*len(t))
plt.subplot(223)
plt.plot(f, np.fft.fftshift(abs(np.fft.fft(rx, 4*len(t)))))
plt.title('RX spectrum. Note the two peaks, one is interference, having frequency over 30Hz')

f2 = np.linspace(-Fs/4, Fs/4, 4*len(t2))
plt.subplot(224)
plt.plot(f2, np.fft.fftshift(abs(np.fft.fft(down_and_filtered, 4*len(t2)))), '-o')
plt.plot(f2, np.fft.fftshift(abs(np.fft.fft(filtered_and_down, 4*len(t2)))))
plt.title('Spectrum of both alternatives. NOte the blue curve has a wrong component')

In this script you have a signal that is composed of two sines of frequency 10 and 50Hz plus some noise. In the FFT of the full signal, both spectral lines clearly occur. Now, what you want to have after your signal processing is only one spectral line at 10Hz (because the 50Hz signal should be filtered out). HOwever, if you downsample first, the 50Hz wave is mirrored to $14Hz=64Hz-50Hz$ and cant be filtered out subsequently. Hence, you need to do the filtering before and then downsample your signal.


Why do I keep seeing this equation for $E_b/N_0$ I keep seeing $$\frac{E_b}{N_0} = \frac{E_\text{s}}{\rho N_0}; ~~ \rho=\log_2{M}$$
But my own calculation is:
$$\frac{E_b}{N_0} = \frac{E_s}{N_0}\frac{1}{k}$$
Where $k = \frac{\rho n}{ms}$, where $\rho = msr$ is the spectral efficiency, $m$ is the modulation efficiency, $s$ is the percentage of carriers that carry data from the actual coded bitstream, $r$ is the coding rate $\frac{k}{n}$, $n = sNm$ is the total number of coded data bits in the symbol and $k$ is the total number of information bits in the symbol and $N$ is the number of subcarriers.
Wikipedia shows the following:

There are multiple issues with this

It uses $\rho$, which is spectral efficiency, and actually refers to it as spectral efficiency but then substitutes it with $\log_{2}M$ which actually is the modulation efficiency
It says that 'this is the energy per bit, not the energy per information bit', to follow up its use of modulation efficiency now instead, but there's a contradictory equation that is using the net bitrate i.e. the information rate $f_b$ in the section above
The contradictory equation, which is correct, produces my formulation, and the formulation that I've seen on a few sources i.e. this one: https://uk.mathworks.com/help/comm/ug/awgn-channel.html

$$\frac{E_b}{N_0} = \frac{\frac{P_C}{f_b}}{\frac{P_N}{B}}  = \frac{P_C}{P_N}\frac{B}{f_b}$$
and
$$\frac{E_s}{N_0} = \frac{\frac{P_C}{F}}{\frac{P_N}{B}} = \frac{P_C}{P_N}\frac{B}{F}$$
Which shows how $E_b/N_0$ differs in that it has a $f_b$ instead of an $F$ term. Therefore $F$ needs to be timesed by something that produces $f_b$, and that thing is $k$. $Fk$, the baud rate times the number of information bits in the symbol, is the information rate.

The equation I keep seeing across multiple sources as well as wikipedia is dividing the energy of a symbol by modulation efficiency rather than the number of bits in the symbol, which does not make semantic sense if the symbol has multiple subcarriers, and if it doesn't, this would only give the energy per gross bitrate transmission bit.

My guess as to what they have done wrong is calling $\rho$ $\log_{2}M$ or using $\log_{2}M$ in these equations at all, because AFAIK $E_b/N_0$ is with respect to the information rate. If they use actual spectral efficiency, then it is correct IF the symbol has one carrier. If it has multiple then they'd have to times the single carrier spectral efficiency $\rho = mr$ by $sN$ to get $\rho sN$ i.e. $k$

It's just a matter of finding out the assumptions behind the equation. Sometimes, unfortunately, those assumptions are not made explicit.
For the first equation you present, the one you have issues with, the assumptions are (off the top of my head):

Orthogonal, ideal sinc pulses are used as pulse shape.
The signals involved are strictly baseband, or, if passband, then the upconversion is suppressed-carrier, and the entire process is ideal.
The transmitted signal is narrowband; there are no multicarriers.
The channel is AWGN; there is no signal distortion, not even attenuation.
There is no coding and the transmitted bits are uncorrelated.
The received signal is matched-filtered. The matched filter is ideal (i.e. an infinite sinc pulse).

Under these assumptions, that equation is correct.
If you change those assumptions, then the relationship changes, as you found. I didn't check every one of you equations carefully but in general it seems like you are on the right track.


Impulse response (general form for linear systems as two-variables function $h(t,\tau)$) applied to Time-invariant systems If a system is linear, then the operator $S$ mapping input signals into output signals - i.e. $y(t)=S\{u(t)\}(t)$ - is the integral of the input weighted by the impulse response:
$$ y(t) = \int_{-\infty}^{+\infty} h(t,\tau) u(\tau) \mathrm{d}\tau $$
where $h(t,\tau):=S\{\delta(t-\tau)\}$, response of the system to an impulse in $\tau$, is a function of two variables: the second one $\tau$ refers to the time at which the input is applied while the first one $t$ refers to the time at which the output is observed.

If this system is also time-invariant, the system operator commutes with the time shift operator:
$$S\{u(t-T)\}(t) = S\{u(t)\}(t-T) $$
that is, time shifting the input, a time shifted version of the output is produced.
Linear time invariant systems have an impulse response function which depends upon a single variable, the difference between $t$ and $\tau$:
$$h(t,\tau)=\tilde h(t-\tau)$$ where the $\sim$ symbol is used just to underline that $h$ and $\tilde h$ cannot be the same, because they have a different number of arguments.
How can I show that the impulse response can be expressed as a function depending only on $t-\tau$ starting from the definition of time invariance?

If the response to $x(t)$ is given by

$$y(t)=\int_{-\infty}^{\infty}h(t,\tau)x(\tau)d\tau\tag{1}$$

then the response to $x(t-T)$ is

$$\tilde{y}_T(t)=\int_{-\infty}^{\infty}h(t,\tau)x(\tau-T)d\tau=\int_{-\infty}^{\infty}h(t,\tau+T)x(\tau)d\tau\tag{2}$$

If the system is time-invariant we require

$$\tilde{y}_T(t)\stackrel{!}{=}y(t-T)=\int_{-\infty}^{\infty}h(t-T,\tau)x(\tau)d\tau\tag{3}$$

For $(2)$ and $(3)$ to be equal for any $x(t)$ we must have

$$h(t,\tau+T)\stackrel{!}{=}h(t-T,\tau)\tag{4}$$

which is equivalent to

$$h(t+T,\tau+T)\stackrel{!}{=}h(t,\tau)\tag{5}$$

From $(5)$ it is clear that for a linear and time-invariant system the value of $h(t,\tau)$ only depends on the difference of its arguments, and not on their individual values. Consequently, the impulse response can be rewritten as a function of $t-\tau$:

$$h(t,\tau)=\tilde{h}(t-\tau)\tag{5}$$


Why do I see ringing in the output of a digital filter with a narrow transition band? I am doing some 'extreme' eq for spectral mangling type effects with audio. Im using brick wall filters, and very narrow band pass and reject filters (vst plugins), and I'd like to know if there's anything I can do about the pre/post 'ring' with the linear phase /minimal phase filters Im using.  Unfortunately I must use steep eq slopes . Im prepared to use minimum phase as it avoids pre-ring.

Specifically, Im wondering:

What exactly causes the oscilliations in the impulse response just after the input, in a a minimum phase filter?
Are these osccilations what causes the audible pre and post 'ringing' sound which is added to the passband with with steep slope filtering?
Are the oscilliations, and thus the ringing frequency always the same frequency, or does the ringing frequency depend in some way on the input signal?

Thanks very much for your expertise. I look forward to any responses.
Dale.

Your observations are an example of the Gibbs phenomenon. When you apply a filter with a very sharp transition band, you will observe oscillations in the filter output (or "ringing") near any sharp transitions in the input signal (e.g. boundaries of pulsed waveforms). The apparent "frequency" of the oscillations is dependent upon the bandwidth of the filter; as you increase the filter's cutoff frequency, the oscillations will become more localized in time (i.e. "higher in frequency"), but the peak overshoot does not change. The Wikipedia article linked above has a good explanation halfway through or so.


Very first values of y[n] in an IIR filter I use a simple biquad IIR filter  :

if (n >= 2)
   y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2]
                    - (a1/a0)*y[n-1] - (a2/a0)*y[n-2]
else
   y[n] = x[n]     // here it probably needs to be modified

(thanks to http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt)

It works very well.

My problem is that the filtered signal sometimes explodes in the first milliseconds, and then returns to normal after a few milliseconds.

How should I initiallize the first values of y[n] ? I have done this

y[n] = x[n]

for n=0,1, but probably the problem comes from here ?

The problem (it explodes in the first milliseconds) is important when I do 2, 3, or more passes of the same filter in cascade.

Thanks in advance.

Typically, you would initialize the filter state (which in your example includes $x[-1]$, $x[-2]$, $y[-1]$, and $y[-2]$) with zeros. This is equivalent to assuming that the filter is causal.


FFT for line detection I am trying use the FFT in a different way then most people ask about. I want to be able to take a picture of a graph with regular repeating vertical lines, and to process the image to determine how far in pixels the lines are apart on average. I have tried canny edge detection and hough line detection and I don't think I can optimize the images enough to accurately detect only the lines I am interested in.

So, my attempt is to scan 10 lines of the picture and to accumulate the pixel values into bins corresponding to the pixel column. What results when you graph it is a very nicely appearing waveform. When I perform a DFT or FFT on this, I can find a peak that I believe should be the frequency of the line repetition. (This may be a faulty assumption)

My question is, what does this number correspond to? i.e. I think I am confused with what my sampling rate would be because it is in pixels. I do think that this is a valid use of the FFT, but am falling right here at the point when i think I should be successful.

As an example. I created a picture that is 300 pixels in width. There are 1 pixel width lines drawn at exactly 30 pixel intervals. I found 2 peaks, one at 75 and one at 225 (which seem symmetrical) for the real component. (I do not believe that the imaginary component should play in??) I know the lines are 30 pixels apart. How does the 75 and or 225 relate?

I am trying really hard to get this, and I am grateful for any help you could recommend. At this point, I am giving up on edge detection, and want to try this approach.

Thank you in advance.

You might look at trying Autocorrelation for this. Here is an SO answer describing how to perform autocorrelation with Matlab using FFTs. This could be extended for two dimensions.

I implemented your test case in numpy as follows:

a = np.zeros(300)
a[::30] = 1
plt.acorr(a, maxlags=50)

This gives the following plot:

As you can see, the peaks pop up at +/- 30.


Is there a penis detection demo similar to face detection? Tutorial for face detection:

Is there a script / tutorial / demo for penis detection? These guys ran into some issues:

Fairly serious question, future of internet memes is at stake. Breast / nipple detection could be handy too.

Is there a script / tutorial / demo for penis detection?
[...]
Fairly serious quesion, future of internet memes is at stake

Yes, there is.
Common Pattern Recognition techniques will be able to spot one even with what would be considered today "traditional" approach (i.e. without "Deep Learning").
There already is a sub-category of the subject with people working on "Adult Content Detection". Here are a few relevant publications.
This is an example for detecting child pornography which contains an overview and references to other techniques along with what the authors are proposing.
There is an increasingly large amount of tutorials out there on how to use a number of popular libraries to build an "Image Classifier", some of them not even requiring writing code. So, getting to a "script" that does what you are after is not a huge amount of work.
Hope this helps


smooth noisy irregularly spaced data containing peaks I've got a set of scans of an object (human body) from different angles, which are being combined to reconstruct a 2D-representation. The raw measurements (blue plot below) contain a fair amount of spikes, which I reduce by using a median filter (result in green plot below). Afterwards I filter the result with a FIR lowpass filter to further remove noise (red figure below).

This produces at least visually appealing results, though I don't know if it is the best possible solution. Are there any pitfalls, when chaining median and linear filtering?

Another problem is, the samples are not evenly spaced. The distribution of samples seems to be somewhat $1/x$-ish, though I don't take this into account when smoothing the data.

Can this be ignored for simple smoothing? If yes, why? If no, how can I factor in the irregular spacing?

The idea of cascading a non-linear then a linear filter is sound in the presence of spikes.

Here are some observations which, if correct, could be exploited for better processing:

Data is circularly arranged, and in fact 2D
Spikes point outward,
The remaining noise seems bounded (the thick blue ring),
The angular sampling is dense, and could be affected by noise.

Some suggestions:

Choose a novel angular sampling $\delta \phi$ that could be used for several datasets. It could be parametrized by a modelling of histograms with a one-sided exponential law, like $a\exp(-\lambda t)$, $t\ge 0$.
Choose an angular width $\Delta \phi$, that will be used as the span of a 2D filter.

Their combination would define a series of overlapping windows on which you would process 2D point clouds.

Inside each angular region defined by $[n\delta \phi - \Delta \phi/2, n\delta \phi + \Delta \phi/2]$, you can for instance:

define a nonlinear filter to eliminate strong spikes: for instance a 2D median, or a 2D weighted median (A fast implementation of two-dimensional weighted median filters). If you use the outward direction of spikes, an order filter, or a rejection of the top $x$% of the radius could be used. The idea behind weighting is to penalize the points that are the farer away from the angular sampling. A bilateral filter for point clouds could be an option too.

Introducing a weighting in a median filter can bridge the gap between pure median and smoothing filters, and kill two birds with one stone: remove spike, reduce noise. As a result, this could suffice. If not, you can then:

perform a 2D linear filtering on the remaining points, in the shape of a center of mass allocation.

One advantage with resampling to a regular angular grid is potentially to allow easiest further processing, comparisons, etc. And using 2D filters, you naturally follow the circular symmetry of your data.


My vocoder sounds bad I'm writing a program to emulate a vocoder and it obviously works, but any speech that comes through is very difficult to understand. I'm going to explain exactly how my implementation works and then maybe somebody with experience with vocoders or DSP can help me.

128 samples of the voice signal go through a Fast Fourier Transform every audio frame. I then take the absolute value of the complex outputs to get the formant of the voice. I chose this rather than multiple band pass filters and envelope followers because it's faster just to do an FFT. I am 100% sure that this is all done correctly.

My carrier wave is run through several band pass filters centered around frequencies corresponding to the FFT output frequencies. To construct the output I multiply the output of these filters by the corresponding amplitudes in the voice formant obtained from the FFT and add all of these products.

There's a few things that I suspect may be causing issues.
The algorithm I used to implement the band pass filters is probably not the best. I feed the output of a low pass filter into a high pass filter. I got the algorithms for the filters under the "Algorithmic Implementation" section of the low pass and high pass filter articles on Wikipedia. I have a strong feeling that these simple filters may not be attenuating unwanted frequencies very well but I can't find a better way to make my band pass filter.
I also suspect that my carrier wave might not be very good. As I understand it, you need a carrier wave rich in lots of different frequency content (so obviously not a sine wave). I've listened to a lot of vocoders and the carrier wave usually sounds like it's got some dissonant notes in it. I've tried additive synthesis and FM synthesis but I just can't get a carrier wave that doesn't sound too "clean".

If somebody could give me a better band pass algorithm (if that's what I need), a good formula for a carrier wave, or point out some other likely problem I would be very grateful.

ALSO, "vocoder" is not an available tag. Maybe somebody with enough permissions could kindly add this tag? Also, there's a tag for signal-processing but not digital-signal-processing or DSP, I'm not sure if it would be appropriate to add one of these tags or not.

First, it seems to me that you are doing the analysis incorrectly. 128 linearly spaced channels is both too high (it won't capture a smoothed spectral envelope - vocoders usually have at most ~30 channels) and too small (it has a very poor resolution in the lowest frequencies and won't discriminate different formants). You really need to replace that by a bank of constant-Q band-pass filters, 2 or 3 bands / octave to get something that is coarse enough to capture the spectral envelope without picking individual frequency peaks, yet discriminates the frequencies of interest.

Your analysis and synthesis filters can be implemented with biquads. Running a dozen of those in parallel is not cost-prohibitive compared to a FFT.

As for the carrier wave, you can start with a sawtooth or narrow pulse - avoid waveforms like triangle, square, or sine whose spectrum have "gaps". Then stack a few of them to play a chord. You can thicken the carrier by adding slightly detuned version of each component or octaves. Note that in order to achieve the best sounding results, band-limited synthesis (for example with minBLEPs) has to be used to synthesize the square or sawtooth waves.


Phase Locked Loop Basics I am trying to better understand discrete-time PLLs. I drew a picture of a simple feedback system here:

I want to be able to describe and understand the system at v(t) and y(t). I found a fairly good reference with the following diagram:

I am still having a hard time and I think if I could see a signal go through this system either visually or computationally, I would be able to understand it better. So for example in my diagram, what would $y(t)$ and $v(t)$ be if $x(t)=5u(t)$ and $h(t)=6{\delta(t)}$?

You're probably better off skipping the $5u(t)$ idea and going straight to a periodic input, probably a sinusoid.  A phase locked loop is for signals that have a phase to lock onto, whereas $u(t)$ doesn't have a phase.  The Stanford reference looks pretty good.  Start by comparing that to your feedback system:  Note that the phase detector is a multiplication, where you appear to have an addition.  Also notice that the VCO isn't actually putting out $\hat \theta$, but $\hat s(t) = \cos\left(\omega_{lo}t + \hat\theta(t)\right)$.

So after the multiplication you end up with the sum of two new sinusoids, one at the sum of the frequencies of $\hat s(t)$ and $s(t)$, the other at their difference.  The loop filter is there to remove the higher frequency one, leaving just the difference, which will be $\sin\left(\theta(t) - \hat\theta(t)\right)$.

Now consider a few different situations.  First, the easy one - the input is running at exactly the same frequency as the local oscillator, so $s(t) = \cos(\omega_{lo}t)$.  Assume there is nothing on the input to the VCO to begin with, and so $\hat\theta(t) = 0$ too.  The output from the phase detector will be only the high frequency term, and that will get filtered out by the loop filter, leaving $e(t) = 0$.  The system is in equilibrium, everyone can go home and have a beer.

But now the frequency of the input signal increases a little.  $\theta(t)$ is now a linearly increasing value (if it were a constant, then the input signal would still have the same frequency, but just be delayed or advanced a bit).  Now the phase detector is also producing a low frequency term, $\sin(\theta(t))$ in addition to the high frequency one.  This is passed through the filter, and since the frequency only increased a little, $\theta(t)$ is small so $\sin(\theta(t)) \approx \theta(t)$.

This means that $e(t)$ is no longer zero, and so the phase of the VCO begins to increase; $\hat s(t)$ runs faster.  If it hasn't caught up to the frequency of $s(t)$, then the phase difference will continue to increase, and since that difference is being integrated by the VCO, the VCO frequency will increase faster.  Eventually the VCO will be running at the same frequency as $s(t)$.  However, note that it will be lagging in phase:  To run at the same frequency, $\hat\theta(t)$ needs to be increasing at the same rate as $\theta(t)$.  For $\hat\theta(t)$ to be increasing linearly, $e(t)$ needs to be a constant, and non-zero, which means $\theta(t) - \hat\theta(t)$ must be non-zero.

If the frequency of the input signal were to decrease, you would get a similar situation, except that now the VCO output would be leading the input signal rather than lagging.

In a real situation, the input signal frequency is usually changing all the time (otherwise you wouldn't need a PLL to track it), and all this is happening constantly.  Hopefully you can see from this description that if the input frequency changes too quickly then the system can be thrown off (remember the small angle approximation I used).

I believe the principles for a discrete-time PLL are similar.  Sampling rates with respect to the input frequency will need to be considered, as usual.

Hope that helps, I'm not sure if you are completely familiar with continuous PLLs and just looking for the differences when you go discrete, or the whole thing is a little fuzzy for you.  Since you've got continuous-time looking variables in your feedback system diagram I've assumed the whole PLL game is of interest to you.


Programming language for sound manipulation Having had little experience in PHP, I moved to C++ in-order to create a spectrogram, I wanted to create probably the simplest/easiest one for my first project, so easy infact I don't need reading/write to files or add settings to it.

But, as far as I have been learning C++ I was able to manipulate audio file, only to read and write, but the spectrogram works, by reading from a mic and this seems kind of hard to do. I have looked every link in Google, with the word spectogram and C++ in it, but I haven't gotten any useful page.

Now, if I am not mistaken, for the project I am working on, it may be good idea to use JAVA, because I read somewhere that JAVA is better than C++  1 in a way it manipulates sound, and 2 it is actually easier to learn than C++

So, I wanted to ask here, if JAVA is better for my project, which I will be doing in Windows7?

For the first few experiments I would recommend using a scripting language like Matlab or Python. They're much easier to understand and much quicker to write than "lower level programming languages" like C++.

Matlab has a signal processing toolbox and can read and write audio files, do windowing, FFTs etc. as well as a very simple playback mechanism. Basic 2D and 3D-plotting tools are also built right into the language itself.

For Python there is the SciPy library that does all kinds of scientific computing, including wavefile handling, windowing and FFT as well. Plus, there is the Matplotlib that does all of the plotting for you.

In C++, you would have to find a UI toolkit that works for your OS, write a drawing widget to plot stuff in and then fill it with data. Plus, obviously all the memory management; wich will probably become the biggest annoyance and time-sink if you're not experienced in doing it.

To implement live playback and recording I suggest using portaudio, a library written in C. There exists a Matlab .MEX utiliy called PlayRec and there are Python bindings called PyAudio.

From my experience, if you want to use portaudio I'd suggest using Python. Compiling the .mex utility is a huge pain while Python seemed to "just work".

If you manage to get it to work but still want to do it in C++ you can then port the codebase (just a couple of lines in both Matlab or Python) while still using portaudio for the audio input/output.

For a superficial comparison of the different code-complexities you can take a look at some work I wrote last year. Each implementation does the same thing, just in a different language:

FTM filter written in C/CUDA (plus some benchmarking stuff)
FTM filter written in Matlab (plus some plotting stuff)
FTM filter written in Python (just the filter, no extras)


How to Calculate the Hann Low Pass Filter / Kernel for Image Processing I'm currently studying for an exam in image processing and stumbled upon a exercise which i could not answer and my professor will not be available before the exam anymore. The exercise goes as follows:

  Create a Hann-lowpass kernel of size $3\times 3$ and calculate the result of the convolution at $(2,2)$ in the image.

The given solution for the filter kernel is shown below, but I do not understand how he got there:

$$
A = \frac{1}{4}
\begin{pmatrix}
\frac{1}{4} & \frac{1}{2} & \frac{1}{4}  \\
\frac{1}{2} & 1 & \frac{1}{2}  \\
\frac{1}{4} & \frac{1}{2} & \frac{1}{4}  \\
\end{pmatrix}
$$

The second part of the question is a simple convolution with the image which I'm able to do. The factor in front of the kernel is given by the constraint $ \sum A = 1$ to achieve a filter gain of one.
I've done literature research but none of the books in the library do explain how to calculate such a kernel. The professor gives a "hint" in his presentation pointing to the generalized cosine window with $A = 0.5$, $B = 0.5$ and $C = 0$:

$$
w_k = A - B\cdot \cos\left(2\pi\frac{k}{K-1}\right) + C\cdot \cos\left(4\pi\frac{k}{K-1}\right)
$$

Yet I cannot derive the final kernel from the given formulas. Can anyone give me a step by step solution for the problem?


What Is the Definition of Linear Predictive Coefficients When the Optimal Value Aren't Unique? Linear predictive coefficients of signal y are defined as the best $k$ coefficients $a_i, i = 1, \ldots, k$, that will approximate $y_n$ by $-\sum_{i=1}^k{a_iy_{n-i}}$. (Best approximation is that which minimizes the sum of squared errors.)

How are LPCs defined when this problem has multiple solutions, e.g. when y$$ is constant zero?

Constant zero sequence is the only finite sequence for which $a_i$ are not uniquely defined.

Indeed, $a_i$ are a solution to the equation
$$
\begin{bmatrix}
R_0 && R_1 && \ldots && R_{k-1} \\
R_1 && R_0 && \ldots && R_{k-2} \\
\vdots && \vdots && \ddots && \vdots \\
R_{k-1} && R_{k-2} && \ldots && R_0
\end{bmatrix}
\begin{bmatrix}
a_1 \\
a_2 \\
\vdots \\
a_k
\end{bmatrix} = -
\begin{bmatrix}
R_0 \\
R_1 \\
\vdots \\
R_{k-1}
\end{bmatrix}
$$
To prove their uniqueness, it is sufficient to prove that the autocorrelation matrix R is nondegenerate.

If $y_0, y_1, \ldots, y_n$ is the initial sequence, $y_0 \ne 0$, and
$$
Y = \begin{bmatrix}
y_0 && \ldots && y_n && 0 && \ldots && 0 \\
0 && y_0 && \ldots && y_n && \ldots && 0 \\
\vdots && \ddots && \ddots && \vdots && \ddots && \vdots \\
0 && \ldots && 0 && y_0 && \ldots && y_n
\end{bmatrix},
$$
then
$$
R = YY^T,
$$
so R is the Gram matrix of the rows of Y. A Gram matrix is nondegenerate iff the vectors are linearly independent, which they obviously are.

To deal with the constant zeros problem, I should probably just trim silence at the beginning and end of the audio.