Frequency synchronization in OFDM Could you please help me understanding how the null sub-carriers can be used to estimate the carrier frequency offset and thus performing frequency Synchronization at the receiver?

That depends, but usually rough frequency correction in an OFDM receiver is done by simply looking at which carriers have significant energy, which don't; you can do that, for example, by comparing the sum of all non-null carriers' energies to that of all null-carrier energies, and shifting that around until the ratio maximizes. Pretty simple idea!


How many bps would be required to reproduce the words in speech recognizably? In Robert Gallager's lecture notes for 6.450 Principles of Digital Communications I, Exercise 2.1 asks an interesting conceptual question:

  Voice waveforms could be converted to binary data by sampling at 8000 times per second and quantizing to 8 bits per sample, yielding 64kb/s. [...] Modern speech coders can yield telephone-quality speech at 6-16 kb/s. If your objective were simply to reproduce the words in speech recognizably without concern for speaker recognition, intonation, etc., make an estimate of how many kb/s would be required. [...] (Note: There is clearly no “correct answer” here; the question is too vague for that. The point of the question is to get used to questioning objectives and approaches.)

My first instinct was to look up number of words used today in English language (~170,000), calculate how many bits would be required to brutely encode each word (~18), and look up, on average, how fast we speak in terms of words per second (~2), and come up with 36 bits per second as my crude answer. I'm not sure if not concerning with the speech waveform itself, and all the signal processing that would entail, is a reasonable way of thinking about the problem. I'm here for feedback and other ways of thinking about this problem.

(This is not a homework assignment. I was studying the notes myself and I'm curious.)

This all depends how you interpret the requirements and what is "good enough".

For example, a really bandwidth efficient  method would to simply encode with a speech to text system, send the text and re-encode with a text to speech system. That would be about 0.075 kb/s (150 words per minute, 6 letters per word, 5 bits per letter).

That would make a pretty terrible speech coding system but whether that meets the "rules" or not is more a semantics than a technical question.


What is the easiest learning curve to develop a Mark Reader software? I am a C# programmer and have a project to create a Mark Reader software. Bubble sheets that have been filled in by students are scanned using a general scanner (as opposed to a special scanner for bubble sheet). The resulting images must be processed using the Mark Reader software in question.

I have a bit knowledge in DSP. Shortly speaking, what is the easiest learning curve to develop the software? Book recommendation is also welcome.

EMGU is a .NET wrapper for OpenCV which is a powerful library for computer vision applications.

You could use something like blob detection to set up the layout of your grid of answers in space, check the pencil shading against a threshold (to prevent erasures from being detected as answers), check that there's one answer per question, and that it's in the right spot.

Get some "toy" problems going with the library first (detect a circle or a square in a field) and build up your implementation from there.


Properties of Discrete-time Sinusoidal Signal I was reading in a book that  below stated are the Properties of Discrete-time Sinusoidal Signal:

Property 1. The discrete time sinusoids whose frequency are separated by an integer multiple of 2$\pi$ are identical.

Property 2. The frequency of oscillation of discrete time sinusoids sequence increases  as $\omega$ increases from 0 to $\pi$. If $\omega$ is increased from $\pi$ to 2$\pi$ then frequency of oscillation decreases.

I was able to understand the mathematical property implementation of above two property but have some basic questions for clarification:

Is it that the above properties are applicable for all Discrete-time periodic Signals?
As per property 1 - if two sinusoids signal has frequency separated by integer multiple of 2$\pi$ - then the signal having the greater frequency will have more cycles than other. Then how they both can be same?
As per property 2 - what is meant by "frequency of oscillation"? How can a frequency of a signal vary (increase) from 0 to $\pi$ and (decrease) from $\pi$ to 2$\pi$ - isn't it should be same from 0 to $2\pi$?

Any discrete-time periodic signal with period $N$ can be written in terms of its DFT coefficients $X_k$:

$$x(n)=\frac{1}{N}\sum_{k=0}^{N-1}X_ke^{j(2\pi/N) nk}$$

So whatever is true for a single sinusoid (or complex exponential) remains true for a general periodic signal.

Two discrete-time sinusoids with a frequency separation of $2\pi k$ are identical because the sine function is periodic with period $2\pi$:

$$x(n)=\sin(\omega n)\\
y(n)=\sin[(\omega+2\pi k)n]=\sin(\omega n +2\pi k n)=\sin(\omega n)=x(n)$$

Finally, the maximum frequency for a discrete-time sinusoid (and consequently for any discrete-time signal) is $\omega=\pi$:

$$\cos(\pi n)=(-1)^n$$

If the frequency is further increased you will observe aliasing:

$$\cos[(\pi+\Delta\omega)n]=\cos[(\pi+\Delta\omega)n-2\pi n]=\cos[(\Delta\omega-\pi)n]=\\=\cos[(\pi-\Delta\omega)n],\quad 0\le\Delta\omega\le\pi$$

So from the maximum frequency of $\pi$ you get a decrease in frequency down to frequency $0$ (for $\Delta\omega=\pi$).


Automatic landmarks on images of hands Problem statement
I am starting a project soon where I need to create a method to automatically create landmarks on a set of images of human hands. I do not have a sample image but the images are taken with the same background and stable lightning conditions. The images probably look similar to this one, (but should be the other side of the hand).

My colleagues have defined some landmarks, but I would like to first just detect the tips of the fingers and the inner most part between the fingers. I was looking into potential solutions and active shape or active appearance models look like a good idea, although they might also be an overkill. I was also thinking about just creating a template of a hand with landmarks, segmenting the hand with foreground/background segmentation (probably thresholding, PCA or even MRFs if the others fail a lot), and then do an affine + non-rigid registration to transfer the landmarks.
What I am asking is:

Are the methods that I mention for solving this task appropriate?
Is there any solution out there that already solves this problem?
Referrals to code to do some of this is appreciated.
Do you have other better suited ideas to solve this problem?

Note
I am going to play around with the images in this dataset, which already have lots of annotated landmarks.
I'm new
Yes, this is my first post, so if you have ideas for other tags to include let me know. Also if you need further clarifications, don't hesitate to ask!

There are two simple approaches.

Find the hand region with thresholding and then use something like convexity defects. Even though you will have many spurious landmarks at the end, you could easily cluster them out, as your scene is well controlled. There are also couple of youtube videos, if you like to get an idea in advance, [1], [2].
Fit a hand model to the image. There are many works in this field, and it's an active academic area. Some are, [1], [2], [3]. Once the model is fit, you could extract any landmark you desire, since the correspondence to the canonical hand frame is determined.

The pdf you provided also gives a simple method for fitting a hand region.


Transform lowpass filter to highpass filter, what happened to unit sample response and difference equation? I know that multiple a lowpass filter signal by $(-1)^n$ can transform it to a highpass filter.
$$
s_{hp}[n]=(-1)^n s_{lp}[n]
$$
But exactly what happened to the unit sample response and difference equation? Assume that the original lowpass filter signal has a difference equation. Can someone explain to me?? Thanks!!!

I know that multiple a lowpass filter signal by (−1)n can transform it to a highpass filter.

Depends a bit what you mean by that. If you assume that this will turn 1 kHz lowpass signal into a 1kHz high pass signal, your assumption is plain wrong. What happens here is that SHIFT your signal frequency spectrum by half the sample rate, $f_s$. Due to the periodicity of the spectrum, this looks like the spectrum has been mirrored at $fs_/4$ (at least for real valued signals), i.e. DC become Nyquist, Nyquist becomes DC.
So this operation flips the entire signal spectrum, regardless of whether it's been high-pass, low-pass, pink, white, etc.
Loosely speaking, if you have signal that's been low passed and you flip it, the spectrum may look a bit like that of a high passed signal, but it's completely different operation than filtering the original signal with a complimentary high pass.

But exactly what happened to the unit sample response and difference equation? Assume that the original lowpass filter signal has a difference equation.

This is NOT an LTI operation, so difference equations and unit sample response do not apply here. Specifically, it's time variant.


Passing a sampled signal through a filter I was wondering why is it wrong to use a band-pass filter on a sampled signal?

If the signal we want to sample has frequencies up to fmax, we sample it with frequency fs = 2fmax (so that Nyquist theorem condition is fulfilled), and the spectrum of the sampled signal will have scaled "copies" of the original spectrum centered around frequencies k * fs (k = +-1, +-2 ...).

What would happen if we filtered that sampled signal with a band-pass filter [fs-fmax, fs+fmax], and why is it wrong?
Wouldn't we end up with the signal we started with?

Not sure what you mean by "wrong".

The Nyquist criteria simply requires you to have "two samples per Hz of bandwidth". It doesn't have to be $[-f_{max},-f_{max}]$, it can be any frequency range that includes at least $2 \cdot f_{max}$ of bandwidth. However, for real signals you need to figure out what to do with the negative frequencies.

For more info google "bandpass sampling".


Meaning of a formula about transient frequency bins I'm trying to understand the following formula in this paper:

Sadly, the paper does not define most of the symbols it is using, so my guesses so far are:

$T_{tr}$ refers to the sampling interval

i.e. $\frac{1}{44100} s$

$\phi$ refers to phase

specifically: $\phi(k_{tr}, (n-1)h)$ refers to the angle of the complex vector representing frequency $k_{tr}$ for the FFT from time point $(n-1)h$ to $nh$

so basically this is saying: out of all FFT frequency bins, only include those frequencies $k_{tr}$ where above equation holds true

I tried googling around to see if there's a standard definition of "transient frequency bins" (like there is for FFT frequency bins) this is referencing, but no success. I also don't understand

Why would one compare radians $\phi$ with a unit of time $T_{tr}$?

Is this maybe a stand-in for 0 to account for sampling errors?

Why isn't $T_{tr}$ just called $T$ as it seems to be a property of the signal and unrelated to the transient bins?

(If this is not the sampling interval, then how would one obtain $T_{tr}$?)

How does this formula define $K_{tr}$?

without reading the paper,

Eqn 4 has the form of a central difference form of the second derivative. It is a numerical second derivative compared to a threshold. It differs from the standard form in that it centers on $n-1$.

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

Whatever $\phi$ is, it is the positive change compared to a threshold.


If I build a low pass gaussian filter in the frequency domain should I also normalize it with the sum of the filter coefficients? I saw two examples of applying a Gaussian low pass filter to a Matlab image: one- when the filter was built in the space domain and was normalized with the sum of the filter coefficients, and another- where the filter was built in the frequency domain and no normalization was made.

Is it consistent that no normalization is needed in the frequency domain and why is it different in this respect from the spacial domain?

It is consistent to apply a normalization depending on what you really want to measure. Either on Fourier transforms, or low-pass linear filters, one customary sees three types of normalizations:

No normalization,
Normalization for a constant input,
Normalization in energy.

To be brief,  for a 3-point uniform coefficient filter, that is a low-pass filter, the impulse response $h$ can be:

No normalization: $h={1,1,1}$,
Normalization for constant output (mean-preserving): $h={\frac{1}{3},\frac{1}{3},\frac{1}{3}}$,
Normalization in energy: $h={\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}}}$.

Generally, the normalization depends on the invariance one wants to preserve. For low-pass filters, one often preserves mean (i.e. a constant image remains constant after filtering) or energy. For high-pass, some preserves energy, or estimated slope (or higher-order derivatives). Very often, the sum of coefficients of a high-pass filter is zero (a constant signal is zeroe out), thus it cannot be used for normalization. Yet, many other choices are possible.

So if you want to preserve the average value of the data, in time you normalize so that the coefficients sum to one. In frequency, you can remember that the amplitude of the $0$th frequency obtained from the Fourier transform is proportional to the sum of coefficients in time ("DC Term" in An Intuitive Explanation of Fourier Theory). So if it is non-zero, you can divise your whole Fourier spectrum by this value to get the same normalization.

The reason why people sometimes do not normalize the intensities with image filtering is that these values are often displayed and compared on a full range in the colormap (in Matlab, with function imagesc). Consequently, a scalar factor on the image is not considered harmful.


How is Weighted Thresholded Histogram Equalization different from Gamma Correction Image enhancement In case of Weighted Thresholded Histogram Equalization we use:-

calculate, the Weighted Probablity Density Function for Image using,

$$P_{wt}(k) = \left( \frac{P(k)-P_l}{P_u-P_l} \right)^r P_u$$

Where $P_u$ is the highest probability of a intensity in a Image.
      $P_l$ is the lowest probability of a intensity in a Image.
      $P(k)$ is the Probability Density Function for the Image,i.e. the number of pixels having the intensity $k$ in a image.

Then perform Histogram Equalization on the image using,
calculating Cumulative Distribution Function

$$C_{wt}(k) = \sum_{m=0}^k P_{wt}(m)$$
$$Im(i,j) = W_{out} C_{wt}\left(F(i,j)\right)$$

For $r<1$ the Power Law Transformation function will give a higher weight to the low probabilities in the PDF than the high probabilities.

How is this procedure conceptually different from the Gamma Correction technique, explain if similar how much similarity is there or if something is different, Explain where does it deviate from Gamma Correction method for Single Image Enhancement.

Does this Procedure Passes for Being Called as Gamma Correction in some Aspect OR Do some Changes are Required.

This Method can very well pass for being Gamma Correction Technique as We are Increasing the Weight of the Less probable Intensities in the Image  by the Exponential Factor r So that the Less Clear Details may be Visible upon the increase in their Weight thus Providing a Enhanced Image , On the Similiar Note the Gamma Correction Function Does this by using Power Law Transformation.


Daubechies Wavelets in mulitresolutional analysis decomposition I have an understanding problem with Daubechies wavelets.

When I use a multiresolutional analysis, I want to approximate the given input Signal $f\in L^2(\mathbb{R})$ on the subspaces $V_i$. By calculating the coefficients I split the signal in the detail part and the rest. Now I am wondering, how the Daubechies wavelet is fitting in this picture. The Daubechies wavelet is defined by coefficients, but in my opinion, these are only usable if I try to perform a wavelet decomposition on input data in $l^2(\mathbb{Z})$.

To specify my question, I have no idea how to calculate the first inner products
$$c_i = \langle f, \varphi_{i,j} \rangle, \text{ and } d_i = \langle f, \psi_{i,j} \rangle $$
if $\varphi$ and $\psi$ are Daubechies scaling function and wavelet, since they have no expression as functions in $L^2(\mathbb{R})$.

Thanks in advance

Matthias

As I mentioned in my comment, the reason seems to be the so-called "wavelet crime" which is taking the sampling data as the inner product. But I found a question, which is similar to mine

Wavelet transform: How to compute the initial coefficients when only samples are available?*

This question was answered with a lot of papers, many thanks to Laurent Duval. In one of the paper

On the Initialization of the Discrete
Wavelet Transform Algorithm

there is a derivation to calculate the initial inner products. One can approximate the inner products as
$$
c_k = \int_\mathbb{R} f(t)\varphi(t-k)dt \approx \varphi(-k)
$$
And one can read the following remark:

"A very good approximation of the time-shape of $\varphi(t)$ can
be obtained as the impulse response of the filter h convolved
with itself a large number of times. One then simply needs
to collect the samples ${\varphi(p) = \varphi(t = p)}_{p\in\mathbb{Z}}$ to perform the
approximate initialization"

Therefore it is sufficient to know the filter coefficients as given in Wikipedia.


How good is the Leech lattice for AWGN coding or MSE quantization of Gaussian sources? Anyone know where to find graphs showing how well it works vs other schemes? thanks

Tables 7.1 and 7.2 from Lattice Coding of Signals and Networks by Zamir are essentially this.

Also the figures in this paper: https://arxiv.org/pdf/1103.0171.pdf


PSD - Is this formula correct? Trying to understand whether or not I am calculating the PSD correctly. Below is the equation that I am using:

$$
PSD_{rms}(f_{m} = m \cdot f_{res}) = \frac{PS_{rms}
}{ENBW} = \frac{2 \cdot \left | y_{m} \right |^{2}}{f_{s}\cdot S_{2}} m = 0....N/2
$$

This equation is applied using the following:

1) Calculate the squared magnitude of the FFT bin:

$ mag = (re*re+im*im) $

2) For each of these magnitude, I calculate the following:

$ PSD(0...NFFT/2) = 2 \cdot (mag[0... NFFT/2] / Fs \cdot S2 $

So it looks like the following in code:

PSD[i] = 2 * (result[i]) / (12000 * S2);

I get the following result:

This looks like the result is being scaled correctly. However, the removal of the noise...

I have been told to use a band-pass filter, but, looking at the source code of Matlab in their spectrogram algorithm they do not pass the results through any kind of filtering algorithm.

My question is: Should the output be this, and, if so, should I therefore apply a filtering algorithm even though matlab does not?

Look at matlab result:

My question is: Should the output be this, and, if so, should I therefore apply a filtering algorithm even though matlab does not?

Without seeing your code I can't determine what's the reason for a difference in the two plots. There's multiple things that come into play like window type, overlap, and FFT size.

I can't say for sure if you should filter your data. It depends on the application. If you're just going to plot it, it may not be necessary. If the out of band energy or noise is really high and leaking into the signal band, then maybe. If this is going to be the input to some other process and you think the out of band energy is degrading some results, then sure, go ahead and filter. Without knowing your application it's tough to say what's right or wrong.

Anyway, in order to filter your signal to remove out of band noise, there's two approaches, but the easiest would be to do it in the time-domain.

Filter the data in the time-domain with either an FIR or IIR filter. Have a look at functions fir1 and fir2 in MATLAB. Once you design your filter, use freqz to quickly plot the frequency response to make sure it's correct. Once you have your filter 'b', compute:

% design filter 'b' using fir1 or fir2
y = filter(b, 1, x); % Filter input x with filter b
% compute PSD of y

Depending on your filter you should see suppression of energy in the stop bands.


Extracting the talking (lyrics) from an audio in python I am new in signal processing and I want to extract the talking (lyrics) of a person from a sound so I can analyze it ;another application would be if that person is talking and there are many sounds behind him the applications will detect his sound only.

I have read many articles about single processing like (analog to numeric,FFT,Gaussian, FTD , Correlation  ..... ) but I can't combine those all information together because i am new in this domain like i said before.

please  someone help me I am totally lost
---> sorry for my English <----

What you are trying to do might be extremely difficult depending on numerous factors. I have read about a technology called blind source de-convolution which tries to separate one audio source from others. In the end rather than recovering every audio source you will most likely recover only the one you are interested in.

I would suggest you start with this survey article which lists many of the attempts that have been done to achieve blind source de-convolution.

If what your are trying to do doesn't have to be too precise, as a very simple start you might just want to filter every frequency component above 4-5 kHz because the main energy of human voice is below that. However if it is a song with a strong bass, for instance hip-hop, pop, drums, etc, this won't help you clean that


Definition of energy of a signal I was wondering why the energy of a signal is defined as

$$E=\int_{-\infty}^{\infty}|x(t)|^2dt$$

I mean, I know that it is a definition, but is there a physical reason for it being defined like that?

As we know $E = Pt$ where $E$ is the energy and $t$ is the time. In the continuous domain $t$ is never ending and its beginning isn't known either.

Now, signals in the continuous domain are voltages and as you know from Ohm's Law $V = IR$ and $P = I^2R$.

This basically tells you that for a given resistance $R$, the energy of a signal over a time interval is given by the integration of the power. So,

$$ E = \int\limits_{t_1}^{t_2} P(t) ~dt$$

which in this case becomes

\begin{align}
E &= \int\limits_{-\infty}^{\infty} I^2(t)R ~dt\\
&= R\int\limits_{-\infty}^{\infty} I^2(t) ~dt\\
&= \int\limits_{-\infty}^{\infty} I^2(t) ~dt~~\text{(For $R = 1$)}
\end{align}

Usually in SP world we like to represent the signal as $x(t)$ giving us

$$E = \int\limits_{-\infty}^{\infty} \left|x(t)\right|^2 ~dt$$

Hope this helps!


Understanding the twiddle factors So I've been looking at this butterfly diagram to try to understand it better:

And I am trying to get a good understanding of the twiddle factors.
The definition is given as:
FFT Twiddle Factor: ${e^{i2{\pi}k/N}}$
and
IFFT Twiddle Factor: ${e^{-i2{\pi}k/N}}$
So k is the index number of the iteration thus $k=0,1...N$ but its $N$ that I am unsure of.
From the image is the first stage N = 8 (since there are 8 butterflies) or is N = 2 since each butterfly only spans two elements? Or is N always 16 every pass?
My guess currently is:
Stage 0 : N = 8
Stage 1 : N = 4
Stage 2 : N = 2
Stage 3 : N = 1

Am i correct in this assumption? If not i hope some one clarifying my misunderstanding here.
Thanks

You got it mostly backwards, but otherwise OK :)
This can be pretty directly answered by writing down what the Cooley-Tukey FFT's twiddle factor $W_N$ is:
$$W_N=e^{-i\frac{2\pi}{N}}$$
and that's it. For example, in your picture, in Stage 0, you're multiplying with $W_2^0=\left(e^{-i\frac{2\pi}{N}}\right)^0=e^{-i\frac{2\pi0}{N}}=e^0=1$.

So k is the index number of the iteration

No, $k$ is the exponent of your twiddle factors. It goes from 0 to half the size of your transform (in each stage).

From the image is the first stage N = 8 (since there are 8 butterflies) or is N = 2

Stage 0 has subtransforms of size N=2, as you can see. (Also, just write down the 2-DFT's formula and compare.)


Why do we always characterize a LTI system by its impulse response? Why do we always characterize a LTI system by its impulse response and not by another response, like the step response? What does the impulse response have that is so special?

Let's try a simple example.  A discrete-time LTI system has unit pulse response
$$h[0] = 1,\quad h[1] = 1,\quad h[n] = 0 ~~\text{for all}~ n \geq 2.
\tag{1}$$
The response $y$ of this system to a signal $x$ given by
$$x[0] = 1,\quad  x[1] = -1,\quad x[2] = 1,\quad x[n] = 0
~~\text{for all}~ n \geq 3 \tag{2}$$ is readily calculated
via the convolution formula,
and we get
$$y[0] = 1, \quad y[1] = 0, \quad y[2] = 0, \quad y[3] = 1,\quad y[n] = 0
~~\text{for all}~ n \geq  4.\tag{3}$$

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

The unit-step response $g$ of the same LTI system is easily found to be
$$g[0] = 1,\quad g[n] = 2 ~~\text{for all}~ n \geq 1.\tag{4}$$

So, suppose that you are given only the step response $g$ in $(4)$
but not the pulse response in $(1)$.  How will you compute the response
of this LTI system to input $x$?  It is known that the answer is $(3)$,
but you are asked to figure out how to get from $(2)$ and $(4)$ to $(3)$.


Given a real signal, increase modulation frequency Lets say I receive the signal

$$x(t)\cos(\omega_c t)$$

and I want to change the the carrier frequency $x(t)\cos(1.001\omega_c t)$

How can this be done? If I modulated it with a complex exponential, increasing or decreasing the modulation frequency would have been simple, however I could not figure out a simple way of doing this.

Is it even possible to do this without filtering?

You can simply multiply by another cosine: $x(t)\cos(2\pi f_ct)\cos(2\pi 0.001f_ct)=x(t)\cos(2\pi0.999f_ct)+x(t)\cos(2\pi1.001f_ct)$ (ignoring scale factors). If $x(t)$ is very narrow in frequency, so that the spectra of the two terms above don't overlap, you can use a high-pass filter to do what you want.

If $x(t)$ is not narrowband enough, you can shift $x(t)\cos(2\pi f_ct)$ two times. With adequate filtering, you will obtain the desired shift.

Simply multiplying by a complex exponential is not going to work, since it will shift the entire spectrum of $x(t)$ up or down, and the resulting signal will be complex. You want to shift the positive frequencies of $x(t)$ "to the right", and its negative frequencies "to the left".

However, you can use the complex exponential in the following way: Filter out all the negative frequencies of the modulated signal (using for instance the Hilbert transform), and up-convert them using a complex exponential. Then, filter out all the positive frequencies, and down-convert the result. Then, add the down-converted and up-converted signals.

So, your options are:

Use one or two cosine signals and a few filters (all real).
Use two complex filters and two complex exponentials.


Shift in Time Domain After DFT Based Convolution To my understanding, multiplying a signal in the frequency-domain is equal to a convolution in the time-domain.

I wrote a small python program, but i always end up with a shift in the time domain. where does this come from? An hints how to get the correct convolution, where the convolved output signal is in place corresponding to the input signals?

import numpy as np
from matplotlib import pyplot as plt

# make simple rectangular signals
r1 = np.zeros(100, dtype=float)
r2 = np.zeros(100, dtype=float)
r1[1:20] = 1.
r2[20:40] = 1.

# fftransform the signals
fft_r1 = np.fft.rfft(r1)
fft_r2 = np.fft.rfft(r2)

# multiply the frequence-spectra of the signals
fft_p = np.multiply(fft_r1, fft_r2)
# ifftransform back to time-domain
conv = np.fft.irfft(fft_p)

f, axarr = plt.subplots(3, 3)
f.tight_layout()
axarr[0, 0].plot(r1)
axarr[0, 0].set_title("r1")
axarr[0, 1].plot(np.abs(fft_r1), '.')
axarr[0, 1].set_title("fft_r1")

axarr[1, 0].plot(r2)
axarr[1, 0].set_title("r2")
axarr[1, 1].plot(np.abs(fft_r2), '.')
axarr[1, 1].set_title("fft_r2")

axarr[2, 0].plot(np.abs(fft_p), '.')
axarr[2, 0].set_title("fft product")

axarr[2, 1].plot(conv)
axarr[2, 1].set_title("conv (ifft of fft_product)")
axarr[2, 2].plot(np.convolve(r1, r2, 'same'))
axarr[2, 2].set_title("conv, np")

Here is an image of the plots that the above code generates:

Note that multiplying the FFTs and taking an IFFT is equivalent to cyclic/circular convolution. This is not the same as linear convolution. However, your signals have sufficient zero-padding for the two types of convolution to be equivalent. The correct result of the linear (and cyclic) convolution is what you obtained via FFT and IFFT, a triangle between indices $20$ and $60$, as shown in the figure. This result is easy to predict using graphical convolution, as explained in this animation (from Wikipedia).

The strange result of np.convolve with the option same can be explained as follows: the result is the center part of length $100$ of the linear convolution of the two signals. Since the linear convolution is a signal of length $199$ ($=N_1+N_2-1$) with $99$ trailing zeros, the same option cuts out the center part between indices $50$ and $149$. So the result is the right half of the cyclic convolution with $49$ zeros added at the end, as shown in the bottom right figure.


What is the rate of expansion of an image as it aproaches the human eye? My intuition about  the scaling (x,y) of the 2d representation of the 3d object generated by our native visual processing system tells me it accelerates as it approaches the visual sensor (our eye).

Is this true?  What is the rate of expansion of objects within our visual field?

The magnification of an object at distance $g$ from a thin lens with focal length $f$ is given as

$$\beta=\frac{f}{g-f}$$

which is generally negative for objects at distances greater than the focal length. This means the image is upside down.

The derivative of the magnification with respect to the object distance is

$$\frac{\partial \beta}{\partial g}=-\frac{f}{(g-f)^2}$$

So for constant $f$ the magnification diverges like $x^{-2}$ where $x=(g-f)$ is the difference between focal length and object distance.

So yes, the magnification grows in an accelerated way and even approaches infinity.

However, the relationship is only strictly true if the image at the image plane is sharp and the aperture is unlimited. So for your eye, at some point the image will become blurry and the magnification will grow in an asymptotic way that is much harder to describe.


Detecting characters in high contrast/oversaturated regions I'm currently tagging insects with reflective tags and filming them in night-mode on my camera which has an infrared light. I have some markings on the tags which I need to catch on film. Unfortunately, when I film about 2 metres back, the the contrast between the reflective tag and the marking is too great, resulting in the white reflection being the only thing visible. Does anyone have any ideas for how to fix this problem? For instance by using filters?

Use your camera's RAW mode, if it has one, and underexpose. If these concepts are not familiar to you, head over to https://photo.stackexchange.com/

Is the first picture a close-up of the second one, because the letter K on the tag is quite legible? If not, you may also want to use a higher resolution camera, and/or not shoot from such a distance.

If you think about it, you are asking us to robustly encode roughly four bits (on letter of the English alphabet), judging by the first picture. You could simply map the binary representation of the letter to five boxes; four in a square and one on the side, so you don't get any confusion due to symmetry. (If you laid them side by side you wouldn't know which side is which.) If you need to encode multiple letters, you could try using QR or a barcode.

Your first order of business has to be getting the exposure under control. I'd investigate the possibility of using aftermarket firmware if the manufacturer's doesn't give you enough control.


Why 1.01 equals -0.75 in Q2 I'm trying to learn Q notation but I don't understand why 1.01 in the Q2 number format equals -0.75. How is that calculation done? I can understand that 1.1 equals -0.5 because the leading 1 means the number is negative and the first digit after the binary point is 1 which means 1*2^(-1) which is 0.5 and therefore the number is -0.5. Similarly, 1.01 should equal -0.25 and not -0.75.

It is 1.11 which is -0.75 or what is my mistake?


Updating FFT algorithm accordingly when upsampling/resampling I have some experimental data from a shaker (vibration exciter) running at sine $10\textrm{ Hz}$ for test purposes (recorded using an ADXL345 accelerometer with $f_s=512\textrm{ Hz}$. I would like to upsample these data for varies ratios, e.g. a factor of $\times 10$.
How do I update the FFT in MATLAB accordingly? As of now, upsampling appears to shift results in the $x$-direction as shown in the figure below (black is the original signal, blue is the resampled signal with cutoff at $1/2\times f_s$):

Peak at approximately 80 Hz (black) appears at approximately 40 Hz (blue).

But I would like to increase the resolution within the existing data without shifting results in the $x$-direction?

MWE:

close all; clear all; clc

fs=512; % sampling frequency

load data.mat

t=X{1}; % time
y=X{4}; % accelelation

y=y./1024; % calibration of data

%% resampling

%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%
ups=10; % upsampling rate
%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%

dns=1; % downsampling rate
fu = fs*ups; % upsampling frequency
% tu = t(1):1/fu:t(end)-1/fu;
ytu=resample(y,ups,dns);

%% FFT and plot

figure;

x=y;
[N,m] = size(x);

freq = 0:fs/length(x):fs/2; %frequency array for FFT
xdft = fft(x); %Compute FFT
xdft = 1/length(x).*xdft; %Normalize
xdft(2:end-1) = 2*xdft(2:end-1);
semilogy(freq,abs(xdft(1:floor(N/2)+1)),'k-')
hold on

x=ytu;
[N,m] = size(x);

freq = 0:fu/length(x):fu/2; %frequency array for FFT
xdft = fft(x); %Compute FFT
xdft = 1/length(x).*xdft; %Normalize
xdft(2:end-1) = 2*xdft(2:end-1);
semilogy(freq,abs(xdft(1:floor(N/2)+1)),'b.')

xlabel('Frequency (Hz)');
ylabel('Accel (g)');
grid on;

Here is the link to data.mat.

Your whole problem is due to the fact that you have made a wrong assumption.

Increasing the sample rate does not increase the frequency resolution. It will increase the Nyquist frequency (the highest frequency you can represents) while keeping the same resolution. And that is exactly the phenomenon you have encountered with Matlab.

In other word, after upsampling by a factor of 10, your spectrum goes up to $2560 \text{ Hz}$, not $256\text{ Hz}$. I think you figured that out as the script you provided takes that in consideration while the figure you have posted does not.

In order to improve the resolution, you have to take a longer acquisition in time.

Your frequency resolution is evaluated as :
$$R = \frac{F_s}{N}$$
Where :

$R$ is the resolution
$N$ your number of FFT samples
$F_s$ your sampling frequency

In the script you have provided, you did increase the sampling frequency by a factor of $10$, but you also increased the number of samples by the same factor. Thus, you kept the exact same resolution of $512/90000 = 0.568 \text{Hz}$.

In order to get a better (smaller) resolution you have to:

Increase $N$ without changing $F_s$. In other word longer acquisition
Decrease $F_s$ without changing $N$. In other word longer acquisition

Hope that helps!


Direct Form II , IIR filter I have read in my textbook and i don't understand this part, that the Direct Form II implementation of an IIR filter is a bit tricky in memory addressing.There is need of clever memory addressing.  Memory is the place that we store partial summations(states of the filter). Anyone that understands the problem and explain to me?

The tricky part tends to be time shifting the state variables. You try to avoid spending extra cycles on explicitly doing this but ideally this can be all kept in registers or just happens during the natural flow of the instructions. The best way to do this will highly depend on your instruction set and/or specific capabilities of the hardware.

In general Direct Form II is a poor choice. The transfer function between the input and the state variables is vine solely through the poles which leads frequently to stability and noise problems. Transposed Form II is almost always the better choice. Only exception would be a fixed point processor with a double precision accumulator in which Direct Form I would be the preferred choice.


Distances In a Single Image With Some Real References Consider following image:

We have some real references here: Stonehenge measurements (stone dimensions and distances), some strait lines at the back of carts, cart wheels that are almost of the same size for both carts. At the same time we don't know the specification of the camera lens.
Question: Is it possible to calculate distance of the cart from stones? Or can we even deduce that the cart belongs to the scene or not?

Well, you can have an approximation. Due to the forms in the image are non regular, completely plain, its hard if not impossible to know the camera distortion. But if you know the size of the wheels and also you have the measures of the stones you can approximate the distance of the cart to the stones by using this equation:

$$
D_w = \frac{S_w * f}{P_w},
$$

where: $D_w$ is the distance from the camera to the wheel, $S_w$ is the real height of the wheel, $f$ is the distance of the focal to the film when the film is in the camera (You need to search for this parameter to use this approach, also you need to know the proportion of the real impressed picture with the film) and $P_w$ is the height of the wheel in the film. The same thing for the stone:

$$
D_r = \frac{S_r * f}{P_r},
$$

where: $D_r$ is the distance from the camera to the stone, $S_r$ is the real height of the stone and $P_r$ is the height of the stone in the film. Finally the distance $D_t$ from the wheel to the stone is:

$$D_t= D_r -D_w$$

I recomend you to use the objects that are marked in red in the next image, because are the nearest to the epipolar line.

Hope this could help, Best.


How to achieve an arbitrary phase shift? To help demonstrate why a time delay appears as a linear phase shift and not a constant phase shift I wanted to plot what a constant phase shift looks like.

I have been looking at the Hilbert transform, which shifts the phase (of positive frequencies) by $-90^\circ$, but I was unable to find a transform or technique that allowed shifting the phase by an arbitrary amount.

How can I shift the phase of a signal arbitrarily?

You can't design a filter that creates a phase shift that's constant with frequency for  real valued input (if that's what you are trying to do).

A Hilbert transformer appears to be doing this. However, the problem is, you can't implement a perfect Hilbert transformer since it's non causal with an infinite length impulse response.

The tricky part is that the DC and Nyquist components of the spectrum have to be real valued and that the phase at these frequencies is constrained to $0$ or $\pi$. Any phase shifter would have to have a transition band from DC to the to target phase and from the target phase to Nyquist.

That's exactly what happens if you design a real-world Hilbert transformer: you trade off the size of the transition bands against the length & complexity of the filter.

The class of filters that only manipulate phase and not the amplitude are called allpass filters. One can easily show that

Allpass filters have a phase that's monotonously decreasing
For an allpass filter of oder $N$ the phase from DC to Nyquist decreases by $N \cdot \pi$
Allpass have zeros that are inverses of the poles,
The z-transform numerator polynomial is the reverse of denominator polynomial

Of course you can always design an "approximation" that's good enough for your requirement. You'd have to specify the target bandwidth, max phase error, max amplitude error and then deploy a suitable FIR or IIR design method

Even an approximate design of absolute phase shift is quite difficult, especially if you have latency or causality constraints. Constant relative phase shift is a lot easier. A popular technique is to run a signal  through two parallel allpass filters which are designed to have a constant phase difference over the frequency area of interest. The outputs have identical amplitude as the input but have approximately constant phase shift with respect to each other (although not with respect to the input)


systematic way to determine best bit mapping policy So in class, we've been asked to find the best bit mapping policy for QAM modulation schemes (4-QAM, 8-QAM, 128-QAM, etc). Best meaning the one that will provide the lowest bit error Probability ($P_{be}$). I think I understand the concept, which is for each point closest to a given point, try to change the fewest number of bits. The absolute best case is gray coding, but that would only be possible for 4-QAM.

For example, I can come up with the best bit mapping policy for the star 8-QAM constellation below. ( the points are numbered in red)

for outside points (points 1,3,5, and 7) I'm making sure that the two closest points to each only 1 bit changes. So point 1 = 000, and point 2 = 100 (only the Most significant bit changed) and points 8 = 001 (only the least significant bit changed). For the 8-QAM, it's small enough that Just pick a random bit combination for the outside points, and then try to only change 1 bit to get to the center, but for a 16-QAM, for even a 256-QAM this starts getting difficult to just go about picking at random.

So my question is this; Is there a systematic way to go about selecting the bit mapping for these larger constellations?

Also, once we find the bit mapping policy, how can we calculate the bit error probability $P_{be}$? I know how to find the symbol error probability $P_{es}$, and i know for gray coding $P_{be} = \frac{P_{es}}{log_2 M}$ where $M$ is the number of symbols (for 8-QAM $M = 8$). but since everything higher than 4-QAM can only be partial gray coding, how do we determine $P_{be}$ in terms of $P_{es}$?

For $2^{2n}$-QAM, Gray code labeling often assigns the odd-numbered bits to one axis (I or Q) and the even-numbered bits to the other axis (Q or I) and uses $n$-bit Gray code labeling along each axis. Alternatively, the first $n$ bits could be assigned to I (or Q) and the last $n$ bits to Q (or I). Depending on the structure of the modulator (e.g. 8 bits entering the modulator on a serial line versus 8 bits entering the modulator simultaneously on a 8-bit-wide bus), one method might be preferable over the other. It is important to note that it is not necessary to convert the data into Gray code; just to use the Gray code labeling.  An extensive discussion of this last point can be found here on Google Groups.


Can we start a wave $\sin(2x-4)$ from position 4 in the horizontal axis? Even though the wave of $\sin(2x-4)$ looks like this and it starts from position 2 out of curiosity I would like to know whether we could start it from position 4 like the wave of $\sin(2x)$ shown here.

You can shift the origin to $x=a$ by replacing the argument $x$ by $x-a$:

$$\sin(2(x-a))\tag{1}$$

Since for $\sin(2x-4)=\sin(2(x-2))$ you have $a=2$, the origin is shifted to $x=2$. If you want to shift it to $x=4$, you get from Eq. $(1)$

$$\sin(2(x-4))=\sin(2x-8)\tag{2}$$


Can you explain how to imagine frequency in 2D? I need to understand how the concept of frequency exists in images The concept of frequency spectrum exists because it is possible to sum sines and cosines of various frequencies together to make any periodic signal. The frequencies of these individual sinusoids make up the spectrum of the final signal. Now when it comes to image processing we treat image as a signal similar to sound for example in that it has a bandwidth and it is possible to filter certain frequencies it to achieve various results.

What I am not clear about is, how did people get concept of frequency in 2D? Sinusiods are 1D signals, how can we extend this concept of frequency into 2D? Is it possible to do the same in 3D and do filtering in 3D? Does anyone do that in any application?

I fully understand low pass, high pass, band pass and band stop filters. Their frequency respons curves are quite straight forward. What I do not understand is, how is this extended into 2D?

Frequency, by its definition, means how frequently something occurs. You can treat images as planes where you have data points, ie. the pixels of the data. In images a frequency can be thought of as how frequently the sinusoidal components of the structure repeat per unit distance (where your unit would be the distance of a pixel).


Why oscillations in PI control? When integral control is added to proportional controller (i.e. PI control), why steady state error becomes 0 and  why oscillations appear , how actually the integrator works ?

For symplicity, consider the SISO linear system
\begin{align*}
\dot x &= ax +bu\\
y &= cx
\end{align*}
with $x$,$u$ and $y$ taking values in $\mathbb R$. Assume that you want to stabilize the system by making $y$ to converge to some fixed set point $\bar y$. Then necessarely there must be a point $(\bar x,\bar u)$ such that
\begin{align*}
0 &= a\bar x+b\bar u,& \bar y=c\bar x
\end{align*}
since $\bar x$ must be an equilibrium and such that $y=\bar y$. Trivially one has $$\bar u = -\dfrac{a}{bc}\bar x$$ One can define the error variables as $\tilde x=x-\bar x$, $\tilde y=y-\bar y$ and $\tilde u=u-\bar u$ that yield the error system
\begin{align*} \dot{\tilde x} &= a\tilde x + b\tilde u,& \tilde y&=c \tilde x
\end{align*}
One can stabilize it by taking $\tilde u = k\tilde x$ such that $a+bk<0$. That in fact ensures that asymptotically $\tilde y=0$ which implies $y=\bar y$. The overall control law is
$$
u = \tilde u + \bar u = k(x-\bar x) -\dfrac{a}{bc} \bar x = kx + \dfrac{1}{c}\left(-k-\dfrac{a}{b}\right)\bar y
$$
Thus as you can se, simple proportional is not enough, you also need a feed-forward control action given by
$$
u_{ss} = \dfrac{1}{c}\left(-k-\dfrac{a}{b}\right)\bar y
$$
If you define an integrator system like $\dot\sigma = y-\bar y = C(x-\bar x)$ then also $\sigma$ will have a set point $\bar \sigma$, and with $\tilde \sigma=\sigma-\bar\sigma$, the extended error system becomes
\begin{align*}
\begin{bmatrix}\dot {\tilde x}\\\dot{\tilde\sigma}\end{bmatrix} &= \begin{bmatrix}a & 0\\c & 0\end{bmatrix}\begin{bmatrix}\tilde x\\\tilde{\sigma}\end{bmatrix} + \begin{bmatrix}b\\0\end{bmatrix}\tilde u, & \tilde y&=C\tilde x\end{align*}
Take $\tilde u=k_p \tilde x + k_i \tilde\sigma$ (proportional+integral), then the closed loop system reads
\begin{align*}
\begin{bmatrix}\dot {\tilde x}\\\dot{\tilde\sigma}\end{bmatrix} &= \begin{bmatrix}a+bk_p & bk_i\\c & 0\end{bmatrix}\begin{bmatrix}\tilde x\\\tilde \sigma\end{bmatrix}, & \tilde y&=C\tilde x\end{align*}
You can chose $k_p$ and $k_i$ to make the closed loop matrix to be Hurwitz, and thus to stabilize the error system and you don't need any feedforward action, in fact the control law is
$$
u = \tilde u + \bar u = k_p(x - \bar x) + k_i (\sigma -\bar \sigma)
$$
Then with $$ \bar \sigma = \dfrac{k_p}{k_i}\bar x$$ the control law to becomes
$$
u = k_p x + k_i \sigma \qquad\qquad (1)
$$
and that's all. Moreover without integral control you need a feed-forward action which depends on the system parameters. If you have uncertainties on the parameters for sure you wil have an error on your set-point. with the integral action your final control law (1) does not depend on the parameters! and as long as the closed-loop system is stable you have null error at the steady state. Oscillations depend on how you chose $k_p$ and $k_i$. In fact if you chose them such that the closed loop matrix has complex eigenvalues then you can have socillations. For MIMO system the reasoning extends trivially.


Auto Correlation vs Cross Correlation vs Convolution and their applications I know from wikipedia that auto correlation in done on the same signal while cross correlation is done on different signals.But what does this actually imply in terms of application.I can always apply cross correlation on same signals and get same output.
And in convolution one signal is reversed.Mathematically i understaand the formulas.

But what does these three mean in terms of applications?

I can tell you of at least three applications related to audio.

Auto-correlation can be used over a changing block (a collection of) many audio samples to find the pitch. Very useful for musical and speech related applications.

Cross-correlation is used all the time in hearing research as a model for what the left and ear and the right ear use to figure out a sound's location in space (this is called sound source localization). In the case of two microphones you would cross-correlate the left channel with the right channel.

Convolution is used in simulating reverberation. A room's impulse response can be determined from measurements and that impulse response can be convolved with any sound source to simulate the reverberant response (at the impulse response recording's exact location).

I know this answer isn't complete but maybe it can give you some idea that there is in fact a practical use for auto- and cross- correlation!

So in general, auto-correlation can be used to extract properties of a signal, cross-correlation can exploit the information between two related signals, and convolution can be used to modify the properties of an incoming signal based on some time, frequency, and phase response specified by the impulse response you're convolving the source with.


Highpass filter applied to noise Suppose the input signal is represented as $$y=x+w$$ where $x$ is the signal to be estimated and $w$ is AWGN. I design a high pass filter $H$. What does $$H^{T}Hw$$ represent? Will it be just AWGN but attenuated or will it AWGN completely distorted. We can also assume that $\sigma_w = 0.2$

The formula represents noise component of the output signal. AWGN stands for additive white gaussian noise, than after passing through the filter the noise will still be additive and gaussian but NOT white.


Image Processing Prior to Feature Detection I have implemented a feature detector based on Harris corners. It works fine most times, but there are cases where it performs poorly. I need to make it work on many different images without configuring it individually.

The problem is with the detector threshold value. If set too low, the detector fires too many times resulting in huge number of features. If set too high, there are too less features.

I have partially resolved this by ANMS (Adaptive Non-Maximal Suppression) to reduce number of features before assigning a descriptor vectors to them.

However, images like this one are the problem:

They have low contrast and I cannot "afford" setting threshold too low for all images. It would cause detector to work on these images, but other images would contain hundreths of thousands of features, which will be slow to filter with ANMS and this will harm overall performance.

I was thinking of adjusting the image before feature detection. Maybe histogram equalization would do the job. This may be a valid operation since global change of contrast have no effect on feature descriptors (they are invariant to change in brightness and contrast).

Maybe working with adaptive threshold or some heuristic would work better.

Any other suggestions?

Do you really have to use Harris corners? There is many features developed after Harris corners, with better properties. A good overview can be found in this article:

A.L. Dahl, H. Aanæs, and K.S. Pedersen. Finding the Best Feature Detector-
Descriptor combination

Based on that article as well as my personal experience, I would suggest switching either to MSER (Maximally Stable Extermal Regions), or even combine them with DoG (Difference of Gaussians) -- the features first presented as a part of the SIFT pipeline.

If the problem really is in low contrast, then MSER features should really make you happy: they are (fairly) invariant to changes in lighting. In short, they are connected regions of the image stable through a series of different threshold binarizations.

The feature extraction process is independent from calculating the descriptors, so it shouldn't be too hard to integrate new ways of feature extraction in to your process.

Also, I've heard of (but never actually worked with) Multiscale Harris corners as an extension to Harris corners. I don't know much about them and personally can not recommend any reading materials on this topic, so I leave article search and picking the most interesting materials to you.

Furthermore, might I suggest that the image you posted might have other problems than low contrast. In my personal experience, vegetation like bushes or possibly the field you have, as well as the lovely bubbly clouds tend to produce "generic features" -- features which tend to have equally similar (or dissimilar) descriptors as a lot of other features.

Practically, this means that when doing feature matching on two images from a different perspective, features extracted from these kinds of surfaces tended to be falsely matched.  I have done a Master thesis that in a large part deals with feature extraction to be used in feature matching further used to calculate a homography transformation between two images when I came across this problem. I didn't find any other articles describing this problem at the time, but my thesis might be helpful for your overall approach.

Lastly, as you have set, thresholds and techniques that work just fine on most images extract to little features in this kind of images, because of its mostly homogeneous areas. This kind of images present problems in feature matching (which can be extended to image stitching), content based image retrieval, and I would presume tracking as well as similar applications. No method currently works quite well on them.

Methods that work good on this kind of images as well as the typical cases are being explored and researched currently, such as an approach I started working on briefly described in this answer.


Discrete Signal Windowing How do we apply windowing of 50% (4 samples) to a discrete signal with the odd number of values as 1234567?
Is it done as: 1234 3456 567(0)
or 1234 3456 67(0)(0)?

Depending on the window function, you may be able to use a DFT-even version of the window function. "DFT-even" means that the periodic extension of the window function is symmetrical. In MATLAB and Octave you can get such a window like this (the first line in the source code):
a = hanning(10, "periodic");
b = fftshift(a);
c = a + b;
plot(a, "x", b, "o", c, "+");

Figure 1. DFT-even Hann window (blue x), shifted DFT-even Hann window (red o) and their sum (orange +).
The window function and the shifted window function have in some sense a 50 % overlap, and they sum to a constant, which I think was the desired quality that the 50 % overlap was though to facilitate.
Not a problem in the above, but if the window function is such that it does not end at a zero value, and the start and the end values are not equal, then one possibility would be to move one half of the value of the last sample to just before the first sample.


How to approach detection of anomalies in a video feed? As part of a project I have a video from a quadrotor that captures the road and the environment around it and I want to identify and detect potholes and speedbumps on the road using image processing. What is the best approach for this?

Here is what I've done until now.

The video feed is converted into grey scale and I use one particular portion of it(marked by the blue rectangle) for processing. In that window, I again get smaller windows of and I calculated their corresponding Grey Level co-occurrence Matrix (GLCM henceforth) of all those windows. I calculated the textural features from the GLCMs like ASM, erratic homogeneity etc and I plotted them as well.

I did this with the hope that there will be some kind of obvious change in the features when the blue rectangle reaches the speed bump(which is slightly ahead in the image). But the results erratic and I cannot make a definite decision.

If there is any other approach, or if there's any change I could implement to get better results in the same method do let me know.

Here is some brainstorming:

You have two interesting variables here, I would say:

spatial changes in y-direction (left to right)
temporal changes (t-direction).

I would not directly include the x-direction (up-down direction) in the first try, since this is correlated to the t-direction.

Try the following:

Calculate the y-gradient
Calculate the t-gradient
Create a binary mask from both using two threshold values (you probably have to figure them out in a test run), so that the whole image is zero where there is no change and one where there is a change.
To suppress spurious noise in the binary mask, you could median-filter it with a 5-by-5 or 7-by-7 median filter (try out the best size)
Use a morphologic dilate filter to close "holes" in the binary mask and also to dilate the masks so that they do have a broader overlap

This way you are directly using the visual features of the image - changes in the gray values. Now a few things can happen with those two masks

Areas of the t and y masks overlap: You have a change in the image in y direction and there is also a change between to consecutive images - you found a localized change that could be a hole or a speedbump
There is a change in y direction, but no or just a small area in t direction: That could be a road marking (the algorithm might fail if your drone tumbles quite strongly and your camera has a low number of frames per second - then the t-gradient might be quite high in this case. To suppress those signal changes from lane markings, you could suppress lines in the binary mask by a Hough-transformation of the mask and kicking out pairs of dots)
There is change in t but not in y - probably a lane mark that started/stopped

If this really works with your data, I do not know. However, starting with simple things can lead you into a quite robust and reasonably fast direction for your image processing task.


Modelling the amplitude and phase through an amplifier I am writing a software simulator to simulate a real physical acoustic interferometer. I now want to write the module to model the microphone amplifier(s). I'd like to be able to specify the frequency response of the amplifier both in terms of amplitude (gain) and phase. So I will feed in a sampled waveform and get a sampled waveform out that is modified according to the properties of the amplifier circuit.

I'm no expert in signal processing but happy to read the background maths. But to save me time going down dead end ideas could someone please tell me the standard method of doing this and what maths I have to read up on to be able to understand it?

I apologize, the following is a bit rough, but I've put together an example in MATLAB illustrating all the different elements I think you need. That is:

a fist-order model for an amplifier (a common approximation for most op. amps.)
how to discretize it to find an IIR filter
one way to implement a realization of the IIR filter (straight forward state-space implementation)
how to generate a transfer-function estimate from measured data
one way (not the most robust, but perhaps the most obvious) to fit the transfer-function coefficients to the transfer-function estimate

To read up on the theory behind all this I would recommend the following text-books:

Chi-Tsong Chen, Linear System Theory and Design (3rd edition). Oxford University Press, 1999.

Petre Stoica and Randolph Moses, Spectral Analysis of Signals. Prentice Hall, 2004.

Lennart Ljung, System Identification: Theory for the User (2nd edition). Prentice Hall, 1999.

The last book is somewhat heavy, but concentrate on "output error" methods, and you should be very well covered for systems where you know the order of the model, or "subspace" methods if you don't. These methods are available in the MATLAB system identification toolbox.

%% Amplifier model
wb = 1e4; % Amplifier bandwidth
ga = 1; % Amplifier gain

% Modeling amplifier as having fist-order dynamics
% If using Octave run these in the command window:
%pkg install -forge control
%pkg load control
H = tf(ga,[1/wb 1]); % Continuous-time transfer function

%% Finding a realization (convert to state-space formulation)
% V/U = ga/(s/wb + 1)
% V*(s/wb + 1) = ga*U
% s*V = -wb*V + ga*wb*U
% dv/dt = -wb*v + ga*wb*u
%
% Standard form:
% dx/dt = A*x + B*u
% y = C*x + D*u
%
% with x = v, y = x
A = -wb; B = ga*wb; C = 1; D = 0;

%% Finding a discretization:
Ts = 1/1e5; % sampling time
% Using Euler integration
Ad1 = 1 + Ts*A; Bd1 = Ts*B; Cd1 = C; Dd1 = D;
% Using exact discretization (use when system has analog-digital converters)
Ad2 = exp(A*Ts); Bd2 = (Ad2 - 1)/A*B; Cd2 = C; Dd2 = D;

%% Computing a step-response to demonstrate one way to implement an IIR filter
u = [zeros(1,100) ones(1,200)]; % input
x = zeros(2,numel(u)+1);
y = zeros(2,numel(u));

t = 0:Ts:(numel(u)-1)*Ts;

% Comupting filter outputs
for k = 1:numel(u)
    x(1,k+1) = Ad1*x(1,k) + Bd1*u(k); % response of Euler integration scheme
    y(1,k) = Cd1*x(1,k) + Dd1*u(k);

    x(2,k+1) = Ad1*x(2,k) + Bd1*u(k); % response of exact discretization
    y(2,k) = Cd1*x(2,k) + Dd1*u(k);
end

figure(1)
plot(t,y(1,:),t,y(2,:))

%% Computing the response to noise for identification of parameters
% The frequency response of a system is in general most reliably found when
% exciting the system with white noise. A pseudo-random binary sequence is
% theoretically optimal when you don't know anything about your system, but
% can in practice cause problems due to the switched nature. A Gaussian
% noise source is usually a good substitute.
un = randn(1,1e4); % input
xn = zeros(1,numel(un)+1);
yn = zeros(1,numel(un));

% Comupting filter outputs
for k = 1:numel(un)
    xn(k+1) = Ad2*xn(k) + Bd2*un(k);
    yn(k) = Cd2*xn(k) + Dd2*un(k);
end
% Remove transient as you want to look at stationary response
un_ = un(201:end);
yn_ = yn(201:end);

% Compute power spectral densitites and cross-power spectral densitites
% using a simplified version of the Welch method
N = length(yn_); % length of original sequence
L = 5; % number of PSD averages
M = floor(N/L); % length of sequence segments
W = 0:(2*pi)/M:2*pi-(2*pi)/M; % normalized PSD frequencies

WIN = hanning(M); WIN = WIN(:).'; % optional windowing (if the segments are short)
Pwin = 1/M*sum(abs(WIN).^2); % window "power" correction (for Welch method)

% PSD for input sequence
Puu = zeros(1,M);
for k = 1:L
    un_segment = un_((k-1)*M+1:k*M);
    un_segment = un_segment(:).';
    un_windowed = un_segment.*WIN;
    Uft = 1/sqrt(2*pi*M)*fft(un_windowed);

    Puu = Puu + (Uft.*conj(Uft));
end
Puu = 1/(L*Pwin)*Puu;

% Cross-PSD between input and output sequence
Puy = zeros(1,M);
for k = 1:L
    un_segment = un_((k-1)*M+1:k*M);
    un_segment = un_segment(:).';

    yn_segment = yn_((k-1)*M+1:k*M);
    yn_segment = yn_segment(:).';

    yn_windowed = yn_segment.*WIN;
    Yft = 1/sqrt(2*pi*M)*fft(yn_windowed);

    un_windowed = un_segment.*WIN;
    Uft = 1/sqrt(2*pi*M)*fft(un_windowed);

    Puy = Puy + (Yft.*conj(Uft));
end
Puy = 1/(L*Pwin)*Puy;

% Transfer-function estimate
Huy = Puy./Puu;

figure(2)
subplot(2,1,1)
semilogx(W,abs(Huy))
subplot(2,1,2)
semilogx(W,unwrap(atan2(imag(Huy),real(Huy))))

%% Fitting a model to the output data

% finding the z-transform of discrete-time model (Euler integration)
% using the standard expression
z = tf('z',Ts);
Hd1 = Cd1*(z - Ad1)\Bd1 + Dd1;

% finding the z-transform of discrete-time model (exact discretization)
% using an alternative method
b0 = Bd2*Cd2;
a0 = -Ad2;
a1 = 1;
Hd2 = tf(b0,[a1 a0],Ts);

% compare the continuous and discrete-time models
figure(3)
bode(H,Hd1,Hd2)

zv = exp(1i*W); % starred transform (between z- and s-domains)

% H = b0/(a1*z + a0)
% H(W) = b0/(a1*exp(i*W) + a0)
% a1 = 1; % in order to have unique coefficients for a0 and b0
% H(W) = b0/(exp(i*W) + a0)
%
% (exp(i*W) + a0)*H(W) = b0
% exp(i*W)*H(W) = -a0*H(W) + b0
% exp(i*Ts*w)*H(W) = [-H(W); 1]*[a0, b0]
% zeta = phi*theta;

zeta = zv(:).*Huy(:);
phi = [-Huy(:), ones(size(Huy(:)))];

theta_true = [a0; b0]
theta_est = phi\zeta


Spectrum Analysis using Windowed FFTs I have a couple of questions regarding windowed FFTs:

Why is the noise floor higher with windowed FFTs (according to Wikipedia's spectral leakage page, anyway), when the whole point of windowing is to reduce side lobes?
I realize that different windows are better for different things, but is there a window that is considered to be the best all-around window for spectrum analysis?  Alternatively, is there a better way to do spectrum analysis than windowed FFTs?  It would have to be a DSP approach (i.e. I can't do an array of analog filters), but within that constraint I am game for different solutions.

Thanks for your time in reading this.

A non-rectangular window will remove "noise" from distant bins at the cost of adding more "noise" to the immediately adjacent bins to a narrow-band spectrum peak.  The sum of both these spectral leakage effects is greater than zero for a non-rectangular window. So if you count the raising of the level total of all adjacent bins as noise, then the S/N ratio is lowered.

Some people don't care about the bins immediately adjacent to a spectrum peak (their spectral peaks are a priori assumed to be widely spaced; and/or they interpret, interpolate, or phase-vocoder adjust the energy out of those adjacent bins back into the central peak bin), so for those purposes, the reduced far-side-lobe energy means less noise.

Another reason for a lower S/N ratio is that windowing of quantized data is an informationally lossy process, and these (re)quantization losses can also be considered a form of noise.

"Best" is relative to some weighting of quality metrics, and different users may have very different weightings.  If you don't have a set of prioritized design goals for which to optimize a window, then you may not have a strong reason to not just use a Von Hann window.

Depending on your data source and your needs, using just some windowed FFTs may not even be a good form of spectrum analysis, much less the best possible.  Or the opposite.


OFDM use a pulse shaping filter or not? I've heard in various literature that band limited channels use pulse shaping filters. By the way, In OFDM it use a cyclic prefix to overcome ISI, so I wonder if OFDM does not use pulse shaping filter.

The purpose of pulse shaping filters is not to overcome ISI as is implied in the OP's question. The only reason for using a pulse shaping filter is spectral efficiency, and in the process ISI can be introduced if not done properly. In order to limit bandwidth, pulse shaping must extend the time domain response for each pulse beyond a symbol boundary, but can be done such that the response passes through zero at the symbol decision location of ever other pulse (Nyquist shaping filter, raised-cosine filter etc).
A rectangular pulse in time is a Sinc function in frequency, with the main lobes of the Sinc that roll-off at only 1/f (20 dB/decade) as shown by the blue spectrum for a single-carrier rectangular pulse modulation in the plot below. The red spectrum shows the significant advantage is spectral containment by pulse shaping:

Notice this about the plot above to see intuitively why we are less concerned with using pulse shaping for OFDM (or when we may want to consider its use). The blue spectrum above is rolling off at 1/f where f is proportional to the symbol rate. If we reduce the symbol rate for a single carrier by a factor $N$, we equivalently increase the rate of frequency rolloff by that same factor. OFDM converts a single carrier modulation into $N$ sub-carriers, each representative of the spectrum above in blue, each with much reduced bandwidth individually and therefore overall the resulting waveform extends much less in frequency beyond its primary occupied bandwidth as given by $NR$ where $R$ is the symbol rate of each of the sub-carriers.
So to recap that with the plot of the blue spectrum above as an example: Here we see the result of a 1 MHz symbol rate (1 MSymbol/sec) single-carrier modulation, and notice within the range of the plot the spectrum rolls off to about -35 dB at 5 MHz offset. Consider if we instead sent the same amount of data by using 64 carriers, all reduced in rate by 1 MHz/64, and spaced in frequency by the same (with is orthogonal carrier spacing). The resulting spectrum would dominantly be in a 1 MHz span of frequency, with the carriers on the ends contributing the most to the spectrum going beyond that range, and specifically we would expect about -35 dBc at a 5 MHz/64 =  78.125 KHz offset. That, without any additional pulse shaping, wouldn't appear very different from the red plot above (which is a single carrier pulse-shaped waveform)! If we wanted to make the spectrum even tighter (with the diminishing return of benefit), we could certainly apply pulse shaping to reduce the adjacent leakage terms accordingly.


What are the options for turning a bit stream into a byte stream? I am creating a system (modulator and demodulator) for broadcasting binary data through an analog signal. I'm using BPSK with great success – at least, I'm generating great graphs.

Green is input signal, red is carrier phase, pink is digital zeroes and ones.

My question is – now that I have that pink stream of binary data, how to delineate bytes? I have been reading about clock recovery, 8b/10b, Manchester, NRZI... And I'm a little confused about what those offer me.

I imagine I could just come up with a system like Manchester, where each 0 is just 01 and each 1 is 10 – and then use a pattern like 0000 to mean "start of next byte" – but this doesn't seem right. Under that system, how would you interpret 10000010? It could mean the last byte ended with a 1 and the next byte starts with a 1 – or it could mean the exact opposite – the last byte ends with a 0 and the next byte starts with a 0.

I know there must be some basic DSP concepts that could help me here, please point me to some resources!

I would use a frame alignment word to synchronize your data.  The basic idea is to insert your FAW at the beginning of every $N$ bytes.  After you demodulate the data you look for the FAW to synchronize the data.  To avoid false positives you should try to have a longish FAW (16 bits should be plenty) and look for repeated FAWs.  To avoid inefficient overhead you should try to make $N$ large compared to the size of the FAW.


Problem in understanding DDFS (Direct Digital Frequency Synthesizer.) I need to generate sine and cosine signals using DDFS (Direct Digital Frequency Synthesizer.) I'm unable to get the idea of the phase increment applied as input which changes the frequency of the quadrature signals.

A Direct Digital Frequency Synthesizer (DDFS) or simply Direct Digital Synthesizer (DDS) usually refers to the combination of a Numerically Controlled Oscillator (NCO) combined with Digital to Analog (D/A) converters for complex quadrature output or a single D/A converter for a real output.

An NCO is actually a fairly simple structure, consisting of a Frequency Control Word (FCW) input to a "phase accumulator" which is nothing more than a simple counter, usually of high precision, and counts in increments of the FCW. Note that an accumulator is a digital integrator (for example in continuous time, the integration of a constant is a ramp; in discrete time, an input of all ones into an accumulator [1 1 1 1 1 ...] would also be a ramp [1 2 3 4 5 ...].  For this reason, if we say the input that is just a digital value represents frequency, and knowing that the integration of frequency is phase ($freq = \frac{d\phi}{dt}$), then we can therefore see that the accumulator output represents phase (so we call it the phase accumulator).  Since the output of the counter is phase versus time, we then need a function to convert that to a sine wave (or sine and cosine wave for quadrature output). This is often done with a lookup table (LUT), that will function as a calculation of the sine and cosine values for the specific phase values at a given time in the phase accumulator. The output of the LUT will therefore be the (digital) sine and cosine functions versus time, with a rate that is set by the input FCW (A higher FCW will result in a higher count rate at the phase accumulator output, meaning more cycles in time of our sine and cosine waveforms at the LUT output.  The LUT will contain exactly one cycle of the sine and cosine waves (effectively, memory compression algorithms are actually used so that less data needs to be stored), so as the phase accumulator rolls over properly ($Cnt_{max}+n = n$), the output sine and cosine waves continue where expected without interruption.

Here are some slides I have that describe the NCO. Also included in the figure is the Phase Control Word (PCW), which is an offset control for the Phase Accumulator Output, providing direct phase modulation control if desired:

NCO Implementation view:

NCO Mathematical view (continuous time equivalent):

In implementation, the phase accumulator is typically quite large in precision (24 to 48 bits), driven by the frequency step size desired. The frequency control word is of a similar size (usually 1 bit less to synthesize all frequencies from DC to $F_{clk}/2$). Since the LUT contains one cycle of the output frequency, counting through every value in the accumulator means counting through every address in the LUT, and this would result in the lowest possible output frequency above DC. Counting by 2 would result in a frequency twice as fast, and so on. From this we easily see the relationship:

$F_{step} = \frac{F_{clk}}{2^{acc}}$

and

$F_{out} = FCW * F_{step}$

Where

$F_{step}$ is the step size in Hz

$F_{clk}$ is the clock rate that the accumulator updates in Hz

$acc$ is the accumulator size in bits

$F_{out}$ is the output frequency in Hz

$FCW$ is the Frequency Control Word in digital counts. [0 to $(2^{acc-1})-1$]

Phase Truncation, SFDR and SNR

If we passed all the phase accumulator bits to a LUT, the digital output for each frequency control word would be perfect, with the only error source being the precision of the sine wave stored in memory (defined by the LUT output bit width). However the memory requirements for such an implementation would be excessive, and given the resulting noise, unnecessary.  Therefore we truncate the LSBs of the phase, leading to truncation error in the frequency output. This results in spurs throughout the frequency spectrum, based on the repetition rates of the truncation patterns that are created. The relationship between the highest spur and the phase truncation comes out very nicely as 6.02 dB/bit where bit is the number of MSB bits that are passed on after truncation.   For example, in the figure shown, 14 bits are passed on after truncation, so the Spurious Free Dynamic Range (due to phase truncation errors) is 6.02*14 = 84.28 dB.  This means, that although there are lots of spurs all across the digital frequency spectrum (from DC to Fs) but the strongest of all the spurs, in this case, is -84 dB below the output signal level. (At least the strongest of the spurs due to the phase truncation that was added).

I had also evaluated the SNR (signal to noise ratio, where in this case the noise is due to the combined power of all spurs from phase truncation) and this came out to be

SNR due to phase truncation:

$SNR = 6.02dB/bit - 5.172 dB$

here, bit is the number of bits into the LUT

So this means for our example with 14 bits after truncation, the combined power of ALL spurs created from phase truncation would be (6.02*14 - 5.17) = 79.11 dB below the output signal.  The strongest of the spurs was 84.28 dB down as calculated previously, so this means all the other truncation spurs increase the total spurious noise power by 5.17 dB. Note that are values for FCW that will result in a spurious free output (spuriuos free of spurs due to phase truncation), all FCW values that keep the truncated bits always at 0 - an example of this specific to the figure below where there are 18 truncated bits, is an FCW of 1+[eighteen zeros] = 262144 decimal) . However for all other values, and sufficiently long data runs, these formulas will apply.

This relationship was derived by making use of the fact that the phase truncation error is a ramp or cyclically sampled ramp, which has a uniform distribution.  This is a useful relationship that can be combined with the well documented quantization noise of a full scale sine wave (referring to the digital output after the look-up table), as:

SNR due to output quantization:

$SNR = 6.02 dB/bit + 1.76 dB$

here, bit is the number of bits out of the LUT (or effective number of bit, ENOB of the ADC if considering the ADC output).

Quantization noise and phase truncation noise can reasonably be considered independent and uncorrelated, meaning the total noise for a composite SNR would sum in power.  Therefore the above two SNR relationships can be very useful in establishing the overall noise performance, leading to the precision requirements on the NCO both at the input and output side of the LUT, while the accumulator size itself is driven from the frequency resolution desired.


Designing a 2x interpolation filter for audio I want to perform 2x interpolation on an audio signal sampled at 44.1KHz by upsampling the signal by adding a zero after each original sample and then using a lowpass filter to interpolate the results.

Assuming an analogue elliptical filter (of order 10 or lower) as the basis for the design. What design specification would constitute a high quality audio interpolator that minimises aliasing?

For instance a passband and stopband at 0.45 * Nyquist and 0.55 * Nyquist respectively (Nyquist = 22.05KHz) with a stopband attenuation of -50dB and passband ripple of -0.1dB would be common for a low cost commercial DAC (I appreciate that digital to analogue conversion is a slightly different scenario). Would there be better values for doing high quality audio interpolation that minimises aliasing and doesn't truncate too many of the high frequencies before 22.05KHz.

For instance, is it a good idea to have the stop band above 22.05KHz for interpolation or does it not matter if the lowpass rolls off a little after 22.05KHz and what is acceptable?

@Keith: Here's the situation as I see it. Your original signal's spectrum is shown in panel (a) below. The solid lines are spectral energy and the dotted lines are spectral replications. After the zero-stuffing (upsampling by two), the stuffed sequence's spectrum is shown in panel (b). Notice the spectral image centered at 44.1 kHz that needs to be removed (attenuated). So after zero-stuffing you need to implement a digital filter whose passband cutoff freq is 0.45*22.05 kHz and whose stopband begins at 0.55*22.05 kHz. (I show that filter as the dashed line in panel (b).) But keep in mind, the filter’s Nyquist rate is 44.1 kHz (88.2/2 kHz) and not 22.05 kHz. The final spectrum you want is shown in panel (c).


What is the phase coded modulation? I'm trying to understand what is the "phase coded modulation" and how does it work.

As far as I know, Although the linear frequency modulated (LFM)
waveform continues to be the workhorse of modern radars, there is growing interest in the use of
phase-coded (PC) waveforms.  Phase coding can be used to reduce radio frequency interference (RFI)
between adjacent radars and make the waveform less vulnerable to exploitation, and more robust to
interference.

and I'm not a familiar with like these radio frequency engineering. So.

Is the Phase modulation same "phase coded modulation"?
I'm confused that because there are several types modulation way in
https://kr.mathworks.com/help/comm/ug/digital-modulation.html#fp45657

  Kind of Modulation
  Amplitude modulation
  Phase modulation
  Frequency modulation
  Continuous phase modulation
  Trellis-coded modulation

But I can't find anywhere about Phase coded modulation.

Would you help me to understand what is the "phase coded modulation" and how does it work with very very simple example for referencing?

Phase-coded modulation is another name for direct-sequence spread-spectrum (DS/SS) modulation. In a typical DS/SS communication scheme using binary phase modulation, a single bit is transmitted not as a sinusoidal pulse of duration $T$ and phase either $0$ or $\pi$ depending on whether the bit being transmitted is $0$ or $1$, but rather as a sequence of successive $N$ shorter sinusoidal pulses of duration $T/N$ with phases $0$ or $\pi$ as defined by a specific binary sequence (called a spreading sequence) $\mathbf b$ of length $N$ or its complement $\mathbf{\bar{b}}$ according as the data bit $d$ is $0$ or $1$. For example, with $N=7$, we would have a phase-coded transmission describable as
$$\begin{array}{ccccc}
d=0 &\leftrightarrow &\mathbf b = (1,1,1,0,1,0,0) &\leftrightarrow &(\pi,\pi,\pi,0,\pi,0,0)\\
d=1 &\leftrightarrow &\mathbf{\bar{b}} = (0,0,0,1,0,1,1) &\leftrightarrow &(0,0,0,\pi,0,\pi,\pi)
\end{array}$$
In the communication context, the interference between two different phase-coded transmissions using two different spreading sequences is determined by the cross-correlation between the spreading sequences. In the radar context, $\mathbf{\bar{b}}$ is not needed, but interference between two radar signals with different spreading sequences is still determined by the cross-correlation between the spreading sequences. The term spreading refers to the fact that the phase-coded signal has a bandwidth $N$ times larger than the non-phase-coded signal. Because the energy is spread over a larger bandwidth, the signal doesn't interfere as much with narrowband signals, and is not easy to detect even by broadband scanners. Thus, covert communication or covert scanning is easier, and signals are not that easy to jam, either by narrowband jamming or by broadband jamming, because the correlation receivers reject both kinds of jamming more easily.


I have a pressure signal and want to do SPL analysis on it Signal

I have an acoustic signal from a Ffowcs Williams Hawkings CFD analysis and would like to convert it to the frequency domain and see the SPL and OASPL. I know I need to use fft() but I am unsure why my plot differs so much from the one from FLUENT.

time = data{1};
p = data{2};

Fs = 1/(time(end)-time(end-1));
L = 0.887*Fs+1;
N = 2^9;
Y = fft(p,N);

for k = 1:N/2
   f(k) = (k-1)*Fs/N;
end
spl=20*log10(abs(Y(1:length(Y)/2)));
semilogx(f,spl,'k','Linewidth',2)

Steps are quite simple:

Apply the window to your signal
Calculate the DFT of your signal
Scale the DFT properly (half of spectra and coherent gain)
Convert the result to decibel scale with a reference acoustic pressure of $2\cdot 10^{-5} \mathrm{Pa}$

I don't have MATLAB with me right now so here is some Python code.
Considering 1 kHz sinusoid with amplitude of 1 Pascal (94 dB SPL). Resulting plot:

Source code:

#!/usr/bin/env python

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt

if __name__ == '__main__':

    plt.close('all')

    # Generate 94 dB (Peak) SPL sinusoid
    fs = 48000
    duration = 1
    npts = int(fs*duration)
    t = np.linspace(0, 2*np.pi, npts)
    f = 500
    ref = 2.0e-5 # 20 uPa
    amp = 1      # 1 Pa -> 94 dB SPL
    s = amp * np.sin((f*duration) * t)

    # Window signal
    win = np.hamming(npts)
    signal = s * win

    sp = np.fft.fft(signal)
    freq = np.fft.fftfreq(npts, 1.0/fs)

    # Scale the magnitude of FFT by window energy and factor of 2,
    # because we are using half of FFT.
    # To obtain RMS values, divide by sqrt(2)
    sp_mag = np.abs(sp) * 2 / np.sum(win)

    # Shift both vectors to have DC at center
    freq   = np.fft.fftshift(freq)
    sp_mag = np.fft.fftshift(sp_mag)

    # Convert to decibel scale
    sp_db = 20 * np.log10( sp_mag/ref )

    plt.subplot(2,1,1)
    plt.plot(t, s)
    plt.xlabel('Time [s]')
    plt.ylabel('Acoustic pressure [Pa]')
    plt.grid(True)

    plt.subplot(2,1,2)
    plt.semilogx(freq, sp_db)
    plt.xlim( (0, fs/2)  )
    plt.ylim( (0, 100))
    plt.grid(True)
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('SPL [dB]')
    plt.show()


ADC sampling rate implications Matlab/Simulink M-PSK modulation I am trying to find a way to evaluate what would be the minimum sampling rate used in a modulated waveform using matlab / Simulink.

The real scenario would be let say a M-PSK waveform of Y bandwidth in baseband reaches the point of the ADC and has to be sampled. Theoretically, the sampling rate should be 2xY and that (in ideal conditions) would allow you to restore what has been sampled.

In matlab/Simulink everything is “digital” meaning that we deal with waveforms produced by using random numbers generators. In such case by trying to introduce the process of sampling can only be realised as effectively down - sampling an existing waveform.

Shall I try to reproduce a more “continuous time” approach to get better estimates?

What is the usual coarse of action for such a problem?

Downsampling (taking every Dth sample and discarding the rest) is an identical process to sampling a continuous time signal; so you approach the continuous time signal as your sampling rate is increased. So you can access the impact of a lower sampling rate by simply downsampling a signal that has a higher number of samples.

Note that down-sampling and sampling does not distort a signal that only occupies a spectrum that is contained within 0 to $F_s/2$ for a real signal and 0 to $F_s$ for a complex signal, as well as any integer multiple of the ranges (where $F_s$ is the sampling rate. Distortion does come from aliasing and quantization noise. Aliasing occurs when there are additional signal or noise components in other frequencies prior to sampling (or downsampling) other than the range of interest, and quantization noise is the error from quantizing each sample directly (or in the case of a sampled/quantized signal that you are downsampling if you further round or truncate this would have the same quantization noise effect). This quantization noise is often well approximated as a uniform white noise process, meaning it's energy is equally distributed across all frequencies: If you increase the sampling rate, the same energy is spread across more frequencies, so the energy per unit bandwidth goes down. If your signal occupies a fixed bandwidth, a higher sampling rate will reduce the amount quantization noise in-band and thus increase SNR if the noise is limited by the quantization noise component (and assuming the out-of-band noise is further filtered digitally after sampling).

To combat aliasing we require anti-alias filtering prior to the sampling process. This would be a low pass filter when the signal of interest occupies the range of $0$ to $F_s/2$ (for the real case) or could be a bandpass filter when the signal of interest is at a higher frequency in which case the sampler also serves as the down-conversion.  As we make the sampling rate lower, we place a stricter requirement on filter design - but given a realizable filter with sufficient rejection of the signals that can alias the quantization noise effect detailed earlier would be the only performance impact.

There are also reasons of convenience to choose a particular sampling rate- for example it is common in many implementations to use a sampling rate of 2 samples per symbol- this works well with certain timing recovery implementations and allows for simple derivation of the symbol clock was the delay offset is corrected.

I further detail performance versus sampling rate in this post: What are advantages of having higher sampling rate of a signal?


Face Recognition / Image Processing For a Face Recognition system which one is better using Neural Network or the Histogram with pattern matching concept ?? which one is more accurate??

I need to integrate the face recognition engine to a Java project so I think it is hard to integrate a math lab project to a java project. So I need to know about Histograms and pattern matching in OpenCv (C++) and Neural Networks. Which one is better??

Well, it's important to understand that face recognition (just like object recognition process) is a two-stage process, the first one is the Computer Vision phase, which is to represent the image (face) in feature space that works well (depending on the task you want to do) so for example, the image itself is a set of ordered pixels (let's say 300x300), which makes each image representative by 90000 numbers, this is the dimensionality of your data points (images) and just FYI, 90000 is insanely huge, if you try to use Neural Networks, Support Vector Machines, or any other classification tool using this high dimensions, you will absolutely learn nothing - this is the curse of dimensionality (you will need to provide tremendously larger set for training, in other words, to learn any classifier you will need many more images than 90000, which is hard and unreasonable thing to do).

So you will need to represent the images in a different space (called feature space), for face recognition, there is something called eigenfaces (look it up), but the idea is that you find the most discriminate features of the faces then for each image, you convert that image to a set of numbers that will use the discriminant features to describe the image, so each image will be represented as a set of numbers (.3, .1, .1, ...) which will be much smaller than 90000. By the way (.3, .1, .1, ...) will mean something like for example .3 there's a difference between hair and face, .1 there's difference between face and background and so... However, this is a really hard step (feature selection) and it may vary based on the kind of images you have (like the basic eigenfaces only work when your images are aligned, are your images aligned?). Of course using Color Histograms is another approach to represent the image/face (personally (this is only my opinion), I think that color histograms won't be so efficient for face recognition, I might be wrong and would love to know if anyone knows)

The second step is the classification step, here you can use Neural Networks (which also needs a lot of parameter tuning), Support Vector Machines (popular), or decision trees (simple and interpretable), adaboost (popular), and some other approaches. - so you will provide your images in the feature space instead of the pixel-representative and hopefully you'll learn a good classifier that will classify with high accuracy.

Hope that helps


Splitting of signal into different components for calculation of PSD In BCI competition 3 dataset 5, It has been stated that we have to take 8-30 Hz frequency and break it with resolution of 2 Hz then it will result in 12 components but I am finding out only 11 components.Why??
The description can be found out here...

http://www.bbci.de/competition/iii/desc_V.html

Hint:

$$f_i=8+i\cdot 2\;\text{Hz},\quad i=0,1,\ldots,11$$

Note that $f_0=8\,\text{Hz}$ and $f_{11}=30\,\text{Hz}$. So you have 12 frequencies in the range $[8,30]\,\text{Hz}$, including both edge frequencies.


Sampling Theorem and Dirac Comb I am reading "The Scientist and Engineer's Guide to Digital Signal Processing" and trying to understand Figure 3.5 below which is about the sampling theorem and aliasing.

I do not understand the picture of the frequency domain. The signal can be represented as
$$f(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(x)e^{ixt}dx$$
where $F(x)$ is the Fourier transform of $f$. This represents $f$ with the "frequencies" $e^{ixt}$ and these can have $x$ negative. How then am I suppose to understand Figure 3-5 which shows the frequency domain for only greater than zero and moreover real and not complex.

Also why does the frequency domain of the impulse train have the original frequency spectrum appear flipped half the time? My current understanding is that the repetition comes from something like
$$\mathcal{F^{-1}} \{f(t)\Delta(t)\} = F(x)*\Delta(x)$$
where $\Delta$ is the dirac comb and probably some kind of constant like $2\pi$ is needed in there somewhere. How does the flipping half the time in the frequency domain arise from $F(x)*\Delta(x)$. The image I think for convolution is that $\Delta$ is flipped and as it is moved to the right by $s$ the "area" of $F(x)\Delta(s-x)$ is the value of the convolution at $s$. This would seem to replicate the original spectrum multiple times but without the flipping half the time as shown in Figure 3.5.

Let $X(f)$ denote the Fourier transform of $x(t)$ where
$$\begin{align*}
X(f) &= \int_{-\infty}^{\infty} x(t) \exp(-j2\pi ft) \mathrm dt\\
x(t) &= \int_{-\infty}^{\infty} X(f) \exp(+j2\pi ft) \mathrm df
\end{align*}$$
which I will denote via $x(t) \leftrightarrow X(f)$.
The following transform pairs will be needed in what follows.
$$\begin{align*}
\delta(t) &\leftrightarrow 1\\
\delta(t-t_0) &\leftrightarrow \exp(-j2\pi f \,t_0)\\
\sum_{n=-\infty}^{\infty}\delta(t-nT) &\leftrightarrow
\frac{1}{T}\sum_{k=-\infty}^{\infty}\delta\left(f-\frac{k}{T}\right)\\
\end{align*}$$

Given a signal $x(t)$, its sampled pulse train (at intervals of $T$ seconds) is
$$x(t) \sum_{n=-\infty}^{\infty}\delta(t-nT)
= \sum_{n=-\infty}^{\infty} x(t)\delta(t-nT)
= \sum_{n=-\infty}^{\infty} x(nT)\delta(t-nT).$$
Since multiplication in the time domain corresponds to convolution in the
frequency domain, we have
$$\begin{align*}
x(t)\sum_{n=-\infty}^{\infty}\delta(t-nT) &\leftrightarrow
X(f)\circledast\frac{1}{T}\sum_{k=-\infty}^{\infty}\delta\left(f-\frac{k}{T}\right)\\
&\leftrightarrow \frac{1}{T}\sum_{k=-\infty}^{\infty}
X(f)\circledast\delta\left(f-\frac{k}{T}\right)\\
&\leftrightarrow \frac{1}{T}\sum_{k=-\infty}^{\infty}
\int_{-\infty}^{\infty} X(f-w)\delta\left(w-\frac{k}{T}\right) \mathrm dw\\
x(t)\sum_{n=-\infty}^{\infty}\delta(t-nT)
&\leftrightarrow \frac{1}{T}\sum_{k=-\infty}^{\infty}
X\left(f-\frac{k}{T}\right) = \hat{X}(f)
\end{align*}$$
Thus, the Fourier transform of the impulse train formed by
sampling $x(t)$ at $T$ second intervals is $\hat{X}(f)$ which
is obtained by repeating $X(f)$ along the $f$ axis at intervals of
$T^{-1}$ Hz and summing the result.
Furthermore, $\hat{X}(f)$ is a periodic function of the frequency
variable $f$ with period $T^{-1}$ Hz.  That is, for all $f$,
$$\hat{X}\left(f + \frac{1}{T}\right) = \hat{X}(f).$$  Note that
all of this holds regardless of what $X(f)$ is: $X(f)$ could be
nonzero for all $f$, and the result would still be valid.

Now suppose that $X(f)$ is zero if $f < a$ or $f > b$.  Then,
$\hat{X}(f)$, obtained by repeating $X(f)$ periodically along
the $f$ axis,  is nonzero only in the intervals
$$\ldots, \left[a-\frac{2}{T}, b-\frac{2}{T}\right],
\left[a-\frac{1}{T}, b-\frac{1}{T}\right],
[a, b]
\left[a+\frac{1}{T}, b+\frac{1}{T}\right],
\left[a+\frac{2}{T}, b+\frac{2}{T}\right],
\ldots
$$
and so if
$$b -\frac{1}{T} < a \Rightarrow b-a < \frac{1}{T},$$
that is, the support $b-a$ of $X(f)$ is smaller than
the repetition interval $T^{-1}$ Hz, then the repetitions
of $X(f)$ do not overlap.  Indeed as the OP's figures show,
for a real-valued signal whose spectrum extends from $a = -f_0$
to $b = f_0$ (support of $X(f)$ is of length $2f_0$), sampling
$x(t)$ at intervals of $T = (3f_0)^{-1}$ (and thus repeating $X(f)$ at
intervals of $3f_0$ Hz on the frequency axis)
leads to no overlap while sampling $x(t)$ at intervals of
$T = (1.5f_0)^{-1}$ leads to overlap of the repetitions
of $X(f)$.

Finally, the OP asks about negative frequencies and the
interpretation of the pictures which show positive frequencies
only.  For a real-valued signal, $X(f)$ has conjugate symmetry,
meaning that $X(-f) = X^*(f)$, and so specifying $X(f)$ for
positive values of $f$ suffices.  In any case, the pictures
he is looking at are of $|X(f)|$ and $|\hat{X}(f)|$ which
are even functions of $f$, and so showing only the positive
axis saves space, though it does make the pictures look a
bit lopsided since only half the lobe is shown at low frequencies.
For complex-valued signals, $X(f)$ does not have conjugate
symmetry and $|X(f)|$ need not be an even function of $f$
and so the whole axis would need to be shown.  But, the general
development above is still applicable, and we still need
to sample at rate exceeding $b-a$ Hz, and it helps to keep in
mind that in this general case, each sample is actually two real
numbers, not one, since we are sampling a complex-valued signal.