Waveform samples and I/Q components I'm quite new to DSP and I guess my problem it's quite silly but I can't find a clear answer. Basically I'm writing some python code for a test. The test requires me to write a memory with the I/Q components of a specific waveform.

As far as I know it's quite easy in hardware to get the I/Q samples. I just need to multiply the input signal ( my waveform ) with $sin(w_ct)$ and $cos(w_ct)$ ($w_c$ should be the sampling frequency) and then sample it with an ADC.

In my case I'm able to create the samples of the waveform ( it's a sawtooth but I would like to understand the general case as well ) but then I'm quite stuck because I don't know how to extract the I/Q components from each sample.

Any idea ?

Cheers,
Stefano.

Assume that signal that comes to you through your antenna is:

$$ s_r[n] = I[n] \cos(\omega_0 n) + Q[n] \sin( \omega_0 n ) $$

which is a bandpass, modulated, narrow band signal where the components $I[n]$ and $Q[n]$ are the lowpass baseband in-phase and quadrature-phase components respectively.

Then in order to obtain the I and Q components from $s_r[n]$ you can apply the following just like the analog case:

$$I[n] = h_{LP}[n] \star ( s_r[n] \cos(\omega_0 n) ) $$
and
$$Q[n] = h_{LP}[n] \star ( s_r[n] \sin(\omega_0 n) ) $$

Where the lowpass filter $h_{LP}[n]$ is will have a cutoff frequency larger than the bandwidth of the signals I and Q and less than any adjacent band (if there is any). You may also apply a scaling such as $2x$ to adjust the magnitudes.

If you would like to form the analytic lowpass signal from the I and Q then you can easily form it like:

$$x_a[n] = I[n] + j Q[n]$$

where it can be shown that the analytic signal is also given by

$$ x_a[n] = A[n] e^{ j\phi(n) } $$

where $A[n] = \sqrt{ I[n]^2 + Q[n]^2 } $ and $\phi(n) = \tan^{-1}(\frac{Q[n]}{I[n]})$

Depending on the application, information in $A[n]$ (the envelope) or in $\phi(n)$ (the phase) is used.

Thanks to Fat32 I'm starting to understand ( I've realized my ignorance on this subject is huge ). I've spent some time to write some python code in order to understand. It helped to make me understand I hope it can be useful to others:

import sys
import os
import math
import argparse
import itertools
import subprocess

import numpy as np

def plot_samples(samples,name, unique = True):

fp = open(name+".gnuplot","w")
  fp.write("set grid\n")
  if unique:
    fp.write("plot \""+name+".dat\" smooth unique\n")
  else:
    fp.write("plot \""+name+".dat\"\n")
  fp.close()
  fp = open(name+".dat","w")
  for sample in samples:
    fp.write(str(sample[0])+" "+str(sample[1])+"\n")
  fp.close()
  gnuplot_pid = subprocess.call(["/usr/bin/gnuplot","-persist",name+".gnuplot"])

def sawtooth_wave(amplitude = 1.0, sampling_freq = 5100.0,frequency = 50.0, num_of_periods = 10, num_harmonics = 50):
  if frequency * num_harmonics * 2 > sampling_freq:
    print("Warning the sampling frequency is less than twice the most high harmonic. This will lead to distortions")
    print("Max signal frequency : "+str(frequency)+" x "+str(num_harmonics)+" = "+str(frequency * num_harmonics))
    print("Sampling frequency : "+str(sampling_freq))
    sys.exit()
  #A/2 - A/pi * sum_k(-1^k * sin(2*pi*k*f*t)/k)
  num_samples = math.ceil(float(sampling_freq)/float(frequency))
  cycle = num_of_periods * 1.0/float(frequency)
  time_axis = np.arange(0,cycle,1.0/float(sampling_freq))
  samples = []

  for t in time_axis:
    harmonics = 0
    for k in range(1,num_harmonics):
      harmonics += math.pow(-1,k)*np.sin(2*np.pi*k*float(frequency)*t)/k
    samples.append([t, -amplitude/np.pi * harmonics])
  return samples,len(time_axis)

def parser():
  parser = argparse.ArgumentParser(description='Shows some stuff')
  parser.add_argument('--frequency', default='50', type=int, help='Frequency of the tone in Hz (%(default)s)')
  parser.add_argument('--amplitude', default='1.0', type=float, help='Amplitude of the tone (%(default)s)')
  parser.add_argument('--sampling_freq', default='2400', type=float, help='Frequency of the sampling (%(default)s)')
  parser.add_argument('--carrier', default='2399', type=float, help='Frequency of the carrier (%(default)s)')
  parser.add_argument('--harmonic', default='20', type=int, help='Number of harmonics for the sawtooth (%(default)s)')
  return parser.parse_args()

def get_fft(samples,sampling_freq):
  samples_f = np.fft.fft([i[1] for i in samples])
  freqs = np.fft.fftfreq(len(samples)) * sampling_freq
  return [[f,np.absolute(i)/len(samples_f)] for i,f in zip(samples_f,freqs)]

def __main__():
  f2q1_8 = lambda v: (int((v * float(1 << 8)) + (-0.5 if (v < 0) else 0.5)) & 0x1FF)
  args = parser()
  samples,num_sample = sawtooth_wave(amplitude = args.amplitude, sampling_freq = args.sampling_freq, frequency = args.frequency, num_harmonics = args.harmonic)
  plot_samples(samples,"sawtooth")

  samples_f = get_fft(samples,args.sampling_freq)
  plot_samples(samples_f,"spectrum_sample")

  # I_Q components
  I_components = [[i[0],i[1] * np.cos(2*np.pi*args.carrier*i[0])] for i in samples]
  Q_components = [[i[0],i[1] * np.sin(2*np.pi*args.carrier*i[0]) * -1.0 ] for i in samples]
  plot_samples(I_components,"sawtooth_I")
  plot_samples(Q_components,"sawtooth_Q")
  I_components_f = get_fft(I_components,args.sampling_freq)
  plot_samples(I_components_f,"spectrum_sawtooth_I")
  Q_components_f = get_fft(Q_components,args.sampling_freq)
  plot_samples(Q_components_f,"spectrum_sawtooth_Q")

  # I_Q mixing
  sample_after = [[i[0], i[1] * np.cos(2*np.pi*args.carrier*i[0]) - q[1] * np.sin(2*np.pi*args.carrier*q[0])] for i,q in zip(I_components,Q_components)]
  plot_samples(sample_after,"sawtooth_after_decomposition")
__main__()

This piece of code basically create some samples of a sawtooth waveforms and print the I/Q waveforms/spectrum. At the moment the model miss the low pass filter and it's a concern because when I put together the I/Q i don't see any evident problem in the shape of the waveform. So I guess there is something I'm still missing.


IFFT-based OFDM system Please explain, mathematically, how can a bank of modulators be replaced with the IFFT in an OFDM system?

Say you want to transmit complex symbols $s_0,s_2,\ldots,s_{N-1}$ using frequency-division multiplexing. In the "bank of modulators" approach, you  would assign a carrier frequency $f_k$ to each symbol and create the signal
$$x(t) = \sum_{k=0}^{N-1} s_k \exp(j2\pi f_kt).$$

Contrast this with the definition of the inverse discrete Fourier transform:
$$x[n] = \sum_{k=0}^{N-1} s_i \exp(j2\pi nk/N).$$

So, if you interpret the symbols as the magnitude and phase of each carrier, and choose the carriers as $f_k=k/N$, then you'll see that you can generate $x(t)$ simply by converting $x[n]$ to an analog signal.

Note that:

I've explained how to transmit one block of symbols; extending this to multiple blocks has a slight complication and requires inserting a cyclic prefix to $x[n]$.
I'm assuming a rectangular pulse shape, which results in sinc-shaped carriers in frequency. The beauty of this method is that the carriers are orthogonal even if they overlap. In practice, this means that you can recover the symbols $s_k$ by performing the DFT of $x(t)$.
It is perfectly possible to use other pulse shapes with this method.

In order to understand the Pilots and Cyclic Prefix insertion you may want to have a look to a simpler OFDM blocks system. I hope this one: http://www.mathworks.com/matlabcentral/fileexchange/42902-lte-16qam in Simulink can help you. Let me know if you have any question about the model :)


Why does the Hilbert-Huang transform not suffer from the uncertainty principle, and what are the implications of this? I have been studying the Hilbert-Huang transform, which I understand is a frequency analysis technique that uses the adaptive IMFs basis components to decompose a signal using EMD.

I understand that frequency analysis techniques such as the STFT (short-time Fourier transform) and Wavelet transform involve the convolution of a signal with a predefined basis in order to extract frequency information.

Several articles from Huang and others state "convolution processes involve integration, which makes the results suffering the limitation imposed by the uncertainty principle".

http://web-static-aws.seas.harvard.edu/climate/pdf/Zhaohua.pdf
https://pyhht.readthedocs.io/en/latest/tutorials/hilbert_view_nonlinearity.html

I also understand that the frequency of each IMF is extracted by finding the time rate of change of phase of the analytic signal for each IMF (which would involve differentiation of the phase with respect to time).

Why does this differentiation mean that the uncertainty principle is avoided?

I have a vague idea that it's due to operators associated with time and frequency that don't commute, but I have struggled to find any sources that approach this from a signal processing/mathematical perceptive rather than in a quantum mechanical sense.

Furthermore, if this kind of signal processing circumvents the uncertainty principle, does that mean time and frequency events can be resolved with arbitrary precision?

(In my head, this, in turn, would have implications in quantum mechanics, but I imagine this would be outside the topic of this site.)

Any insights and/or articles are gratefully received.

I'm not that familiar with the Hilbert-Huang transform, so I won't comment on that.

You seem to be under the impression that the DFT "suffers from the uncertainty principle".  This is not true.  Don't feel bad about thinking that as it is a widespread misunderstanding.  This is a quote from an email to me from a well known expert in the field after I tried to get them to look at my exact frequency formula for a single (non-integer) pure real tone in a DFT:

  An "exact" computation of frequency based on FTs is not possible for
  some very simple reasons:  The Fourier Transform is itself an
  estimator, and the time-frequency uncertainty principle means
  frequency can't be determined exactly without an infinitely long
  observation window.  So essentially all frequency computations of time
  series or signals are estimates, it's just that some are more accurate
  or more efficient than others.

  So claiming that you have an exact solution just isn't going to get
  the attention of everyone interested in the field.

The first paragraph is entirely hogwash.  The last sentence was prophetic, hence I started my blog so no one could dispute the math.  You will find my original exact solution here:

Exact Frequency Formula for a Pure Real Tone in a DFT

I have discovered, and written about, several more.  You will find them among my articles.  There are also a few in the pipeline, including a correction to Macleod's three bin formula which makes it exact also.

The FT does apply to the uncertainty principle in physics, but in a different manner.  There they are talking about the transform of a Gaussian curve which happens to be an eigenfunction of the FT (but not the DFT, close though) and how it is the "narrowest function" in the transform.

I'm not an expert in Physics either, but I dabble.  I consider this to be my greatest mathematical discovery ever (exclipses the simple math of the frequency formulas by a mile):

Off Topic: Refraction in a Varying Medium

This is the vector form of Snell's law in a continuously varying index of refraction medium.  It is a theoretical particle model so its applicability to light is still a subject of study for me.  I believe it has the potential to illuminate subatomic behavior (i.e. the quantum world) but I haven't gotten there yet.


Adding Bass to a signal How can bass/low end be added to a signal?

With a low-pass filtered sub-oscillator running at an octave or two below the signal, there still doesnt seem to be enough bass.

Is there a mathematical function to easily add some warmth or boom to the sub oscillator? Or possibly adding two in unison, slightly detuned?

Let's say the audio is given by $ x $.
Let's define Low Pass Filter by $ h $.

Whet you need to do is:

$$ y = x + \alpha (h \ast x) $$

Where $ \alpha $ is the amplification level of the Low Pass data.
Once you set which frequencies you want to amplify it is easy to design.

You could use a low-pass filter and add this to the original signal to boost the bass. What kind of speakers are you using? Regardless of how much you boost the bass signal, you need to have speakers that can put out the bass you want.


How to understand Kalman gain intuitively? The Kalman filter algorithm works as follows

  Initialize $ \hat{\textbf{x}}_{0|0}$ and $\textbf{P}_{0|0}$.

  At each iteration $k=1,\dots,n$

  Predict

  Predicted (a priori) state estimate    $$ \hat{\textbf{x}}_{k|k-1} =
 \textbf{F}_{k}\hat{\textbf{x}}_{k-1|k-1} + \textbf{B}_{k}
 \textbf{u}_{k} $$
    Predicted (a priori) estimate covariance $$
 \textbf{P}_{k|k-1} = \textbf{F}_{k} \textbf{P}_{k-1|k-1}
 \textbf{F}_{k}^{\text{T}} + \textbf{Q}_{k}$$
  Update

  Innovation or measurement residual     $$ \tilde{\textbf{y}}_k =
 \textbf{z}_k - \textbf{H}_k\hat{\textbf{x}}_{k|k-1}$$ Innovation (or
  residual) covariance   $$\textbf{S}_k = \textbf{H}_k
 \textbf{P}_{k|k-1} \textbf{H}_k^\text{T} + \textbf{R}_k$$ Optimal
  Kalman gain      $$\textbf{K}_k =
 \textbf{P}_{k|k-1}\textbf{H}_k^\text{T}\textbf{S}_k^{-1}$$
    Updated (a posteriori) state estimate      $$\hat{\textbf{x}}_{k|k} =
 \hat{\textbf{x}}_{k|k-1} + \textbf{K}_k\tilde{\textbf{y}}_k$$
    Updated (a posteriori) estimate covariance     $$\textbf{P}_{k|k} = (I -
 \textbf{K}_k \textbf{H}_k) \textbf{P}_{k|k-1}$$

The Kalman gain $K_k$ represents the relative importance of the error $\tilde{\textbf{y}}_k$ with respect to the prior estimate $\hat{\textbf{x}}_{k|k-1}$.

I wonder how to understand the formula for the Kalman gain $K_k$ intuitively? Consider the case when the states and outputs being scalar, why is the gain bigger, when

$\textbf{P}_{k|k-1}$ is bigger
$\textbf{H}_k$ is bigger
$\textbf{S}_k$ is smaller?

Thanks and regards!

I found a good way of thinking intuitively of Kalman Gain $K$. If you write $K$ this way

$\displaystyle \quad\ \bf{K_k} = \bf{P_k^-\, H_k^{\rm T} (H_k P_k^-\, H_k^{\rm T} + R_k)^{-1}}
  = \bf{\frac {P_k^-\, H_k^{\rm T}}{H_k P_k^-\, H_k^{\rm T} + R_k}}$

you will realize that the relative magnitudes of matrices ($R_k$) and ($P_k$) control a relation between the filter's use of predicted state estimate ($x_{k}⁻$) and measurement ($ỹ_k$).

$\displaystyle \quad\
  \lim\limits_{\bf{R_k \to 0}} \bf{{P_k^-\, H_k^{\rm T}} \over\
  {H_k P_k^-\, H_k^{\rm T} + R_k}}\
    = \bf{H_k^{-1}}$

$\displaystyle \quad\
  \lim\limits_{\bf{P_k \to 0}} \bf{{P_k^-\, H_k^{\rm T}} \over\
  {H_k P_k^-\, H_k^{\rm T} + R_k}}\
   = \bf 0$

Substituting the first limit into the measurement update equation

$\displaystyle \quad\
  \bf{\hat x_k} = \bf{x_k^-} + \bf{K_k}(\bf{\tilde y_k}-\bf{H_k}\bf{x_k^-})$

suggests that when the magnitude of $R$ is small, meaning that the measurements are accurate, the state estimate depends mostly on the measurements.

When the state is known accurately, then $H P^⁻ H^T$ is small compared to $R$, and the filter mostly ignores the measurements relying instead on the prediction derived from the previous state ($x_k⁻$).

I am working on the Kalman Filter (KF) algorithm. I observed that the kalman gain deals with convergence of algorithm with time, that is, how fast the algorithm corrects and minimizes the residual.

Coming to the equation choose an initial kalman gain value and vary it from low to high, that can give you an approximated one.


Calculating the energy of trapezoid wave Here is the trapezoid wave,although i know we can use $\int^{4}_0 |x(t)|^2dt$ to calculate its power,but i want to divide this trapezoid to right triangle$(t=0$~$2)$,and a square ($(t=2$~$4)$).

I know the the triangle effective voltage is $V_{eff}=\frac{V}{\sqrt{3}}$,and i want to use $P=\frac{V_{eff}^2}{R}$ to calculate the power ,  but i am not sure about which  triangle is correct,i mean,for which triangle(a or b or c),i can use $\frac{V}{\sqrt{3}}$ to calculate its power,or are these three wave's power the same,but different in energy?

it seems that there power are the same,but there energy are not the same ,because energy=power $\times$ time,and time are different,so their energy are not the same too


Signals - question due to sampling I don't have an idea how to start it. Help pls :)

$$m = z(t) + \frac{\operatorname{Sa}(t-T)}{4}$$

We have some signal $z(t)=\left(\operatorname{Sa}(\frac{t}{T})\right)^2$. How often the m signal should be sampled to be sure of correct signal reconstruction

This is not an answer.  I am just trying to restate the question so that the equation is dimensionally homogeneous:

I don't have an idea how to start it. Help pls :)

$$m(t) = z(t) + \tfrac{1}{4}\operatorname{sinc}\left(\tfrac{t-T}{T}\right)$$

We have some signal $z(t)=\left(\operatorname{sinc}(\tfrac{t}{T})\right)^2$. How often should the $m(t)$ signal be sampled to insure correct signal reconstruction?

this question can be meaningfully answered.

let's define:

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

$$ \Pi(u) \triangleq \begin{cases}
1 \qquad & |u| < \tfrac12 \\
\tfrac12 \qquad & |u| = \tfrac12 \\
0 \qquad & |u| > \tfrac12 \\
\end{cases} $$

$$ \Lambda(u) \triangleq \begin{cases}
1 - |u| \qquad & |u| \le 1 \\
0 \qquad & |u| > 1 \\
\end{cases} $$

So in order to reconstruct a signal without loosing information you have to sample with $\frac{1}{T}=f_\mathrm{s}>2B$.  This is known as the Nyquist rate.

Therefore in order to determine the minimal sampling rate we need to examine the function in the frequency domain.

$m(t) = \operatorname{sinc}\left(\tfrac{t}{T}\right)^2 +\dfrac{\operatorname{sinc}({t-T})}{4}$

So
The Fourier transform would be

$M(f) = (T\operatorname{\Pi}(T\pi f) \circledast T\operatorname{\Pi}(T\pi f)) +  \dfrac14\operatorname{\Pi}(\pi f)e^{-j2\pi fT}=T\Lambda (\pi f T) +  \dfrac14\operatorname{\Pi}(\pi f)e^{-j2\pi fT}$

in order to faithfully reconstruct the signal :
$ f_\mathrm{s}=\max(4T,2)$

The exponent doesn't matter it's bounded by the $\operatorname{\Pi}(\cdot)$


Regarding the choice of cost function in adaptive control - squared error vs absolute error I did search the question database regarding this question, and although one or two questions came close, they didn't really address my specific question.

In adaptive control based on minimizing tracking error (e.g. between plant and model), the designer is free to chose a cost function. More often than not the cost is selected as a function of the squared error.

But I've found in some practical applications that I can achieve a more robust controller by using absolute error. I understand that absolute error provides a more uniform weighting on the size of the error, and I suspect the squared error tends to 'wind-up' the adaptive controller with initially large errors. But I'm not sure how to show this in a generalized way. So I have two questions regarding this:

Is there perhaps a simple analysis that can demonstrate the stability characteristics between absolute and squared error choices in the cost function?
Any references on the matter?

This is an interesting question since both squared error and absolute error are convex functions, so they are both going to give the optimal solution when minimized. My intuition is that the $\ell_2$-norm (sum of squared values of error) converges to zero more quickly than the $\ell_1$-norm (sum of absolute values of error) when the search direction is right. For example, if the error is $\epsilon = y - Ax$. If $y$ is very close to $Ax$ and $0 < |y-Ax| \leq 1$, then $\epsilon^2$ is much smaller than $|\epsilon|$. Similarly, if $|y-Ax| > 1$, then $\epsilon^2 > |\epsilon|$. So the $\ell_2$ norm penalizes large errors more, and small errors less than the $\ell_1$ norm. This can also be understood by looking at the graphs of the two functions when $\epsilon \in \mathbb{R}$.

Of course the continuity of the derivative is another factor. The derivative of the absolute value technically does not exist at the minimum.

You should use squared error as versus absolute error if you expect to have both positive and negative error values in your population.  By squaring the error, only the magnitude of each error value matters (although it is squared).  If you want to recover the actual magnitude of error in the population, just take the square root at the end.  This value is known as the root mean squared error (RMSE).

If, instead, you were to use absolute error values, the positives and negatives in your population would cancel each other out giving you a false metric.


Unexpected peaks in spectral analysis I am doing the spectral analysis of some antenna recorded signals by computing PSD (logarithmic absolute value of the FFT) and the spectrogram of the signal. The problem is that I get some unexpected spikes in the signal that appear through all the signal, as it can be seen in the next figure:

The two peaks that appear are the spikes I am speaking about, and they are the same in all the signal. The center spike is not located in $0\textrm{ Hz}$, so it is not a DC offset of the signal. This is giving me some problems, so I would like to know what can they be or how can they be removed (maybe I am not computing the spectral analysis correctly).

I have also performed Welch PSD estimations, and the peaks still appear. I have been told that it can be some internal interference form the electronics recording the data or some issue with ADCs.

The large spike at ~9MHz could be some form of interference from electronics near the receiver or internal to the receiver (maybe from using a poorly conditioned power source).  My best guess for the spike at DC is a DC bias from quantization.  If the digitizer truncates instead rounds, the resulting digital signal will have a non-zero mean and therefore a DC offset.  We had this happen in an FPGA downcoverter implementation at a previous job.  Once they changed the FPGA to round, the DC offset went away.  If you have ever used an Ettus USRP radio, you see the same DC offset in many of those boxes.


FFTW - Inverse FFT input data I am trying to use FFTW function for IFFT:

fftw_plan_dft_c2r_1d(size, out_cpx, out, FFTW_ESTIMATE).

My goal is to input desired frequency response data and obtain related coefficients using Inverse Fourier Transform. Since out_cpx is complex, I calculated magnitude (0,1) and phase, taking into account it is linear phase filter. Since I know the phase I calculated real and imaginary parts as Acos(w) and Asin(w) and with these pairs I populated out_cpx[i][0] and out_cpx[i][1]. When I perform Inverse Fourier Transform using fftw_plan_dft_c2r_1d() I assume I am right away getting FIR coefficients h[n] as the result. With these coefficients i go to MATLAB and get response using freqz(C,1) where C holds my coefficients. But I am getting a lot of ripples in my graph for response and for phase. Am I doing something wrong here? It is not clear to me how to specify phase of the desired frequency response.

The ripple you see is due to Gibbs phenomenon and is caused by your filter having a step response. You will have to pass the frequency response through a window function first.


How detect and filter noise from signal This is the spectrum of the signal:

And this is the signal:

On the FFT spectrum you can notice the high amplitude peaks and I think these frequencies are the main ones in the signal:

I think that the "low" harmonic must be clean, but I have the problem that low amplitude components are placed in every frequency interval.  The main problem is that I haven't got any model of the signal for a hypothesis about noise, and as result, I must detect noise components.  After this step, I must design a filter or use another method of filtering for denoising the signal.  Which methods can I apply for (1) detecting noise and (2) filtering?

Just "eyeballing" your signal, it looks like the interesting things are all down below 1Hz.  The peaks in your signal browser diagramm seem to be spaced at every 5 to 10 Seconds, which would be 0.1 to 0.2 Hz.  Your spectrum plot goes from zero to 25 Hz, and doesn't show anything really interesting.

Try zooming in on the spectrum and see what lies below 1 Hz.  Yo may well find that the signal itself lives down in the very low frequencies, and that separating the noise from the signal is as simple as using a low pass filter.

You might try to window the FFT so that you do not introduce harmonics. Then filtering in the FD might be easier.


What is the physical meaning of the convolution of two signals? If we convolve 2 signals we get a third signal. What does this third signal represent in relation to the input signals?

There's not particularly any "physical" meaning to the convolution operation. The main use of convolution in engineering is in describing the output of a linear, time-invariant (LTI) system. The input-output behavior of an LTI system can be characterized via its impulse response, and the output of an LTI system for any input signal $x(t)$ can be expressed as the convolution of the input signal with the system's impulse response.

Namely, if the signal $x(t)$ is applied to an LTI system with impulse response $h(t)$, then the output signal is:

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

Like I said, there's not much of a physical interpretation, but you can think of a convolution qualitatively as "smearing" the energy present in $x(t)$ out in time in some way, dependent upon the shape of the impulse response $h(t)$. At an engineering level (rigorous mathematicians wouldn't approve), you can get some insight by looking more closely at the structure of the integrand itself. You can think of the output $y(t)$ as the sum of an infinite number of copies of the impulse response, each shifted by a slightly different time delay ($\tau$) and scaled according to the value of the input signal at the value of $t$ that corresponds to the delay: $x(\tau)$.

This sort of interpretation is similar to taking discrete-time convolution (discussed in Atul Ingle's answer) to a limit of an infinitesimally-short sample period, which again isn't fully mathematically sound, but makes for a decently intuitive way to visualize the action for a continuous-time system.

In simple terms it means to transfer inputs from one domain to another domain where we find it easier to work with. Convulation is tied with Laplace transform, and sometimes it is easier to work in the s domain, where we can do basic additions to the frequencies. and also as laplace transform is a one to one function we are most likely to not corrupt the input.
Before trying to understand what the general theorem of convulation means in physical significance, we should instead start at the frequency domain. addition and scalar multiplication follows the same rule as Laplace transform is a linear operator. c1.Lap(f(x)+ c2.Lap g(x)= Lap (c1.f(x) + c2.g(x)).
But what is Lap f(x).Lap g(x). is what defines the convulation theorem.


How to find frequency response, stability, and causality of a linear system? I have the following transfer function:

$$H(s)=\frac{s}{(s+1)(s+2)}$$

How can I find the gain and phase response of the above system? I know the first step has something to do with substituting $s = j\omega$ into $H(s)$. How can I find whether the system is stable and whether it is causal?

NOTE: This answer is community wiki: If you feel there is a need for
  correction -please update the answer right here.

1. Frequency response:
For any system - the (discrete time or continuous time) Fourier Transform of the impulse response is same as the frequency response of it. You can replace $s$ by $j\omega$ to work on this.

2. Stability:
System is stable if, the bounded input produces bounded output. The definition of can be found here in the wiki link.  For this, the necessary condition to prove is to see if the impulse response be absolutely integrable, i.e., its L^1 norm exist.

For a rational and continuous-time system, the condition for stability is that the region of convergence (ROC) of the Laplace transform includes the imaginary axis.

3. Causality:
Any system who's impulse response at sample $s_i$ doesn't require to know any samples presented $s_{i+1}$ or after wards, than the system is causal.

EDIT:
Now given the above i am trying to put the calculation as per your system.

$$ H(s)=\frac{s}{(s+1)(s+2)}  $$

Now we can replace $s$ with $j\omega$ to get the frequency response in terms of magnitude and phase response.

$$ H(j\omega) =  \frac{j\omega}{(j\omega+1)(j\omega+2)} $$
$$ H(j\omega) =  \frac{j\omega}{2 - \omega^2 + 3j\omega } $$
$$ Magnitude =  20 Log_{10} \sqrt { ( H(j\omega) )^2 } $$
$$ Magnitude =  20 Log_{10}  \sqrt { \frac { (j\omega)^2 }{(2 - \omega^2)^2 + (3j\omega)^2} }  $$

NOTE: square is applied separately on real terms than on imaginary terms! I got this understanding from this document here.

Replacing $-j^2$ by $-1$

$$ Magnitude =  20 Log_{10} \sqrt { \frac{-1\cdot\omega^2}{(2 - \omega^2)^2 - (3\omega)^2} } $$

$$ Magnitude =  20 Log_{10} \sqrt { \frac{-1\cdot\omega^2}{(4 - 5\omega^2 +\omega^4)} } $$

$$ Phase = arc tan (real / imaginary)  $$

$$ Phase = -1 \cdot \tan^{-1}(\frac {-3\omega}{(2-\omega^2)} )   \ $$

Please correct me if i am wrong.

Another literature that can help you.


Is there a way to separate 2 signals recorded on the same channel? I have 2 low frequency signals on the same line that look like this:

I would like to separate them so that I obtain the two signals that can be seen here:

The sampling is not perfectly intercalated, so it's not as simple as taking one point each in alternation.

I looked at ICA but my understanding is that I would need 2 recordings of the 2 mixed signals in order to separate them.

Any ideas?

You should look into a generalization of ICA called Independent Subspace Analysis (ISA). Practically it is Principal Component Analysis (PCA) followed by ICA.

There are threads on dsp.stackexchange that discuss this as well as github repos (and here) that provide various implementations.

It would appear that the two signals are multiplexed (naive time-division multiplexing), that is to say, even-numbered samples belong to signal 1 and odd-numbered samples to signal 2 - in which case the separation is trivial. Is that the case?


Using large DFT for OFDM to avoid the guard interval In OFDM, we usually insert the guard interval to avoid the inter-symbol interference at the expense of achieved data rate. that guard interval can be either zeros or CP guard interval, which it must be longer than the channel delay.
I am wondering, what's about using long symbol to avoid that guard interval ? for example instead of having number of subcarriers $N= 64$ or $N = 1024$, we can use $N = 5120$, what's the impact of using a high number of $N$ ??

Sure,  you can !
Practically, some systems are using OFDM system with number of subcarriers 4096, but the PAPR will be an issue in that case, and if your channel is fast-fading, the time-domain equalizer will be high-complex since one-tap frequency-domain equalizer will be useless.
Usually, choosing the length of the symbols depends on the channel delay, e.g., if your channel delay is 128,  you must use at least length of symbol 1024 and so on.

No, that can't work.
The guard interval has the purpose of "swallowing" the previous' symbol's impulse response.
Your longer DFT distributes energy all over the symbol duration. This doesn't solve the problem.


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

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


Should I use CV_HAAR_SCALE_IMAGE while using LBP CascadeClassifier? I trained a cascade-LBP to detect Lollipops with 1000 images, now I'm trying to "adapt" the openCV HAAR-Cascade example to use my LBP trained .xml but I'm not sure about the "HAAR_SCALE_IMAGE" flag, I read that optimizes the downscale to increase the HAAR features is it right to use it with LBP?

//-- Detect Lollipops
  lollipops_cascade.detectMultiScale( frame_gray, lollipops, 1.1, 2, 0|CV_HAAR_SCALE_IMAGE, Size(30, 30) );

An example of imput image:

For more info:
http://docs.opencv.org/2.4/doc/tutorials/objdetect/cascade_classifier/cascade_classifier.html

The OpenCV documentation suggests using CASCADE_SCALE_IMAGE when performing a detectMultiScale with LBP features.

lollipops_cascade.detectMultiScale(frame_gray, lollipops,1.1, 2, 0|CASCADE_SCALE_IMAGE, Size(30, 30));

Even though, it works perfectly with  CV_HAAR_SCALE_IMAGE flag instead of CASCADE_SCALE_IMAGE.

According to EMGU, which is an .NET wrapper for OpenCV, and sometimes has way better documentation than OpenCV,

  DO_CANNY_PRUNING
  If it is set, the function uses Canny edge detector
  to reject some image regions that contain too few or too much edges
  and thus can not contain the searched object. The particular threshold
  values are tuned for face detection and in this case the pruning
  speeds up the processing

  SCALE_IMAGE
  For each scale factor used the function will downscale
  the image rather than "zoom" the feature coordinates in the classifier
  cascade. Currently, the option can only be used alone, i.e. the flag
  can not be set together with the others

  FIND_BIGGEST_OBJECT
  If it is set, the function finds the largest
  object (if any) in the image. That is, the output sequence will
  contain one (or zero) element(s)

  DO_ROUGH_SEARCH
  It should be used only when
  CV_HAAR_FIND_BIGGEST_OBJECT is set and min_neighbors > 0. If the flag
  is set, the function does not look for candidates of a smaller size as
  soon as it has found the object (with enough neighbor candidates) at
  the current scale. Typically, when min_neighbors is fixed, the mode
  yields less accurate (a bit larger) object rectangle than the regular
  single-object mode (flags=CV_HAAR_FIND_BIGGEST_OBJECT), but it is much
  faster, up to an order of magnitude. A greater value of min_neighbors may be specified to improve the accuracy.

Source


MATLAB : How to filter a data signal I have a data signal (column vector) with N values taken over at specific sample rate (dx) at a fixed frequency. After taking the Fourrier transform of this signal and representing it in the Fourier Domain (with fftshift), I need to eliminate the parts situated in the origin and in the positive side. The final signal I need is the one with the negative shift. Therefore, after filtering the unwanted parts I will need to perform an inverse Fourier transform to recover the signal in time domain and center it in the origin.
Finally, I do not understand how to filter the spectrum when there is no frequency variation and when I’ve tried specifying a sample frequency and then applying a matlab filter I didn’t get any results. Can anyone help me to understand how to make this work please?
PS: the spectrum is in the figure below

This task its not so hard, you only should create a vector, using ones, of the same size and set to zero the positions that correspond with the frequencies that you don't want. Apply an element-element multiplication between the vector and your signal. Finally apply the inverse Fourier transform in the obtained vector. The response you will obtain will be a complex one, so you need to graphic only real or imaginary part if you want to appreciate it in a graphical way, you also can graphic the absolute of the signal in order to observe the magnitude of it.

This type of filter is known as ideal filter and is impossible to apply in real world. In Fourier domain, filter is just a multiplication while in space domain to filter is necessary to apply a convolution between the signal and a kernel.

In general, removing the energy at DC and all the positive frequencies is equivalent to subtracting the mean and dividing by two.  The FFT will always be two-sided and symmetrical.


Obtaining power spectrum from ACF, FFT using Matlab and FFTW I am using Matlab R2012b 64-bit on Windows 7 in order to estimate the power spectrum of a simple signal that is:

$$\cos(10t) + \sin(20t) $$defined in the time interval from 0.0 to 10.0

Here is what I did:

1) There was no predefined value of points so I took 1001 points by forming a time array and sampled the function at these points.

2) I calculated the autocorrelation function of the signal with the built-in function "xcorr".

3) I performed FFT directly on the autocorrelation function with no additional options.

4) The resulting data set had imaginary parts(I think this should not be the case since the autocorrelation function is even).

5) I took the absolute value of the FFT of ACF and plotted it as the power spectrum from the ACF after having it shifted with fftshift.

6) I took the FFT of the sinusoidal signal after padding it with 1000 (n-1) zeros.

7) I calculated the power spectrum from the FFT of the signal as the square of the absolute value of the frequency coefficients.

8) I plotted both of them and compared them by assigning an array to difference between them.

9) They turned out to be the same.

Everything seems to be working correctly but my concerns are as follows:

How can this data be interpreted? I expected the sum and difference frequencies of the individual sinusoidal functions in the signal.

-How can I normalize the output and represent negative frequencies properly since they have no physical significance?

-Does Matlab use FFTW libraries? My main task is to achieve the same results in a C program using these libraries.

I have benefitted from the following URL:
Efficiently calculating autocorrelation using FFTs in organizing the post and achieving these results.

I am also confused about the terms power spectrum, power spectral density and energy spectrum. Are they used interchangeably? I would be grateful for an explanation of this phenomenon and normalization procedure.

Thanks in advance

P.S : The code snippet I have used is as follows and it can reproduce my results:

t=transpose(0:0.01:10);
n=size(t,1);
f=cos(10*t)+sin(20*t);
acf=xcorr(f,f);

FFTacf=fft(acf);
PSacf=abs(FFTacf);

FFTp=fft([f;zeros(n-1,1)]);
PSf=FFTp.*conj(FFTp);

figure;
subplot(2,1,1);
plot(fftshift(PSf));
title('PS from FFT');
subplot(2,1,2);
plot(fftshift(PSacf));
title('PS from ACF');

d=PSacf-PSf;
fprintf('Max error = %6.2f \n', max(abs(d)));

I am also confused about the terms power spectrum, power spectral density and energy spectrum. Are they used interchangeably? I would be grateful for an explanation of this phenomenon and normalization procedure.

EDIT with @jojek`s helpful comment:

power spectrum is NOT the same as power spectral density despite the fact that terms are often mixed up. Units of power spectral density (PSD) are /Hz, e.g. /Hz. Power Spectrum (PS) is equal to PSD*(Equivalent Noise Bandwidth of used window funktion) and has units of  . You often normalise it with N, the number of data points (for discrete time series). For physical processes the PSD and PS are symmetric, thus negative frequencies contain no additional information and only positive frequencies are plotted, negative frequencies are discarded.

energy spectrum is a term of spectroscopy: it describes the intensity of a particle beam as a function of particle energy. What you maybe meant is called energy spectral density, which describes the distribution of energy over frequency; it has units of J/Hz

A good overview is here.

I want to add the term amplitude spectrum/amplitude spectral density, which is the square root of the power spectral density.


Normalization factor in the convolution theorem Maybe it's a trivial question, but I couldn't find any explanation for this.
According to the convolution theorem, in the continues case we add normalization factor, i.e.
$$
\mathcal F\left\{h\star g\right\} = \sqrt{2\pi}\cdot \mathcal F\left\{h\right\}\cdot \mathcal F\left\{g\right\}
$$
Do we still need to add a normalization factor also in the discrete case? If we do, what should its value be?

The choice of the normalization factor is just a matter of convention. Note that the specific correspondence between convolution in the time domain and multiplication in the frequency domain with a scaling of $\sqrt{2\pi}$, as shown in your question, applies only to the unitary definition of the Fourier transform with angular frequency as the independent variable in the frequency domain:
$$F(\omega)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(t)e^{-j\omega t}d\omega\tag{1}$$
If we use the non-unitary definition (as is common in the field of signal processing)
$$F(\omega)=\int_{-\infty}^{\infty}f(t)e^{-j\omega t}d\omega\tag{2}$$
then the correspondence is simply
$$\mathcal{F}\{(f\star g)(t)\}=F(\omega)G(\omega)\tag{3}$$
The same is true if we use frequency instead of angular frequency.
For the discrete-time Fourier transform (DTFT), we have exactly the same choices.
For the discrete Fourier transform (DFT), where we deal with samples in both domains, the same holds for the normalization factor $N$, where $N$ is the DFT length. In the field of signal processing, the DFT is commonly defined as a non-unitary transform:
$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{4}$$
The inverse transform is then given by
$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{5}$$
The unitary DFT definition is
$$X[k]=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{6}$$
and the corresponding inverse transform is
$$x[n]=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{7}$$

There is no normalization factor in either the continuous or the discrete case. Please check your source. A normalization factor ONLY occurs for the Fourier Series with $2\pi$ periodic functions and circular convolution.
https://en.wikipedia.org/wiki/Convolution_theorem


error variance of frequency and phase estimation using DFT What is the error variance of estimating the (frequency) and (initial phase) of a sinusoid using DFT.
please provide references if possible.

Since I don't have sufficient reputation to comment i will quote and add to Dsp_user's answer.

Summary: The error variance will depend how well you can estimate the frequency of your signal. This in turn will be limited by the amplitude of the DFT peak vs noise level. Spectral leakage will likely also skew your estimates.

  the latter can't be assigned a physical meaning (please someone correct me if I'm wrong on this point).

this is incorrect, calculated phase from the DFT is linearly correlated with phase shift of a sinusoid. When generating a cosinus curve and shift it between 0:2*pi i get errors in the range of (10^-15) between signal phase shift and phase estimate. Phase estimates from DFT can thus be said to compare a signals shift to a non-shifted signal. Simulated curve is shown bellow:

In my opinion, phase estimation error should only be dependent on how well the "sinusoidal-fit" from the DFT fits the actual signal. One of those limitations (frequency resolution/spectral leakage) was explained by Dsp_user. If your signal is between frequency bins you will estimate the wrong signal amplitude (correlated with how close the frequency is to the middle between two bins). Since phase estimation by DFT is linerly dependent on signal frequency i would guess that you will introduce a higher/lower derivative when trying to estimate real time shifts between signals but I admit here that I am not 100%.

Added: In case of spectral leakage, if you have other sinusoids close to your peak of interest, their leakage could pour onto that and introduce errors.

Your noise floor (if it is white) will in real cases likely be your largest source of error. It will cause additive errors to your phase estimates whose size will depend on DFT-Peak Amplitude vs noise floor amplitude (noise power spread out over the signal length). Bellow you can find a figure of this as well as the code to generate both in my answer

t=0:1000^-1:.999;
phi=0:1000^-1:2*pi-1000^-1;
signal = cos(2*pi*10*t+phi');

FFTSig= fft(signal,[],2);
figure; plot(phi,squeeze(atan2(imag(FFTSig(:,11)),real(FFTSig(:,11)))),'linewidth',2)
xlim([0 2*pi])
ylim([-pi pi])
xlabel('Signal Phase Shift')
ylabel('Phase shift estimate from DFT')
title('Phase Shift Estimate - No Noise')

signal2 = cos(2*pi*10*t+phi');
for iter=1:6283
   signal2(iter,:)= squeeze(signal2(iter,:))+2*randn(length(signal(1,:)),1)';
end

FFTSig2= fft(signal2,[],2);
figure; plot(phi,squeeze(atan2(imag(FFTSig2(:,11)),real(FFTSig2(:,11)))))
xlim([0 2*pi])
ylim([-pi pi])
xlabel('Signal Phase Shift')
ylabel('Phase shift estimate from DFT')
title('Phase Shift Estimate - Additive White Noise')

I can't find any references but I'll try to give an intuitive explanation for these two points. First, while a DFT operation
produces both magnitude and phase information, the latter can't be assigned a physical meaning (please someone correct me if I'm wrong on this point).
The reason for this is obvious--phase has a physical meaning only when comparing 2 or more sinusoids etc but you'll get phase information even for
a single sinusoid so it's probably best to view phase as an inherent property of DFT (which is necessary to fully reconstruct the original signal). EDIT: Even when dealing with a single sinuosoid, DFT will still compare it to sinusoid basis functions, so we can think of phase (and magnitude) information as a match between the measured signal and the individual basis functions.

So, let's just focus on the frequency part of your question. The frequency resolution (bandwidth) can be computed as

bin frequency resolution = sampling rate / frame length;

This bin frequency resolution also determines the maximum error or variance, which is simply computed as half the frequency resolution.

frequency error/variance = bin frequency resolution / 2;

As an example, let's say that our signal is sampled at 8000 Hz and we use frames 512 in length for our DFT. We get

bin frequency resolution = 8000 / 512 = 15.625 Hz

and

frequency error/variance = 15.625 /2 = 7.8125 Hz

This is just the maximum frequency error/variance for a particular bin but the actual magnitude that is present in the bin can even be the result of contributions from adjacent bins (due to the
spectral leakage so the true frequency may not even be in this bin--for this see http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/)

UPDATE

  The more off-center it is in frequency from the bin frequency, the
  larger this phase offset will be. So, what we see here is that while
  signals whose frequency is exactly on a bin frequency have the same
  phase offset in each frame, the phase of signals that have a frequency
  between two bin frequencies will have an offset that is different with
  each frame. So, we can deduce that a difference in phase offset
  between two frames denotes a deviation in frequency from our bin
  frequencies.

This quote is from the article found in the link given above (the context is finding the true frequency in STFT).


Can an Delta Dirac function (not response) ever be sampled? This is purely theoretical question but can an actual delta-dirac function ever be sampled? To make it clear again, I'm not talking about an impulse response from the impulse function. Lets say we have a truly instantaneous impulse as defined with delta-Dirac that the impulse with infinite(or finite) height and no width, an instantaneous pulse. This pulse is to be detected by sampling, but to find this instant signal, wouldn't this require an infinitely fast sampling or infinite sampling rate? Thinking this in terms of frequency domain, I believe delta function is 1 or has infinite spectrum. So would delta function ever be sampled unless we have infinite sampling rate?

A Dirac delta impulse cannot be sampled because it is a distribution and as such it does not make sense to talk about the value of $\delta(t)$ for e.g. $t=0$ because that value isn't defined. However, sampling extracts values of a function for certain values of $t$. The Dirac impulse is no function in the conventional sense. It's only under an integral that $\delta(t)$ makes sense.

Furthermore, the mathematical description of sampling is a multiplication with a Dirac comb. However, the expression $\delta^2(t)$ is undefined, so sampling a Dirac impulse can't be described mathematically in any meaningful way.

Note that the discrete-time unit impulse $\delta[n]$ can be sampled without any problem, because it has well-defined values for any $n$, unlike its continuous-time counterpart.

Yeah, I think your conclusion is right. According to Nyquist' sampling theorem, the sampling rate should be at least twice as many as the maximum frequency of the target signal. And the maximum frequency of delta function is infinite. If sampled, aliasing must occur. In fact, the discrete delta function $\delta(n)$ is a sampled version of continuous delta function $\delta(t)$. The spectrum is overlapped while sampling, of course.


Voice recovery with high-noise, low-gain I have a recorded interview with two speakers, in which one voice's gain was much higher than the other (by three orders of magnitude). Increasing the gain does not help, because speaker 2's voice is so quiet that it is hitting the digitization limit (I believe). Audio quality is not important, but it would be very helpful to have an approximate technique so that the words may be intelligible.

If helpful, I can try to post an audio example of the specific distortion here. The cadence of the voice is there, but it's so mangled the words are unintelligible. Basic interpolation and filtering (both long shots) were not successful in recovering the voice.

Any suggestions for analysis techniques and/or software would be much appreciated.

I would suggest dynamic range compression as a first port of call. Set the threshold low and the ratio quite high. Then apply make-up gain to an appropriate level. This functionality is built into many free audio editors such as Audacity. Matlab or Python implementation is also reasonably simple.

I might try Audiophonic, a web-based podcast post-production service. It does exactly the type of processing you're wanting automatically. It's free for the first two hours a month of audio (as of this writing). It includes an "intelligent leveler" and "loudness normalization". If the volume difference is too extreme, a fully automated solution might choke. But it's worth a shot. Or you could do some of the tricks mentioned in other comments to get it in better shape and let Audiophonic do the rest.


QAM Amplitude Levels I'm working on a project to stream video between 2fpgas by using QAM. While mapping symbols with different amplitude pulses(PAM) a question popped up in my mind. In real life how does communicating wireless devices know which amplitude levels to agree on for digitally modulating-demodulating bitstreams? Is there a standard for that?

n real life how does communicating wireless devices know which amplitude levels to agree on for digitally modulating-demodulating bitstreams?

You can't "agree" on levels, since the received amplitude depends on the attenuation on the channel, and that depends on distance and a lot of other factors.
So, your receiver depends on detecting something that it knows the right amplitude (or power) of – for example, a preamble, or an average squared symbol, or something like that.
How that is done depends on the system you look at – generally, this is part of the synchronization and equalization efforts every system has to make. And it really differs! For example, some continuously transmitting system will not insert any specific symbols for that, but depend on the receiver doing statistics. Others have a "silence" period followed by known symbols; other systems are packetized and depend on preambles with specific autocorrelation properties (but not necessarily 100% defined content). Other systems depend on a pilot tone, and there's many more options.

Digital communications are defined by standard. Here some of them:

3GPP for mobile phones
ETSI with DVB for satellite broadcasting and terrestrial TV
IEEE 802.11 for Wifi

Inside each of those standards, there is a part which describes modulation schemes and constellation diagrams.
All wireless vendors implement them as it is written so they are compatible between themselves.


Determine whether a signal is periodic or not and get fundamental period Here is my signal

Cos(n/2)*cos(pi*n/4)

cos(n/2) has period 4pi and cos(pi*n/4) has period 8

Now, the question is will the signal be periodic for fundamental period 32pi ?

$\cos\left(\dfrac{n}{2}\right)*\cos\left(\pi*\dfrac{n}{4}\right)$

Now finding out period of $\cos\left(\dfrac{n}{2}\right)$
$\cos\left(\dfrac{n}{2}\right)=\cos\left(2*\pi*\dfrac{n}{N1}\right) \implies \dfrac{1}{2}=2*\dfrac{\pi}{N1} \implies N1=4*\pi$

Now finding out period of second function $\cos\left(\pi*\dfrac{n}{4}\right)$
$\cos\left(\pi*\dfrac{n}{4}\right) = \cos\left(2*\pi*\dfrac{n}{N2}\right)\implies \dfrac{\pi}{4}=2*\dfrac{\pi}{N2} \implies N2=8$.

Now for total function period is LCM of $(N1,N2)$ i.e, $\mathrm{LCM}(4*\pi,8) \implies32*\pi$.
Therefore, the period of the function is $\boxed{32*\pi}$.

You can try using this:

(1) cos(a+b) = cos(a)*cos(b)+sin(a)*sin(b)
(2) cos(a-b) = cos(a)*cos(b)-sin(a)*sin(b)

Then you can do (1)+(2) and get that

cos(a)*cos(b)=1/2(cos(a+b)+cos(a-b))

Then you can identify the period of cos(a+b)-cos(a-b) by using this for example
http://fourier.eng.hmc.edu/e101/lectures/Fundamental_Frequency.pdf

Hope this helps


concept of image rotation i am trying to understand how does the image rotation occur? i studied and came to know that we just multiply the image coordinates with rotation matrix which gives us the value of new coordinates. ok that is fine i get this. But main thing which is confusing me is the pixel values. i read that the pixel coordinates which we get after rotation can be non decimal hence interpolation is needed to get value at that point for this we rotate image back to original coordinate and find value by interpolation. But if we will rotate it back to original coordinate the the pixel will regain its original coordinate and at that coordinates value w are already having the value then why to do interpolation.please somebody clear my confusion.

When you multiply the image matrix with the rotation matrix, you are essentially multiplying an integer matrix with a real number matrix. The result is a real matrix, which has to be quantized back to a discrete matrix.

What is being quantized? The coordinates themselves. The integer coordinates of the previous position of the pixel is now a real valued coordinate, which must be converted back to an integer in order to be placed on the screen.

If you were asked to rotate it back to the original angle you could have been right in questioning the futility of interpolation. However, the user might not want to rotate it back, so you have to perform interpolation or quantization to be able to use the rotated image. If you are concerned about image fidelity, you might want to save a copy of the original real valued matrix in memory or file until or opertions are finished and the user exits program.

Here is a useful link you might want to peruse.

Image is a matrix.You can see the Jacobian matrix  in linear algebra to understand how the shifting of the Axes occurs


Expectation of deterministic signals Maybe this a stupid question, but I came along with this point and I would like to confirm that I am not wrong.

I found some similar posts in StackExchange, like:

Expectation and Auto-correlation of an Output from a Matched Filter.

In one of the answer, it is stated "as S(t) is determistic[...]" and then the deterministic signal is taken out of the expectation.

and we know that for instance if I have a constant value, I can take it out of the Expectation operator, e.g.:

Where lambda is a random variable and "a" is a constant. Ok, so, if "a" is a constant, then fair enough, but what happens when "a" is a deterministic signal, but it is not a constant, for instance, what happens if "a" is a squared subcarrier?

Can I still safely say that the expectation will be the product of the deterministic signal multiplied by the expectation of the random variable?

Thanks, guys!

If $x(t)$ is a deterministic signal, then you have $E[x(t)]=x(t)$. Furthermore, if $Y(t)$ is a random process you have $E[x(t)Y(t)]=x(t)E[Y(t)]$.

A hand-waving example might be helpful, let's assume $x$ is Gaussian with mean $m$ variance $\sigma$

$$
E\{ x\}=m = \frac{1}{\sqrt{\pi}\sigma}\int_{-\infty}^{\infty} x \exp(-\frac{1}{2} \left( \frac{(x-m)}{\sigma} \right)^2 ) dx
$$
$m$ doesn’t depend on $\sigma$ as long as $\sigma>0$. If we left $\sigma \Rightarrow 0$
$$
\frac{1}{\sqrt{\pi}\sigma} \exp(-\frac{1}{2} \left( \frac{(x-m)}{\sigma} \right)^2 )\; \Rightarrow \; \delta(x-m)
$$
where $\delta(t)$ is the Dirac Delta, You can think of a deterministic variable as a special case of an infinite SNR random variable or perhaps that a deterministic variable has a pdf like:
$$
p_x=\delta (x-m(t))
$$
So going back to your factoring out your constant, You should be ble to show without too much trouble, If $x$ and $y$ are independent,
$$
E\{xy\}= E\{x\}E\{y\}
$$
then
$$
E\{a \lambda\} =  E\{a\}E\{\lambda\}= aE\{\lambda\}
$$
The trickiest part of Expectations is knowing when an expectation is a deterministic variable or a random variable. The conditional expectation,
$$
E\{x | Y\} \; \text{where}\; Y \; \text{is a random variable, the expectation is a random variable}
$$
$$
E_{y}\{E_{x} \{x | Y\}\} \;\text{ is a deterministic quantity.}
$$
Unfortunately most books don't carry the subscripts of which quantity the expectation is operating on.


Simulation of a pressure field developed by a song Greetings!

What I am doing:

In short, I want to write a program that is fed a song as a wav file, and spits out an animated contour plot indicating the relative pressure/loudness/quantity indicating amplitude at various points in a plane. Each speaker is modeled as a point source a set distance from each other hanging in an infinite void. My weapon of choice is Matlab.

How I am doing it:

Poorly. Song is read and split into channels, then cut into little slices. Each slice is FFT'd and the amplitudes and phases are split. For each point in the simulation space every frequency is taken, wavelength calculated and used in combination with the phase to determine the contribution of that frequency to the dimensionless loudness/pressure/sound intensity indicating quantity in that cell (the exact unit isn't that important - explained in next section). Do this for each slice of song, draw contour plots for each slice, and animate!

What my problem is:

The above seems kinda clumsy, and my preliminary results tell me I'm likely neglecting something huge or making some kind of fundamental error. Please note the idea here isn't to make a simulation for speaker setup or any kind of real physical use, I want to draw pretty pictures of what songs do with some basis in reality (specific application: psytrance). This is why the units are not that important, whether I use the normalised voltages in the wav file or convert it to dB or pressure... apart from looking a bit different in lin or log space the relative difference should be the same, at least for my pretty-picture-drawing-purposes. Once the program works as planned some number-fudging will occur to make it look better.

My Questions

1) Can someone suggest a better way of doing this? I feel like I'm doing this in a too complex manner, but I don't even know where to begin looking for a different method to accomplish my goals.

2) Am I right in doing phase corrections with the values given in the FFT? Gut says the error is here - Fourier coefficients are returned as complex numbers a+bi - for the use I have described here, how would I use this to yield "amplitude" and phase angle? It's a simple question... but humour me, self-kicking might follow your answer.

Edit: On second thought this may not be the most appropriate place to post this (at least for part 1 of the question), if someone can suggest a better place, please do.

For a perfect isotropic point source the pressure or air displacement at a point is just the pressure/displacement stimulating signal delayed by the time it took the disturbance to travel to that point.

That means you get spherical shells of constant displacement, and the displacement is the signal at the time $t-d/c$ where $d$ is the distance from the source and $c$ is the propagation speed.

Addendum:

If you do not require the full resolution of the pressure field that would allow you to actually reconstruct the sound from the field then you may want to average the pressure over spatial volumes of constant size. That is not equivalent to a time average of the power, because the influence of spatial coarse graining on the pressure history depends on where you are. Only in the far field where the radius of the constant pressure sphere is much larger than the typical coarse graining scale, time averaging ans spatial averaging become equivalent.

So it depends on what you want exactly and what you need it for. The physics involved is very simple for a perfect point source, but gets very complicated as soon as you have participating objects. Maybe you can give a more detailed description of what you are trying to achieve with this.

I think your method is kind of over kill.  All you really need is the instantaneous power for each point (or each block) in your music.

You can use the hilbert function to get the instantaneous amplitude of your signal, which will do as a representation of the instaneous power.
Split your signal into blocks (cells, as you called them) then do abs(hilbert(x)) where x is a cell.  Average if the length of a cell is all the time resolution you need, or use each value in the sequence as the amplitude at a single sample.

Now, for your display.

A flat surface, with a center point that show the amplitude for the current sample (or cell,) with rings around it representing the amplitudes for the last samples (cells) - each ring is one sample (cell) into the past.

Alternatively, use the distance from your point source to the surface to determine how far into the past each point displayed should be.  Straight ahead (angle 0 degrees) is shortest, and represents now.  At ten degrees off center, you would have 1/(cos(10)) = distance of 1.015 times the straight ahead distance.  Use the speed of sound and your simulated distance from speaker to the surface to figure out how much time (in samples or cells) that relates to.  That will produce a ring on your surface representing the amplitude for a certain time in the past.

For 1 meter, the ten degrees would be 3 milliseconds or 135 samples at 44100 samples per second (typical CD sampling rate.)  70 Degrees at one meter would be about 8.8 milliseconds or 389 samples.

Simulating the distance from the speaker to the surface will tend to produce relatively smooth gradients in the center that get sharper the further out you go.
Simply plotting rings based on cells would give you an image with more jumps in the gradient.
Either method gives you room to play around and make things more interesting.


Errors when converting from digital filter object to transfer function I am trying to use the designfilt tools in MATLAB to create a bandpass filter, but to be compatible with other code I need the resultant filter in transfer function coefficient form with B and A. To start with I took an example from the filter design gallery, and just changed the sampling rate to match that of my real data.

Fstop1 = 150;
Fpass1 = 200;
Fpass2 = 300;
Fstop2 = 350;
Astop1 = 65;
Apass  = 0.5;
Astop2 = 65;
Fs=16384;
% Fs=1000;

d = designfilt('bandpassiir', ...
  'StopbandFrequency1',Fstop1,'PassbandFrequency1', Fpass1, ...
  'PassbandFrequency2',Fpass2,'StopbandFrequency2', Fstop2, ...
  'StopbandAttenuation1',Astop1,'PassbandRipple', Apass, ...
  'StopbandAttenuation2',Astop2, ...
  'DesignMethod','butter','SampleRate', Fs);

This filter seems to be behaving itself, isstable returns 1 and the $z$-plane seems acceptable:

However, when I try to get the coefficients using [B,A]=d.tf; the resultant filter is no longer stable and the $z$-plane looks totally different:

The values for B are incredibly small around e-26, so I am imaging that it might be to do with gains of the coefficients? If I set Fs back to the original $1\textrm{kHz}$, this two $z$-planes plots are the same.

Can someone help me with sorting this problem? Is the fact that the transfer function form is unstable speak to the stability of the filter when it is in the digitalfilter object form?

edit

The documentation for the butter function even describes these errors here

Generally you should never use the Transfer Function syntax. Look how your poles are closely placed next to the unit circle. With high order Butterworth filters you will see the effect of round-off numerical errors. What is the high order? Well 5 is already high and your filter has 30.

Instead you should use SOS or Zero-Pole-Gain representation. For example try to visualise the same filter but with SOS matrix:

fvtool(d.Coefficients) % or: zplane(sos2zp(d.Coefficients))

You will see that there are no visible numerical errors.


Output of KISS FFT interpretation I'm receiving PCM data trough socket connection in packets containing 320 samples. Sample rate of sound is 8000 samples per second. I am doing with it something like this:

int size = 160 * 2;//160;
int isinverse = 1;
kiss_fft_scalar zero;
memset(&zero,0,sizeof(zero));
kiss_fft_cpx fft_in[size];
kiss_fft_cpx fft_out[size];
kiss_fft_cpx fft_reconstructed[size];

kiss_fftr_cfg fft = kiss_fftr_alloc(size*2 ,0 ,0,0);
kiss_fftr_cfg ifft = kiss_fftr_alloc(size*2,isinverse,0,0);

for (int i = 0; i < size; i++) {
    fft_in[i].r = zero;
    fft_in[i].i = zero;
    fft_out[i].r = zero;
    fft_out[i].i = zero;
    fft_reconstructed[i].r = zero;
    fft_reconstructed[i].i = zero;
}

// got my data through socket connection

for (int i = 0; i < size; i++) {
     // samples are type of short
     fft_in[i].r = samples[i];
     fft_in[i].i = zero;
     fft_out[i].r = zero;
     fft_out[i].i = zero;
 }

 kiss_fftr(fft, (kiss_fft_scalar*) fft_in, fft_out);
 kiss_fftri(ifft, fft_out, (kiss_fft_scalar*)fft_reconstructed);

 // lets normalize samples
 for (int i = 0; i < size; i++) {
     samples[i] = rint(fft_reconstructed[i].r/(size*2));
 }

After that I fill OpenAL buffers and play them. Everything works just fine but I would like to do some filtering of audio between kiss_fftr and kiss_fftri. Starting point as I think for this is to convert sound from time domain to frequency domain, but I don't really understand what kind of data I'm receiving from kiss_fftr function. What information is stored in each of those complex number, what its real and imaginary part can tell me about frequency. And I don't know which frequencies are covered (what frequency span) in fft_out - which indexes corresponds to which frequencies.

I am total newbie in signal processing and Fourier transform topics.

Any help?

What you might want to investigate is FFT fast convolution using overlap-add or overlap-save methods. You will need to expand the length of each FFT by the length of the impulse of your desired filter. This is because (1) FFT/IFFT convolution is circular, and (2) each index in the FFT array result corresponds to almost all frequencies (a Sinc shaped response), not just one (even if mostly concentrated near one), so any single bin modification will leak throughout the entire frequency response (except certain exact periodic frequencies).

One problem might be that the  kiss_fftr() function needs a regular kiss_fft_scalar array and not a kiss_fft_cpx array as the input argument.


Why is random noise assumed to be normally distributed? From residual in the linear regression to noise in signal processing are assumed to be normally distributed? By considering them as normally distributed we are kind of telling the pattern in the noise but shouldn't noise be considered random. This seems contradictory to me as on one side it is random then on the other side their distribution is considered normally distributed. Shouldn't the noise distribution be just random?

I believe there is some lacking in my understanding of the concept of statistical distribution which has lead me to this confusion, or I am looking at it all wrong.

One more example- when one augment data by adding Gaussian noise then it is not expected to change the overall distribution of data, why?

Starting at an even more basic level than the other (much smarter) answers, I'd like to pick up on this part of the question:

  This seems contradictory to me as on one side it is random then on the other side their distribution is considered normally distributed.

Perhaps the issue here is what ‘random’ means?

To be clear: ‘random’ and ‘normally-distributed’ do not contradict each other.  ‘Random’ simply means that we can't predict exactly what the next value will be.  But that doesn't mean we can't make probabilistic statements about it.

Consider two experiments:

If you throw a (fair) die, then it could show any number from 1 to 6.  We can't tell which number will come up, but we can say that all numbers are equally likely (i.e. the distribution is uniform).
If you throw two dice and take their sum, that can be any number from 2 to 12.  Again, the sum is still random — we can't predict what it will be — but we can say that those values are not equally likely.  (For example, 7 is six times more likely than 12.)  So in this case it has a non-uniform distribution.  (You can plot all the probabilities; they take on a peaked shape a bit like a normal distribution.)

So there's no contradiction: both cases are random and have a known distribution.

In fact, most things that are random tend to have a non-uniform distribution: electrical noise, weather, the wait for the next bus, voting patterns…  Being able to make general statements about them without being able to predict the exact values is one of the strengths of statistics.

(As for why you often end up with a normal distribution, that's a result of the Central Limit Theorem, which says that when you combine many independent random variables, the result tends towards a Gaussian (normal) distribution.  So you see that crop up a lot.)

Noise is not random. It is fractal in nature.

Mandelbrot discovered that while working at IBM. And knowing that led to the improvement of dial-up modems, among other things. Before that, 9600 baud was out of reach.


How does spectral density of a power signal change when it is multiplied by $\cos(2\pi f_ct)$? Suppose we have a power signal $n(t)$ whose power- auto-correlation function is: $R_{n}(\tau)$.

Now we multiply the signal by a cosine as
$$\cos(2\pi f_ct) \cdot n(t)$$
whose auto-correlation is $R_{n_i}(\tau)$.

How is $R_{n}(\tau)$ mathematically related to $R_{n_i}(\tau)$ ?

$$R_{n}(\tau) = E(n(t)n(t-\tau))$$
$$R_{n_i}(\tau) = \cos(2\pi f_ct)\cos(2\pi f_c(t-\tau))E(n(t)n(t-\tau))=\cos(2\pi f_ct)\cos(2\pi f_c(t-\tau))R_{n}(\tau) $$
Using $\cos(a-b) = \cos a \cos b + \sin a \sin b$, we get
$$R_{n_i}(\tau) = \cos(2\pi f_ct)[ \cos(2\pi f_ct ) \cos(2\pi f_c\tau)+ \sin(2\pi f_ct) \sin(2\pi f_c\tau)]R_{n}(\tau) $$
or
$$R_n(\tau) = \frac{R_{n_i}(\tau)}{\cos(2\pi f_ct)[ \cos(2\pi f_ct ) \cos(2\pi f_c\tau)+ \sin(2\pi f_ct) \sin(2\pi f_c\tau)]}$$


Limitations of Arithmetic coding? I have been working on an assignment for university in which I have to compress an 80K message consisting of only two symbols, 'a' and 'b', having probability values of 0.6 and 0.4.

I have successfully written my code in Matlab and it is working efficiently for 10-20 symbols. However when I pass my 80k symbol message through my code it stops updating at the 1057th value due to the Range parameter, which achieves a value of 10^-324 (the theoretical limit of type double float).

My question is, is this the limitation of arithmetic coding when the range becomes so small that we cannot work with it or am doing something wrong?

I am working on arithmetic encoder, so no answers related to range decoding will work.

You should not be using floats or doubles to store your message. You need to write your compressed message in binary.


Discrete Fourier Transform by hand I have an assignment where I'm given the DFT of a sequence $x[n]$ as $X[k]=\{4,3,2,1,0,1,2,3\}$ and also
$$y[n] = \left\{ \begin{array}[cc]
xx[n/2] & \text{if n is even} \\
0 & \text{otherwise}
\end{array} \right. $$

and I'm supposed to find and sketch the DFT of $y[n]$.

So $y[n] = \{x[0], 0, x[1], 0 ... x[7], 0\}$ and it's not complicated to find $Y[k]$ if we know $x[n]$

I know how to use the definition of the DFT and IDFT to calculate $x[n]$ but it's a tedious task to do by hand, especially when the sequence is longer than a few items. Is there a quicker way to calculate the DFT and IDFT by hand without using a program like Matlab?


Is there any alternative basis for a fourier-like transform? Fourier transform of a continuous signal is just the projection of the signal on the sinusoidal family for imaginary part and the same family with phase offseted by a quarter of period for the real part.
Beside the magical Euler's formula $e^{ix} = \cos x + i\sin x$ that ease the analysis of the transform and the fact that sinusoidal functions happen to be the solution of a prefect harmonic oscillator, I always wondered where the Fourier supremacy came from. Sinusoidal signals aren't so common in everyday's life. Even with human speech, vocal tract function isn't sinusoisdal.
For example, one transform basis that would seams interesting to me would that with a saw tooth function as signals with fast attack and slow release aren't uncommon...
I guess that the transform need to be inversible and that might not be the case for every basis, but aren't there any alternative basis that provide a bijective forward and inverse transform ? (And if so, what would be the conditions ?)

I will elaborate a little. I'm looking for a basis family of periodicals functions named $B(t)$ for example with B of unitary periods such as : $$ S(f) = T\{s(t)\}(f) = \int s(t) B\left( t * f\right) dt + i \int s(t) B\left((t - \frac{1}{4})*f\right)dt$$ with T inversible in some way...

The most important reason for the "Fourier supremacy", as you called it, has been pointed out by Jazzmaniac in his comment. The complex exponential $e^{st}$ with $s=\sigma+j\omega$ is an eigenfunction of the convolution operator:

$$e^{st}*f(t)=\int_{-\infty}^{\infty}e^{s(t-\tau)}f(\tau)d\tau=
e^{st}\int_{-\infty}^{\infty}e^{-s\tau}f(\tau)d\tau=e^{st}F(s)$$

where the eigenvalue $F(s)$, if it exists, is the Laplace transform of $f(t)$, and for $\sigma=0$, i.e. for $s=j\omega$, $F(j\omega)$ is the Fourier transform of $f(t)$ (again, if it exists). Even though there are several important differences between the Fourier and the Laplace transform, they are deeply related, and for answering your question their differences are not important. By analogy, the same is true for the discrete-time Fourier transform and for the $\mathcal{Z}$-transform in the discrete domain. So when in the remainder of my answer I mention the Fourier transform, you might as well replace it by Laplace transform or $\mathcal{Z}$-transform, and everything will remain valid.

So why should we be concerned with eigenfunctions of the convolution operator? Because convolution describes the input-output relation of linear time-invariant (LTI) systems, and these systems are often very good approximations of practically important systems such as frequency selective filters, integrators, differentiators, etc.

Apart from being an invaluable tool for the analysis of LTI systems, the Fourier transform also has an important physical interpretation as the spectrum, i.e. the frequency content, of a signal. Furthermore, the FFT is a highly efficient algorithm for the computation of its discrete version for finite length signals.

So why don't we just say we're doing fine with the Fourier transform and its siblings? Because there are applications where we are looking for other properties of a transform. One example would be time localization in signal analysis, since the basis functions of the Fourier transform have no time localization at all because of their infinite extent. In such cases other transforms have become widely used, such as the short-time Fourier transform and the wavelet transform. For discrete-time signals, filter banks provide very general signal expansions. They can be designed to match the desired properties of the signal expansion.

Since you mentioned the saw-tooth function as a possible basis function I would like to point out to you the Slant transform (mentioned in this answer), which is used in image coding. Furthermore, there is the related Slantlet transform, which is a special type of discrete wavelet transform with piecewise linear basis functions. You can see the basis function in Fig. 5 of this paper.

There's plenty. Laplace transform is just one example. There's uncountably many bijective transformations. So the question would be if there are other useful transforms. Turns out, there are some.


What type of circuit is responsible for convolution in the classic analog telephone? I'm interested in learning how telephones work, so I did a little bit of reading about signal processing. When I came up with the word convolution, I quickly realized the importance of this term.

To begin with, I want to know how classic analog telephones worked. The apparent simplicity of their design appeals to me.

What kind of circuit was responsible for convoluting the microphone impulses in a classic telephone receiver, so that this signal could be transmitted over wires?.

Convolution is a mathematical abstraction describing how a linear, time-invariant system affects a signal going through it.

Sometimes one explicitly designs a system to convolve a signal by a predefined impulse response (for example when building a digital filter); but more often than not, convolution is used to model various physical processes involved in a system. These physical processes can be transmission delays or dispersions, the limited bandwidth of an amplifier, transmission medium or electro-mechanical system, the explicit use of passive R, L, C network to achieve some filtering, etc... Thus, it might not be out of place to find convolution used in the description of some elements of a telephone system, but it's a modeling tool, not an actual process.

Here is an analogy: the trajectory of a cannonball is a parabola, but it makes little sense to ask which mechanical device in a cannon is responsible for computing $y = ax^2 + b$.

In a conventional 20th Century telephone, the "transmitter" (microphone) consisted of a small vial containing graphite powder, with two electrical contacts poking through the otherwise insulated vial wall.  Vibrations from sound would push the graphite particles closer together or pull them farther apart, creating a variable resistance (as "seen" between the two contacts).  Applying a voltage to the receiver would thus result in a variable current that was the "signal".


How can I automatically identify and remove a fundamental frequency and all of its harmonics? I have some microscopy data that is contaminated by a heartbeat artifact that I'd like to remove. The data consists of a large time series of images captured at ~60Hz.
Here's a small example clip in GIF form:

I have taken the average pixel intensity over time, and computed the periodogram using Welch's method:

As you can see there is a sharp peak at ~1.8Hz which is likely to correspond to the heart rate (~108 beats/min). There are also a bunch of harmonic peaks at integer multiples of 1.8Hz. The exact heart rate is likely to vary from dataset to dataset, but I can specify a biologically plausible range as shown by the shaded area on the periodogram.
What I'd like to be able to do is:

Automatically detect the fundamental frequency corresponding to the heartbeat, and all of its harmonics
Filter the data so as to remove the fundamental and all harmonics.

At the moment I can solve point 1 very crudely by finding the largest peak in the periodogram, then multiplying it by $1, 2, ..., N$ where $N$ is the estimated number of harmonic peaks, but I'm sure that there must be a better method than this hack.
Regarding point 2, I came across this question which mentions using a comb filter to remove a fundamental and all of its harmonics. Is this the best method to use? One important consideration is that I will have to apply the filter to each pixel timeseries in a large array, so a computationally efficient method would be highly desirable.
Example data

I've uploaded a 1D timeseries of example mean
pixel values in .csv format here.
I've also added a sample of 1000 frames (spatially decimated 2:1) in a compressed .mat file here

Your method isn't that bad for a first attempt.

However, the following method tends to work better:

Search for local maxima
Check for maxima that are close (2 or 3 bin spacing) and merge them
Create a few hypotheses about the fundamental frequency. You currently assume the highest peak is the fundamental frequency, which is one hypothesis. You should also check for the possibility that the highest peak is the first harmonic, i.e. that there's a smaller peak at half the frequency of the main peak. You may also have other cases to consider, using your knowledge of the problem at hand (Interference? Irregular heartbeat?).
Assuming each of those hypothesis, find the ground frequency by fitting parabola to each harmonic peak. Each peak will produce a slightly different estimate due to noise, but these errors are uncorrelated and average out. One of the hypotheses will lead to much better fits, pick the ground frequency predicted by this hypothesis.
Using the ground frequency you found in step 4 as a given, refit the parabola around each peak to estimate the height of the peak. Note that the peaks probably fall between two bins.
You now have the location and strength of the fundamental and its harmonics, but not the phase. It's probably easiest to find the phase of the fundamental, subtract that, find the phase of the first harmonic, etc.

The core reason that this works better is step 4. Any wrong hypothesis for a ground frequency will fail horribly as you're trying to fit peaks around the location of predicted harmonics. Say you've got a peak at 2 Hz. This could be the ground frequency or a first harmonic. When you test the "first harmonic hypothesis", i.e. whether the ground frequency is actually 1 Hz, you fit parabola's to the data around 1,2,3,4,5 ... Hz. If that hypothesis is wrong, you get garbage around 1,3,5 Hz. If it's right, you might find peaks near 1.1 Hz, 2.2, 3.3, 4.4 and 5.6 - which suggests that the actual ground frequency is 1.12 Hz.

1. Automatically detect the fundamental frequency corresponding to the heartbeat, and all of its harmonics

You can take a local average across the DFT to find the heartbeats. If a certain point in this group is greater than threshold and the greatest value in some range around it then it is a heartbeat or harmonic.

2. Filter the data so as to remove the fundamental and all harmonics.

You could just use the previously located indices of heartbeats and stitch across them linearly. If that's not too crude for you.

What is the purpose of this part of the operation? Is this a medical analysis that requires the absence of the heartbeat to look for other (perhaps) hidden data? Or is this par just for aesthetics?

I imagine a comb filter will be somewhat uncontrollable. It will require careful adjustment of the feedback parameter. It could be viable.

EDIT: You require a pitch detection algorithm? I wrote one of these a few years ago by ignoring the (close to) DC terms and finding the first local maximum. Then using the bins on either side I could use a quadratic interpolation and find the location of the maximum frequency more accurately than the resolution allowed by the frequencies of the bins themselves.

Would it be more accurate to utilise harmonics? It would probably just add an opportunity for errors. Though you can try getting a value for the fundamental with this method and finding the maximum at double, triple etc and then using a similar method as before:

  find local maximum 2. quadratic interpolate to find inter-bin value for frequency.

If you veto the quadratic step and just take the local max you will not get an accurate frequency and moving up the harmonics would help.


What is meant by 20 dB signal-to-noise ratio? I'm reading a paper in which a discrete signal

$$x(n) = s(n) + w(n)$$

is considered. $s(n)$ is a known, deterministic series, and $w(n)$ is white noise with zero mean. The authors write that

  signals were generated with an SNR of 20 dB

What does this mean?

What does signal energy mean? There seems to be several ways of defining this, but no attempt to do so is done in this paper.
What does 20 dB SNR mean? I suppose that either $20=10\log{E_s/E_w}$, or $20=20\log{E_s/E_w}$, but am not familiar enough with these concepts to be sure.

A discrete signal is often interpreted as an analog voltage signal. The signal x(n) should be considered as a sampled version of x(t) with a unknown sampling frequency.

Q1:

So when we're talking about the energy of the signal, it is presumed that the voltage signal is obtained in a circuit with an 1Ω resistor in series, so that v(t)^2/R = P.

Q2:

Signal-to-noise ratio is used in many areas and is defined as the logarithmic power ratio of the signal and noise, i.e. SNR = 10*log(Ps/Pw). The other SNR = 20*log(Ps/Pw) is not true, because the Ps and Pw are power values. The latter is used when the signals are represented as "voltage" signals, i.e. SNR = 20*log(s(n)/w(n)) = 10*log(s(n)^2/w(n)^2).

Even though the 20*log(ratio of arbitrary quantities) form is almost ubiquitously used when defining SNR in the DSP realm, I maintain it is incorrect, and the 10*log(ratio) should always be used. The definition of dB is 10*log(ratio of powers). In the analog realm, power is V^2/R, where V and R are voltage and impedence across that voltage respectively, so we have dB = 10*log((V1^2/R1) / (V2^2/R2)) = 20*log(V1/V2) + 10*log(R2/R1). Especially at high frequencies (radar, communications, etc), it is generally desirable to match impedences for maximum transfer of power, meaning that R2 = R1. If that is the case then the second term vanishes and we're left with this definition of dB: 20*log(V1/V2). Perfectly acceptable when dealing with analog circuits. But when working with software in the DSP realm, what is an impedence? What is a voltage? Those concepts have no meaning, they do not exist in software, they are HARDWARE concepts. Once we have converted an analog signal to digital, the concepts of voltage and impedence are meaningless. So we can hardly say that impedences match (what would that even mean? the impedence of the circuit sourcing the ADC matches the input impedence of the ADC? LUDICROUS!!), meaning we have no basis for saying that the second term vanishes, and thus no basis for defining a dB as 20*log(ratio of arbitrary quantities). Therefore, it makes most sense to revert to the definition of a dB and use the 10*log(ratio of arbitrary quantities) form.


Find frequency/phase response of a microphone I am trying to find the frequency/phase response of two cheap microphones bought off Ebay. Basically, I want to graph the frequency responses and phase responses of both microphones on the same axis so I can compare their sensitivity in MATLAB.

One approach I have read is downloading and playing a file that consists of pink noise through a speaker and simultaneously recording this using both microphones, then using the two recorded signals and the original pink noise file to get the bode plots. This is achieved by first taking the FFT of the recorded signals and of the original pink noise signal, and then using the bode function in MATLAB. The frequency/phase response of the speaker is unknown, however I'm more interested in the differences between the two microphones, so my thoughts are this wont affect my analysis. Any thoughts on this?

Is there an easier approach to finding the frequency/phase response of two microphones that I am not aware of using MATLAB?

I would say skip MATLAB and learn more about the published frequency responses by looking up your mics here: http://recordinghacks.com/microphones. There are hundreds of mic models with responses as such:

You want to find the transfer function of a microphone. This is a hard thing without some prerequisities. You'll need an anechoic or hemi-anechoic chamber. A flat spectrum excitation source. Don't know if your source is relatively flat? You are out of luck. You need to average also for a smoother response. The command you're looking for in MATLAB is freqz(). Go with the published response when you don't have these things at hand.


What's the ideal FFT window for measuring a group of signals of differing amplitudes but close in frequency? In my system, I need to measure the amplitudes and phases of four lasers, each of a different, known wavelength, chopped at a 50% duty cycle via TTL at frequencies spaced evenly (or not) around some fundamental frequency f0. The lasers are attenuated by some analyte medium, and the resulting laser power (which is measured by a single detector) tells me what I need to know about the medium.
The signals are recovered by performing I/Q demodulation. The application is similar to the one in this question. The signals' individual frequencies can be controlled by me, but should not be more than 6-7 Hz away from f0.
In single-signal applications of my setup, I've used the Exact Blackmann window, which is known to be good for single-frequency measurements, at least according to National Instruments. My question is, is there a more robust choice for recovering four extremely close signals?
Other details: sample rate is 8192 or 16384 Hz for two seconds, f0 is in the neighborhood of 1230 Hz. Typically I have the frequencies evenly spaced, at f0±6 and f0±2, so they're all 4 Hz apart. The input signal is demodulated and windowed four times, once at each of the frequencies of the lasers. Followup question: could spacing them unevenly help?

As for ideal windows for FFT frequency resolution and dynamic range, many assume that the Gaussian window would be best since it has the minimum time-bandwidth product, and due to this it is the filter most often used in spectrum analyzers. However, this is only true when the time domain extends to infinity. For finite time durations the window that actually has the minimum time-bandwidth product is the Discrete Prolate Spheroidal Sequence (DPSS, or Slepian) window which maximizes the energy concentration of the main lobe. The Kaiser window very closely approaches this optimum and both have an controlling parameter to trade frequency resolution (ability to discern two closely spaced frequencies) and dynamic range (ability to discern two frequency tones of different amplitudes). For further details on the comparison of the two, see this link by Julius Smith: https://ccrma.stanford.edu/~jos/sasp/Kaiser_DPSS_Windows_Compared.html
That said, the signal after IQ demodulation could be windowed and then FFT'd to provide the phase and amplitude (relative to the demodulation signal used) of closely spaced modulation frequencies. For this application I would use a Kaiser window with Beta = 8 which provides 80 dB dynamic range of signals captured over a 2 second duration while being able to resolve two signals that are 4 Hz apart.

FFTs seems not the best way of doing it since you need very long time sequences to get enough frequency resolution and you have to deal with significant spectral leakage. Windowing will reduce spectral leakage put also reduce frequency resolution so it's a game of whack-a-mole.
If you only have a few known frequencies a heterodyne detector should work very well:

Build a local oscillator at the desired frequency
Multiply your input with the local oscillator  (complex or cos & sin)
Lowpass filter the output

That will give a complex signal that represents the spectral amplitude over time. The bandwidth of the low pass filter determines the bandwidth of your spectral analysis.
If you don't know the frequency exactly or it can drift, you can wrap a PLL (Phase Locked Loop) around the  local oscillator to track your spectral component.


How to calculate resolution of DFT with Hamming/Hann window? The frequency resolution of a DFT with a rectangular window of size $N$ is given by $f_s/N$. However, when using other window functions like a Hamming or Hanning window the resolution gets worse. How can I calculate the frequency resolution of a DFT when using a Hamming/Hanning window of size $N$?

I found some numbers regarding "effective noise bandwidth (ENBW)" but am not sure what exactly it means. The numbers were 1 for a rectangular window, 1.5 for Hanning and 1.37 for Hamming window. Does that mean the frequency resolution using a Hanning window is $1.5 \cdot f_s/N$?

The frequency resolution of a DFT is indeed significantly affected by the choice of window. Two key parameters to consider when selecting a window are the frequency resolution and dynamic range. The frequency resolution is the ability to discern two closely spaced frequencies at similar power levels (therefore a metric of what that spacing is would be the frequency resolution), while dynamic range is the ability to discern two signals at different frequencies and significantly different power levels; specifically the ability to see a low level signal in the presence of a much stronger signal (therefore a metric of what this power difference is would be the dynamic range). Typically a higher dynamic range is achieved by using appropriate windowing at the expense of frequency resolution. The highest frequency resolution that can be achieved is with a rectangular window (and the lowest loss; any other window will also introduce a processing loss).

The Equivalent Noise Bandwidth (ENBW) is a commonly used metric for frequency resolution. The ENBW is the bandwidth (often given in FFT bins) of a brickwall filter that would result in the same noise power as the DFT "filter" for white noise (the DFT can be described as a bank of filters). For the rectangular window, the ENBW is one bin, but the poorest dynamic range (sidelobes with peaks rolling off at a $1/f$ rate in magnitude.  Any other windowing used would have a higher ENBW.

For a given length N window function $w[n]$, the ENBW in bins is computed as follows:

$$ENBW = N \frac{\sum \left(w[n]\right)^2 }{\left(\sum w[n] \right)^2}$$

See this classic paper by fred harris tabulating the ENBW as well as other paramters of many commonly used windows. fred harris windowing

Below is a slide I have further depicting the ENBW for a rectangular window:

Below compares a N = 30 rectangular window (denoted the Dirichlet kernel) to a Blackman window of the same length. The decrease in level shown of -7.83 dB is the coherent gain for this particular window used. The coherent gain is simply the average of the window weights. The result below is is consistent with the coherent gain of 0.42 as listed for the Blackman window tabulated in the fred harris reference given; $20Log_{10}(0.42) = -7.54$ dB.  What is very clear from this picture is the significant increase in both dynamic range and frequency resolution.

The reason for not being an exact match to the tabulated results is the window length. As the window length approaches infinity the physical result approaches the tabulated, as shown in the figure below for coherent gain. For this reason caution must be applied when using the tabulated results for any smaller length windows.

The frequency resolution does not depend on the type of the window. It is the length that matters. The frequency resolution is

$\qquad\text{frequency resolution} = \displaystyle\frac{\text{sampling frequency}}{\text{window length}}.$

For example, say you have an audio signal with $f_{s} = 44100$ Hz. If the window length is $2\text{ s} = 2\cdot 44100 = 88200$ samples then the frequency resolution will be

$\qquad\displaystyle\frac{44100\text{ Hz}}{88200} = 0.5\text{ Hz}.$

Just to clarify, what I am describing here is the resolution of the DFT.