Time Alignment of 2 Sensors Sampling the Same Signal with Different Hardware Delays I have two radios sampling simultaneously, and I am trying to time-align their two signals (they will always start sampling at slightly different times with random delays due to hardware differences). I have a circuit that automatically switches between the antennas and nothing at a set interval, and I'm trying to use the transitions between signal and silence to determine how to time-align the two signals. I tried using cross correlation to find the shift between the two signals, but it's not getting the correct alignment (it's very off). I suspected it was due to the signals being too noisy, so I smoothed them both and tried to make their amplitudes similar, but it still didn't get the right alignment.

Here's a screenshot of an example of the two signals. The image is zoomed in to the region where the circuit switches the radios from nothing to the antennas, and I would like to align the two signals such that these transitions match up.

EDIT:
Here are text files of two more signals that I'm trying to align:
s0.txt   s1.txt. When they're aligned, they look like this:

(note that these are not the same signals in the previous picture). I found the optimal delay to be 1069 for the first signal. I used other methods to find the delay (described in a comment below), but I would like to reproduce this result using cross correlation as I believe it will prove to be more reliable / precise.

If the above is a good representation you should just try to infer when there is energy in the signal to align them.
As it seems they start with nothing (Zero value).
Then all needed is just to find where "Something" happens. This could be easily done with high resolution (Few samples).

Regarding Cross Correlation, try to normalize both signal to have the same Maximum Value (Normalize both to have amplitude of 1).

In the case above you'll be able to say where the alignment point is.
Pay attention that if the whole signal is harmonic then the cross correlation is harmonic as well.


Frequency detection in audio and what parameter or set of parameters help/limit frequency detection I'm trying to learn audio and voice/speech processing techniques out of interest. I intend to learn, experiment with algorithms and design, implementation of such systems, frameworks using programming languages like C/C++, Python, Matlab. Right now, to begin with, I'm only interested in audio or speech detection and not in speech recognition. I did some basic reading and research into this. Would like to understand a few basic concepts:

How to detect the presence of human speech (voice frequency of male being between 85 Hz to 180 Hz, and that of female voice frequency between 165 Hz to 255 Hz but maximum range being 20 Hz to 20 kHz), differentiate them from other frequencies, identify noise etc?
Being aware of the sampling theorem and the limitations of using lower sampling rates, how to know if any error was introduced during digital transmission or during reconstruction of the signal in analog domain from the digital domain (i.e., error detection) ?
How to detect tones (produced digitally or externally recorded analog tones) and music notes? Do they follow similar techniques as speech? If no, how are they different?

I seek guidance, if not direct answers, on how to proceed with the learning the concepts and reading material. Any help in this regard will be grateful.

I'll give a few answers with the qualification that this isn't really my specialty area:

1) If you are using a DFT (FFT), then you should keep your sample frames relatively short.  Speech utterances are of short duration and the DFT is best with steady tones across the entire frame.

2) Your sampling rate will likely be 44.1kHz (CD quality) or 8kHz, either of which is sufficient that the Nyquist frequency is nowhere close to the ranges you are talking about.  In other words, this shouldn't be a concern.

3) Music tones tend to be of longer duration and a bit steadier.  But long drawn out vowels fit the same pattern.

My advice would be get some software, like Audacity, that has a good spectrogram function and look at speech and music samples that way.  Or write your own.  Anyway, that is a good starting point to build up some familiarity with what various sounds and voices look like in the frequency domain.


Why is the Fourier Series a special case of the Fourier Transform and not the other way around? I was reading a text book on the frequency domains in signal processing and my understanding is that the Fourier Transform considers signals that are a-periodic in time while the Fourier Series considers periodic signals in time (such as a sum of sinusoids).
But the book clearly mentions that the FS can be seen as a special case of the FT and that doesn't make sense to me.
My understanding is that the FT is a special case of FS where periodicity is just infinity. But the book says its the other way around without providing its own reasoning.
It makes no sense and I can't find answers online

That's more of a philosophical question. It makes no practical difference if you choose one answer over the other.
Personally, I would agree with the author. I find it  mathematically a lot easier to derive the Fourier Series from the Fourier Transform than the other way around. It also makes more sense to me to say "periodic signals are a subset for all signals" as compared to "all signals are periodic, some just have a period of infinity". That kind of defeats the purpose of having periodicity as a specials case.


How to quantify asymmetry of a signal? I'm looking for a suitable measure in order to quantify the asymmetry of a signal. I've read this: Detect signals degree of imbalance and this: How to measure the agreement between to curves? but I can't see how they apply to my case.

So this is the part of the signal that I am interested in:

The red line shows the original signal and the black line shows the right part of the signal mirrored about the time=0 axis. What I need to do is find a way to quantify the asymmetry of the red line, i.e. a measure of its difference from the black line.

What I have tried doing is taking the square root of the sum of the squares (SRSS) of the differences of each point. That is:

$$
SRSS = \sqrt{\sum\limits_i (y_{black,i} - y_{red,i})^2}
$$

Unfortunately, this measure is affected by the number of points that I take for each curve. So for example if the resolution is half, SRSS will decrease even if I'm comparing the same lines.

Is there any other measures of asymmetry that you can suggest? Ideally, they should be a single number and not be affected greatly by the resolution of the lines.

Assuming you're looking for symmetry around $x=0$, you can decompose any function into a symmetric and an antisymmetric part:

$$ f(x) = \frac{1}{2}\left(f(x)+f(-x)\right) + \frac{1}{2}\left( f(x) - f(-x) \right)$$

Calling $f_+(x):=\frac{1}{2}\left( f(x) + f(-x) \right)$ and $f_-(x):=\frac{1}{2}\left( f(x) - f(-x) \right)$ you can easily verify that $f_+$ is symmetric while $f_-$ is antisymmetric.

Now you can measure symmetry by taking any suitable norm on your functions. For example the symmetry of $f(x)$ can be measured as

$$S[f]:=\frac{||f_+||}{||f_+||+||f_-||}$$

which ranges from $0$ (fully antisymmetric) to $1$ (fully symmetric).

A possible useful norm on functions is defined as

$$||f|| = \sqrt{\int_\mathbb{R} f(x)^2 dx} $$

which can be discretised to

$$||f|| = \sqrt{ \sum_i f(x_i)^2 \frac{x(i+1)-x(i-1)}{2} }$$


How to convert DFT to DCT I've been trying to find a fast 16 point dct2 and dct3 implementation, however I could only find it in literature expressed as mathematical equations, which honestly I couldn't understand.

However I did find a code generator that outputs dft transforms. The main problem I have is the DFT and Inverse DFT don't have the same numbers going in and out.

// Before DFT              // After DFT and IDFT
inputArray[ 0] = 12;       outputArray[ 0] = 184;
inputArray[ 1] = 12;       outputArray[ 1] = 194;
inputArray[ 2] = 12;       outputArray[ 2] = 178;
inputArray[ 3] = 14;       outputArray[ 3] = 198;
inputArray[ 4] =  8;       outputArray[ 4] = 155;
inputArray[ 5] = 10;       outputArray[ 5] = 141;
inputArray[ 6] = 12;       outputArray[ 6] = 164;
inputArray[ 7] = 12;       outputArray[ 7] = 149;
inputArray[ 8] = 12;       outputArray[ 8] = 138;
inputArray[ 9] = 12;       outputArray[ 9] = 121;
inputArray[10] = 12;       outputArray[10] = 107;
inputArray[11] = 12;       outputArray[11] = 90;
inputArray[12] = 12;       outputArray[12] = 74;
inputArray[13] = 12;       outputArray[13] = 55;
inputArray[14] = 12;       outputArray[14] = 37;
inputArray[15] = 12;       outputArray[15] = 19;

I realized the first 5 or so indexes do equal the inputs when divided by 16, however this trend doesn't continue as you go down.

Is this the expected behavior? Or is there something else I need to do the get a proper conversion?

Also I did find an 8 point dct that works well and gives the proper results, is there anyway to expand that into a fast 16 point dct?

Edit:

The reason I want to find a fast 16x16 DCT is because I'm working on a javascript video codec that supports transparency.

So after inter frame prediction is finished I'm left with a lot of residue. Because it's in javascript, and it's processing 24 frames a second I need the fastest way possible to compress and decompress this residue.

It depends on which DCT (there are more than four of them), the idea is to take your original data and append to it a mirror-reflection copy of that data.  Now your data is twice as long, and you perform an FFT on twice the length.  but you have some symmetry to it that allows you to toss half of the data in the result of the DFT.

exactly how, depends on if it's DCT I, DCT II, DCT III, DCT IV, or the "MDCT" or whatever is the DCT flavor of the month.

looks like something useful here.  it may need further exploration.


Is step detection the correct approach to this problem? what if not? I'm looking for some advice on where and what to start reading for learning to solve this.

I've the time series of the position coordinates (x,y) of an animal in an open field (just a cage). I want to detect the time instants when the animal stops or starts "walking" (i.e.: moving from one place to another, not in the same place).
I  thought it couldnt be so difficult, but I'm finding some trouble solving it.

So, I'm looking for the times from which the signal will be flat, or the instant where the signal starts changing from flatness:

I think this thing may be related with step detection but I'm not sure.
I would start with just one coordinate for simplicity. Is seems that step detection could be what im looking for, but:
a. my signal "baseline" will be different every time. the animal will move and stop somewhere else, and so.
b. the signal can be VERY "noisy" since the animal could and will move a lot while staying in the same place (ie: grooming)
c. this changes can be whether slow or fast, and I need both.

Firstly, i ll be glad if i can solve this problem for just one coordinate, although I will have to look for this changes in the two coordinate system.

so, my questions are:
1. is step detection a good aproach to this problem? what else if not?
2. any suggestion in doing this for both (x,y) coordinates?

thanks in advance

EDIT: I get (x,y) coordinates by acquiring an overhead image with a camera and tracking a led attached to the animal's head. Tracking is done by color filtering.

EDIT2: copy of the data:

https://www.dropbox.com/s/oph33szu891rgrl/data_sample.txt

https://www.dropbox.com/s/idulfdr965eeh7i/data_sample.mat

data format is (x,y,t)

EDIT3:
I've been trying smoothing the data but it is not really what I expected. I need to clean not high frequencies but low amplitudes. The movements I look for can be either fast or slow, but with big amplitude.
here, an example of signal with matlab function smooth()
note that I'm looking for the green moments, which I lose when smoothing

First off, I would start by applying a smoothing/averaging filter to your position data. This will get rid of a lot of the noise from your stationary motion.

To answer your first question, step detection is the way to go, but it depends on how you're detecting the steps. The best way in my mind would be to look at the derivative, the speed, and detect periods of non-zero speed. If you've smoothed out your position data, this should result in periods of time where the animal is moving in one direction, rather than back and forth.

This approach makes two dimensional analysis fairly simple. Once you've calculated your x and y speeds, you can combine them to get an absolute speed, independent of direction. You can use the same detection approach to find the periods of motion, which should give you the answer you're looking for.


How to prove this norm inequality? $$
\left( \sum |z_i|^a \right)^b \geq \left( \sum |z_i|^b \right)^a, \quad \forall a\in [1,2], b\in [2,+\infty)
$$

Actually this is equivalent to
$$
\|z\|_a \geq \|z\|_b \quad \forall a\leq b
$$

Although the question could belong to SE.math, mastering inequalities for $\ell_p$ norms (for $p\ge 1$) or quasinorms (for $0<p< 1$), and their norm ratios and powers, is quite important in  modern signal/image processing: starting from Cauchy-Bunyakovsky-Schwarz, to power means, Rogers-Hölder, Young, Minkowski, Dresher-Gini etc.)  There are scalar, additive and convolutive versions, and reverse forms, that are used in data restoration and deconvolution, especially when using signal sparsity.

If $a_n$ and $b_n$ are positive finite sequences of length $N$, $u,v,w>0$, and $1/u+1/v\le 1/w$, the generalized Rogers-Hölder's inequality gives:

\begin{equation}
\left( \sum_{n=1}^{N} (a_n b_n)^w \right)^{1/w} \le \left( \sum_{n=1}^{N} (a_n )^u \right)^{1/u}
\left( \sum_{n=1}^{N} (b_n )^v \right)^{1/v}\,.
\label{eq_rogers-holder}
\end{equation}

You can find a proof in P. S. Bullen, Handbook of Means and Their Inequalities, Kluwer Academic Publishers, 2nd edition, 2003, p. 187, Corollary 8.
22,

If you now set $a_n=|z_n|^p$, $b_n=1$, $w=p$, $v=q$ and  $1/u=1/p-1/q$, you get
\begin{equation}
\ell_{q}{(|z|)} \leq \ell_{p}{(|z|)}  \leq \ell_{0}^{\frac{1}{p}-\frac{1}{q}}\ell_{q}{(|z|)}   \,.
\end{equation}

$\ell_{0}$ here, not a norm, nor a quasinorm, is the count index, ie the number of nonzero terms in $z$, or the sparsity. In other words, there is  an inclusion in $l_p$ spaces: if $p\le q$, $l_p\subset l_q$.
Your answer is incorrect, because you should have $a\le b$ instead.

One of the basic block is alluded to by  @robert bristow-johnson. If $0<c\le 1$ ($a_n$ positive as above):

$$ \left(\sum a_n\right)^c \le  \sum a_n^c \,.$$

Using additivity, prove it first with two terms:

$$(a_1+a_2)^c \le a_1^c+a_2^c\,.$$

By factoring by the largest of $a_1$ or $a_2$, you end up to prove that for $0<t<1$:

$$(1+t)^c\le 1+t^c\,.$$

Just differentiate the function $t\mapsto 1+t^c- (1+t)^c$ as $c\left(t^{c-1}- (1+t)^{c-1}\right)$, and verify the latter is positive.


Splitting a Complex FIR Filter into Real Filters I'm working on a software defined radio project on a PIC32 MCU. I'm using a CMX 973 quadrature frontend to provide I and Q. I'm independently low passing I and Q with carefully matched analog filters to reduce the overall bandwidth to 192 kHz, then digitizing it with an audio codec.

I'm going to be receiving a couple of closely-spaced channels so I'll need to filter them separately (bandpass), mix, and downsample. I'll also need to do some filtering for pulse shaping and PLL. These operations are all on complex-valued signals (I-Q).

The trouble is the PIC32 DSP library only provides an FIR filter implementation for real-valued filters. I can write my own complex filter implementation relatively easily, but it won't come with the benefit of optimization done by Microchip - I'll have to do it myself. I can do it, but I'd like to avoid it if I can.

Mathematically, is there any way to split a complex signal into I and Q and perform filtering using several real-valued filters which would have the same response as a complex-valued filter (particularly for bandpass filtering)? If so, how, and what are the limits to this approach?

I'd say there are two ways to look at this problem.

If you're receiving a conventional $M$-QAM signal with AWGN noise, then the I and Q streams are actually two independent, real, $\sqrt{M}$-PAM signals. You can filter and process them independently. The same is true over the wireless channel, as long as the fading is flat and you do the appropriate detection in the receiver.
If $s(t)=s_I(t)+js_Q(t)$ is a complex signal and $h(t)=h_I(t)+jh_Q(t)$ is a filter's complex impulse response, then

$$s(t)\star h(t)=s_I(t)\star h_I(t)-s_Q(t)\star h_Q(t)+j(s_I(t)\star  h_Q(t)+s_Q(t)\star h_I(t)),$$

where all the filters involved are real. This follows from the convolution's properties.


Phase error correction for Fourier transform basis vectors I was wondering what is the best way to account for phase noise in a set of basis vectors?

Example:
I have a measured signal, say $f'(x)$ that should be related to a desired spectrum $f(\nu)$ through the following relation (which is a Fourier transform for x-symmetric, band-limited signals, as pointed at in comments below):

\begin{equation}
f(\nu) = \int_0^{x_\text{max}} f'(x) \cos(2\pi x \nu)\,dx \tag1
\end{equation}

This can be expressed with vectors and matrices as $f'(x) = M\cdot f(\nu)$.  Here, the rows of $M$ correspond to an orthogonal set of sinusoidally varying basis vectors.

In my particular problem, the system has constant phase error that results in an actual relation of:
\begin{equation}
f(\nu) = \int_0^{x_\text{max}} f'(x) \cos\left(2\pi x \nu + \phi(x)\right)\,dx \tag2
\end{equation}
When I calculate this particular calibration matrix $M'$ (s.t. $f'(x) = M'\cdot f(\nu)$) for this system and perform the Moore-Penrose pseudo-inverse ($M'$ is non-square), I find very significant errors resulting in my spectra calculated via $f(v) = M'^{-1}f'(x)$.  I have also confirmed that spectral reconstruction via the pseudo-inverse without the phase noise works well (as I would expect).

Assuming I know apriori what my phase errors $\phi(x)$ are, is there a better way for me to reconstruct my spectrum $f(v)$ than through just the pseudo-inverse of $M'$?  Any additional insight or ways of looking at this are appreciated.

Edit:
Assume the matrix $M'$ is measured with a known, single frequency input so that the rows of $M'$ are indeed sinusoidally varying basis vectors, and the spectrum is known for $k$ points, and the signal is measured for $m$ $x$-points, so $M'$ is a $k\times m$ matrix ($k>m$). I could simply measure fewer spectral points, so $M'$ is square, though I had thought oversampling might provide some additional information.

This can be expressed with vectors and matrices as $f'(x) = M\cdot f(\nu)$.  Here, the rows of $M$ correspond to an orthogonal set of sinusoidally varying basis vectors.

Woah, there, that's a big step you're doing. So, you have your band-limited signal; it has to be band-limited, otherwise the $x_\text{max}$ limit in $(1)$ would be illegal.

So, since it's limited in Fourier ($\nu$) domain, it has to have infinite support in time domain ($x$). (Imagine you take an infinitely long signal, and window it with a rectangular function: You need to convolve the Fourier transform of that spectrum with a scaled $\text{sinc}$ function, and that has infinite support)

Now, you're right to consider the Fourier transform as a linear operator $M$. What is very hard to do is justify that $M$ is actually a matrix with discrete rows of discrete frequencies. That integral there doesn't give you discrete, separate sinusoids, it gives you a continuum of frequencies!

Basically, the whole of the sampling theory is to prove that you can do that. Short idea: Thanks to Nyquist we know that we can represent any band-limited signal by sampling that signal at a sufficiently high rate, giving us, indeed, a discrete vector. The Time-Discrete Fourier transform of a signal collapses to a sampling rate-periodic spectrum, and we can show that this is also representable by a discrete vector, leading us to an understanding of the Discrete Fourier Transform (DFT), which actually is just a Matrix with complex numbers in it, of size $N\times N$, $N$ being the length of the time observation vector.

But: that discretization into an actual matrix over $\mathbb C$ requires you to first convert to discrete signals. Then, $M$ must be square, and is indeed (easily!) invertible. So, your problem lies in that somehow, you ended up with a non-square $M$, and that can't happen if you adhere to the sampling theorem.


Fast fourier transform for an audio file I need to do FFT on an audio file(to be exact it is a real number array of audio samples). But the output of FFT should be complex numbers. I have a hard time understanding how to convert this.
For example: input is int input[N], and what to do next to get the output.re[] and output.imp[]? Should I treat input as a complex array only all the imaginary components equal to 0?
Thanks!

Your question is so "basic", so "fundamental", that I may be misunderstanding its true meaning. In any case, hopefully you have some software routine available to you that performs FFTs and will accept your 'int input[N]' sequence as an input. I have no idea what software 'package', or software language, you're using but if variable 'x' is your input sequence then there should be some available command (or subroutine) like

Spec = fft(x)

you can use to compute your FFT spectral results. And there should be no need for you be concerned about defining your input sequence as complex-valued. Of course the length of your input sequence (number of samples) must be an integer power of two. Your FFT results will automatically be in the form of a sequence of complex numbers, so there's no need to worry about that.


Can someone explain waveshaping to me? Can someone explain waveshaping to me?

I only know it shapes waveshapes by passing through some functions, but I don't yet understand e.g. what the plots (such as the following in Melda Production MWaveShaper VST-plugin) mean.

Both axes have dB values so I read the axes to be input and output. However, how is the distortion "visualized" in the graph? How does one know how the graph alters audio?

In the audio domain, waveshaping is simply applying a memoryless nonlinear function to an input signal.

$$ y(t) = g\big( x(t) \big) $$

The waveshaping function, $g(x)$, is most often a continuous function that goes through the origin: $g(0)=0$.

Sometimes $g(x)$ is an odd-symmetry function: $g(-x)=-g(x)$, but it doesn't have to be.  Sometimes you want 2nd harmonic distortion and then $g(x)$ does not have odd symmetry.

I am more interested in the case where the waveshaping function, $g(x)$ is a polynomial of finite order $N$:

$$ g(x) = \sum\limits_{n=0}^{N} a_n\,x^n $$

Note that if $g(x)$ has even symmetry ($g(-x)=g(x)$), then all odd terms are zero ($a_n=0$ for $n$ odd).  If $g(x)$ has odd symmetry ($g(-x)=-g(x)$), then all even terms are zero ($a_n=0$ for $n$ even).  I guess, if you want $g(0)=0$, then $a_0=0$.

Now, if the input is sinusoidal:

$$ x(t) = A \cos(2 \pi f_0 t) $$

then the output is periodic, having harmonic frequency components, with an upper limit to the frequency content.  The highest frequency coming out will be $N f_0$.  So waveshaping with polynomial mapping functions guarantees a limit to the highest frequencies generated and that can be useful in deciding an oversampling ratio necessary to guard against aliasing.  As I alluded to in this answer, if the order of the polynomial is $N$, to insure against aliases folding back into your original passband, one must oversample with an upsampling ratio of at least $\frac{N+1}{2}$.  Oversampling by 4x sufficies for a 7th-order polynomial.

Also note that an even-symmetry polynomial will generate only even harmonics and an odd-order polynomial will generate only odd harmoncs.

$$\begin{align}
y(t) &= g\big( x(t) \big)  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\big(x(t)\big)^n \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\big(A \cos(2 \pi f_0 t)\big)^n \\
\\
 &= \sum\limits_{k=0}^{N} b_k \cos(2 \pi k f_0 t) \\
\end{align}$$

where the $b_k$ coefficients depend on the $a_n$ coefficients and the input amplitude of $A$.  This relationship can be worked out using the Euler identity and the binomial theorem:

$$\begin{align}
y(t) &= g\big( x(t) \big)  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\big(A \cos(2 \pi f_0 t)\big)^n \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}(e^{j 2 \pi k f_0 t}+e^{-j 2 \pi k f_0 t})\right)^n \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n\left(e^{j 2 \pi f_0 t}+e^{-j 2 \pi f_0 t}\right)^n \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} (e^{j 2 \pi f_0 t})^m (e^{-j 2 \pi f_0 t})^{n-m}  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi m f_0 t} e^{-j 2 \pi (n-m) f_0 t}  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \tfrac12 \left(\sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi m f_0 t} e^{-j 2 \pi (n-m) f_0 t} \\ + \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi m f_0 t} e^{-j 2 \pi (n-m) f_0 t} \right)  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \tfrac12 \left(\sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi m f_0 t} e^{-j 2 \pi (n-m) f_0 t} \\ + \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi (n-m) f_0 t} e^{-j 2 \pi m f_0 t} \right)  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \tfrac12 \left(\sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi (2m-n) f_0 t} + \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi (n-2m) f_0 t} \right)  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \tfrac12 \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} \left( e^{j 2 \pi (n-2m) f_0 t} + e^{-j 2 \pi (n-2m) f_0 t} \right)  \\
\\
 &= \sum\limits_{n=0}^{N} \sum\limits_{m=0}^{n} a_n\,\left(\tfrac{A}{2}\right)^n  \frac{n!}{m!(n-m)!} \tfrac12 \left( e^{j 2 \pi (n-2m) f_0 t} + e^{-j 2 \pi (n-2m) f_0 t} \right)  \\
\\
 &= \sum\limits_{n=0}^{N} \sum\limits_{m=0}^{n} a_n\,\left(\tfrac{A}{2}\right)^n \frac{n!}{m!(n-m)!} \cos \big(2 \pi (n-2m) f_0 t \big)  \\
\\
...
\\
 &= \sum\limits_{k=0}^{N} b_k \cos(2 \pi k f_0 t) \\
\end{align}$$

Okay, to get an expression for $b_k$, the only way I can figger it out is for three cases.  We let $k=n-2m$ and observe that since $2m$ is always even, only even $n$ terms contribute to $b_k$ when $k$ is even.  Likewise, only odd $n$ terms contribute to $b_k$ when $k$ is odd.

The expression $\lfloor N/2 \rfloor$ is the floor() function applied to $N/2$.  just round $N/2$ down to the nearest integer.

Case 1, $k=0$ (the DC term)
Only the even terms of $n$ will contribute to $b_0$.  And the only term of the inside summation that contributes is that when $m=\tfrac{n}{2}$.  So letting $n=2i$:

$$ b_0 = \sum\limits_{i=0}^{\lfloor N/2 \rfloor} a_{2i}\,\left(\tfrac{A}{2}\right)^{2i} \frac{(2i)!}{\big((i)!\big)^2} $$

In the following two cases, keep in mind that the $\cos()$ function is even symmetry and these two terms are equal:

$$ \cos(2 \pi (-k) f_0 t) = \cos(2 \pi k f_0 t) $$

Case 2, $k$ even, $k>0$ (even harmonics)
Only the even terms of $n$ will contribute to $b_k$.  The only two term of the inside summation that contribute are those where $m=\tfrac{n-k}{2}$ and $m=\tfrac{n+k}{2}$ .  So letting $n=2i$:

$$ b_k = \sum\limits_{i=0}^{\lfloor N/2 \rfloor} a_{2i}\,\left(\tfrac{A}{2}\right)^{2i} \frac{2(2i)!}{(i-\tfrac{k}{2})!(i+\tfrac{k}{2})! } \qquad k>0 \text{ even} $$

Terms with $i<\tfrac{k}{2}$ will result in $\infty$ in the denominator because factorials of negative integers are infinite.  That effectively modifies the bottom limit:

$$ b_k = \sum\limits_{i=k/2}^{\lfloor N/2 \rfloor} a_{2i}\,\left(\tfrac{A}{2}\right)^{2i} \frac{2(2i)!}{(i-\tfrac{k}{2})!(i+\tfrac{k}{2})! } \qquad k>0 \text{ even} $$

Case 3, $k$ odd, $k>0$ (odd harmonics)
Only the odd terms of $n$ will contribute to $b_k$.  The only two term of the inside summation that contribute are those where $m=\tfrac{n-k}{2}$ and $m=\tfrac{n+k}{2}$.  So letting $n=2i+1$:

$$ b_k = \sum\limits_{i=0}^{\lfloor (N-1)/2 \rfloor} a_{2i+1}\,\left(\tfrac{A}{2}\right)^{2i+1} \frac{2(2i+1)!}{(i-\tfrac{k-1}{2})!(i+\tfrac{k+1}{2})!} \qquad k>0 \text{ odd} $$

Terms with $i<\tfrac{k-1}{2}$ will result in $\infty$ in the denominator because factorials of negative integers are infinite.  That effectively modifies the bottom limit:

$$ b_k = \sum\limits_{i=(k-1)/2}^{\lfloor (N-1)/2 \rfloor} a_{2i+1}\,\left(\tfrac{A}{2}\right)^{2i+1} \frac{2(2i+1)!}{(i-\tfrac{k-1}{2})!(i+\tfrac{k+1}{2})!} \qquad k>0 \text{ odd} $$

.

(i wouldn't mind if someone else checked my math.  these results tell us something about how a polynomial waveshaping function will sound, at least what the amplitude of harmonics will be.)


Finding peaks in an Audio spectrogram I'm using gist.github.com/sixtenbe/1178136#file-peakdetect-py to find maximum peaks in my spectrogram.

I used

_max, _min = peakdetect.peakdetect(array, lookahead=300, delta=0.3)

I get the following peaks. x_axis = time, y_axis = frequency

As you can see in the image, the peaks are not evenly spread out. All the points are skewed towards frequency axis. My array is in this format: [1, 2, 3 , 0 , 1, 3, 4]. It's 2D time+frequency array flattened into 1D array.

What can I do to show the image like this where the peaks are spread out evenly.

Try flattening your 2D array along the Y or frequency axis instead of along the X or time axis before doing the peak detect.


Retrieve frequency value of an 2D frequency spectrum How do you retrieve the values of the frequency component occurring in a 2D freq spectrum?

I have this, and if I wanted to know the values of the white dots.  what should i do ?

I apolgize beforehand for asking this question.

You have a signal in a 2D-frequency domain. You know the maximum frequencies in both directions, I assume. Hence, you know which bandwidth is represented by one pixel and then... just count the pixels from the center to the point of interest along the axes and your done. Depending on your application, taking the Euclidean distance from the center to the peak is also an option, which does end in one instead of two frequencies.

If you would like to do a more sophisticated approach, you could fit a Gaussian distribution (or another function, if you have one in mind that is more appropriate in your case) to the peak. That way you could get a sub-pixel accuracy of the location of the peak.


How to choose a phase for the deconvolution of an autocorrelation? Say I have a function, $C=C\left(x\right)$, whose fourier transform is denoted by $c=c\left(k\right)$, i.e. $C\left(x\right)=\sum_{k=-\infty}^{\infty}c\left(k\right)\chi\left(x\right)$, where $\chi\left(x\right)$ is some complex valued basis function.

I know that the autocorrelation is given by $M\left(\Delta x\right)=\sum_{k=-\infty}^{\infty}m\left(k\right)\chi\left(\Delta x\right)$, where $\Delta x = x_2-x_1$ and by the convolution theorem I know that $m\left(k\right)=\left(c\left(k\right)\right)^*c\left(k\right)$. This means that knowing the coefficients of the autocorrelation immediately gives me the squared magnitude of the original coefficients. This implies that any two functions $C'\left(x\right)=C\left(x-x'\right)$, that differ only by a translation, will possess the same autocorrelation, $M\left(\Delta x\right)$, and so if I am trying to obtain the $c\left(k\right)$ from the $m\left(k\right)$, then I can choose an arbitrary phase for the $c\left(k\right)$. However, I don't think I can choose the phase of $c\left(k\right)$ independently for each $k$, it seems like there needs to be some way to consistently choose the phase for all $k$ (which corresponds to specifying a translation for $C\left(x\right)$).

My question is: How do I choose a consistent phase for each of the coefficients for the original function, i.e. for all $k$?

For a real-valued signal $C$, the autocorrelation function $M$ is a real-valued even
function and the power spectral density $m$ (Fourier transform of $M$) is a real-valued
nonnegative even function.  Now, $m(k) = c(k)c^*(k)$ where $c$ is the Fourier transform
of $C$, as you correctly assert, but given only $m$ and no other
information about $c$ (or $C$),
it is not possible to determine $c$ at all.  Note that $m(k)$ also equals $[c(k)e^{j\theta_k}][c(k)e^{j\theta_k}]^* = c(k)c^*(k) e^{j\theta_k}e^{-j\theta_k}$ for
arbitrary choice of $\theta_k$, and the inverse Fourier transform of $c(k)e^{j\theta_k}$
is also a real-valued signal if we are careful to preserve conjugate symmetry so
that we choose $\theta_{-k}$ as $-\theta_k$. In other words, there are infinitely many
quite different signals that share the same autocorrelation function, and it is not
the case that only $C$ and its time-delayed versions have autocorrelation function $M$.
As an example, all PN signals, that is, pulse trains generated from a maximal-length
linear feedback shift register (LFSR), have an inverted thumbtack autocorrelation,
and the distinct PN signals (from different LFSRs)
are even two-level signals, a property not shared
by all of the infinitely many real signals that have inverted
thumbtack autocorrelation.


Determining how similar a signal is to a sine wave HeartMath Institute calculates a coherence score for heart rate variability by giving a score to how sinusoidal wave-like the instantaneous heart rate signal is.

http://store.heartmath.org/s.nl/ctype.KB/it.I/id.585/KB.623/.f

What DSP techniques can be used to calculate how sinusoidal a signal is, in a manner similar to score used above?

Assuming that the description I made of the FreezeFramer process is correct, then I would try determining the SINAD rating of the peak frequency in the interpolated data.

You have a sampling rate of 2Hz, giving a Nyquist of 1Hz.  You will be looking for a sinusoid with a frequency well below 1Hz.  The SINAD value is (basically) the amplitude of the peak frequency divided by the sum of all the other frequency amplitudes.  There are better descriptions available, as well as implementations for many of the more common DSP frameworks.


Reduce Signal Size to Compare Them I have multiple experiments and each of them produce several ($k$ for example) binary signals; some artificial example next:

I have a metric to compare experiment results but I need vectors of equal size to do it.

The problem is that signal length can differ from one experiment to another so I need to reduce the signal dimension somehow.

What I tried/ideas:

Calculate area under each signal. Result is vector of dimension $k$. The problem is that this method cannot differentiate signals with the same area.
Split signal on $m$ bins and calculate an area in each bin. Result is a vector of dimension $k \times m$.

Take Fourier series coefficients. Result is vector of dimension $(a+b)\times k$. The problem is that the spectrum is too large and Fourier cannot work nice (am I correct?)
Some wavelet transforms (like Haar for example). But I don't understand how to do it correctly.

I will be grateful for advices or any method that can help me.

It seems your needs matches the Dynamic Time Warping algorithm.
You should try it as a metric to compare them.

The idea of changing the size is trickier as you may loose / change data which is important for your comparison.

Maybe what you should do is alter the "Distance" metric in the Dynamic Time Warping algorithm to something which suits your needs.


Why use a 1-bit ADC in a Sigma Delta Modulator? When looking at the discrete model of a Sigma-Delta Modulator as shown below, we can see that the quantizer is modelled as a white-noise source $e[n]$. From this model, we can derive the noise shaping property of the modulator.

$$
Y(z) = X(z) + (1-z^{-1})E(z)
$$

In a typical Sigma-Delta Modulator, the quantizer is realized by a 1-bit ADC, i.e. a comparator. However, the discrete model does not require 1-bit. So why only use an ADC with a resolution of one bit, and not e.g. two or three bit? The noise is shaped to be mostly present at higher frequencies, yes, but reducing the noise even further by choosing an ADC with greater resolution is an advantage anyhow.

First of all, because it's easy to build a 1-bit ADC. It's a comparator. It's literally the easiest ADC you can build. The $\Delta\Sigma$ ADC was invented (or, rather, published) in 1962¹ !
The 2-bit ADC is more than twice as complex as that,  you need some window decision: so if you have the choice of making your 1-bit ADC run faster or building a somewhat exact 2-bit ADC, there's a solid chance you'd go for the 1-bit ADC, simply because there's less analog semiconductor design to do, and, more importantly, to go wrong!
Imagine the 1-bit ADC doing a sign decision ("is the analog voltage > 0V?"). No matter how you scale the analog voltage (multiply it with a factor $\alpha$), the result will always be the same.
Thus, the 1-bit ADC-based $\Delta\Sigma$-ADC can be scaled to fit any range simply by changing (scaling) the threshold (which effectively is the analog part of your 1-bit ADC) in the forward chain (picture from wikipedia):

It's absolutely not obvious how you'd change above diagram to fit a 2-bit ADC! But, above diagram is very close to the original publication¹, so abstraction yields that the threshold is just a single-bit quantizer.
However, it's been done, there's multiple patents² covering higher-quantization  methods. You'd have to limit the bandwidth of the forward chain³, and, depending on what you want to achieve (class-D amplifier, faster ADC, lower power consumption…), you'd then have a bit of clever logic that translates the ADC output into different kinds of feedback "pulses" and different kinds of counter increments.
So, your original claim "we don't use more than one bit in the quantization within the loop of a $\Delta\Sigma$ ADC is wrong; it is only right for the very classical implementation of the original inventors.
To totally disprove your point: There's a lot of $\Delta\Sigma$ ADCs on the market that actually do higher-order quantization; for example, the Analog Devices AD9267 uses a nine-level quantizer⁴ (3.something bit, which was a pretty impressive feat on its own to integrate, considering the speed of 640 MS/s of that quantizer).

¹  Inose, H., Yasuda, Y. and Murakami, J., 1962. A telemetering system by Code Modulation – $\Delta$-$\Sigma$ Modulation. IRE Transactions on Space Electronics and Telemetry, (3), pp.204-209.
² WO2008028142, which already cites multi-bit feedback as "prior art".
³ otherwise, you'd always only get the lowest ADC output before you could get any of the higher ADC outputs, and with infinite bandwidth, that would mean you'd already emit a pulse to "reset" the integrator, and you'd never use any but 1-bit of the ADC. That bandwidth-limiting typically happens "involuntarily" because electronics are limited by physics.
⁴ AD9267 Datasheet, p. 13, "Theory of operation"


Double check z-transform ROC of $a^nu[n]$ For the z-transform ROC of signal $a^nu[n]$, it has been computed to be $|z|>a$. For example (as I have found on Wikipedia), the signal $(\frac{1}{2})^nu[n]$'s ROC will be $|z|>\frac{1}{2}$, as
$$\sum_{n=0}^{\infty}(\frac{1}{2})^nz^{-n}=\sum_{n=0}^{\infty}(\frac{0.5}{z})^n=\frac{1}{1-0.5(z)^{-1}}$$
To this point, it said for such equality to hold, $|z|>0.5$. But I am thinking is it possible that $|z|$ can also be less than $0.5$, as long as $z\neq1$ and so the denominator won't go to $0$. However, I then realize that if $|z|$ is less than $0.5$, this term $\sum_{n=0}^{\infty}|(\frac{0.5}{z})|^n$ will go to $\infty$. Could I know if my "realization" is correct: $|z|$ cannot be less than $a$ simply because the term $\sum_{n=0}^{\infty}|(\frac{a}{z})|^n$ will go to $\infty$. Thanks in advance.

You're right, the sum doesn't converge for values $z$ with $|z|\le \frac12$, that's why the equality only holds for $z$ inside the region of convergence ($|z|>\frac12$).


signal reconstruction error in compressed sensing Does the signal reconstruction error in compressed sensing using $l_1$ norm minimization depends on the amplitude of non-zero coefficients and their location ?

Yes and No. If you consider your signal, noise free, meaning that all of the zero elements are really zero, then your answer is No, however though if the input noise is taken into the account, then the higher amplitude of the non-zero elements results in better signal reconstruction (high input SNR gives high output SNR).


rtlsdr sample rates and nyquist rate I have a rtlsdr dongle which apparently has a max sample rate of about 2.4 MHz.  I'm wondering how it's possible to capture signals at say 900 MHz with this device.  Doesn't the sample rate need to be 2x900 or 1.8 GHz to capture without aliasing?

Sorry, I'm so new to this stuff!

Thanks!

The RTL-SDR, as most software-defined radios, has an analog front-end that includes a tunable mixer. This mixer can down-convert the frequency band you're interested in to baseband, where it is sampled.

For example, let's say that you want to receive a signal centered around 900 MHz and with bandwidth $2B$ -- that is, the band covers the range from $900-B$ to $900+B$ MHz. You tune the mixer to 900 MHz and (after filtering) you end up with a signal centered around DC, covering the band from $-B$ to $B$ MHz. In general, this is a quadrature signal, that can be sampled at rate $2B$ complex samples per second.

As you see, with this setup, the sampling rate of the device only limits the bandwidth it can capture, without regard for the value of the carrier frequency (within the limits of the filters and mixers, of course).


resources for learning more about pedestrian tracking I would like to build a computer vision application capable of doing robust pedestrian tracking in potentially crowded environments from a single video source. To this end I would like to find resources in order to learn as much about the current state of the art as possible. This could be:

recent survey papers
papers outlining novel approaches that have been demonstrated with success
open-source applications or algorithms
online courses or even course materials such as slide decks comparing different approaches
...

Is there anything along these lines that you could point me to? Thanks.

What I have so far:

A Survey of Human-Sensing: Methods for Detecting Presence, Count, Location, Track, and Identity; Teixeira, Dublon and Savvides
I just found the question People Detection and Tracking on SO which also has some good leads.


Modelling an existing DSP filter by analysing the differences between its inputs and outputs A DSP filter produces always the same output.

$$\boxed{\textrm{wav file}}{\longrightarrow}\boxed{\textrm{DSP}}{\longrightarrow}\boxed{\textrm{wav file(modified)}}$$

The DSP is a black box as I don't have access to its code, but I have access to original and modified wav files.

Questions:

Can the DSP filter be mathematically modelled (say in MATLAB) by comparing the original and modified wav files?
If so, would the resulting modelled DSP, when applied to any wav file, generate a wav equals to the one that will be produced by the original DSP?

If I understood correctly, you have a signal $x[n]$ and an unknown discrete-time LTI filter, so that you can look at the output of the filter $y[n]$. Now you are looking for the impulse response of the filter $h[n]$.

The output follows a convolution rule
$$y[n]=x[n]*h[n]=\sum_{m=-\infty}^{\infty}x[m]h[n-m]$$
The process to find $h[n]$, given $x[n]$ and $y[n]$ is called deconvolution. There are some approaches to do that.

For example, if the Fourier transform of the signals exist, using convolution theorem, the convolution can be converted to product in the frequency domain: $$Y(\omega)=X(\omega)H(\omega)$$
Therefore, $$H(\omega)=\frac{Y(\omega)}{X(\omega)}$$
You should however be careful about the zeros of $X(\omega)$. If $x[n]$ can be an option, then you should choose it such that it can excite the system at all frequencies such that you can probe the desired frequency range that you are interested in $H(\omega)$. Alternatively, you can excite the filter in repeated number of times at different frequencies, which can be interpreted as sampling the frequency response. This will give you a sampled spectrum, so if you know the filter does not have singularities between the sample frequencies, then you can model it by an LTI filter with your desired order (depending on your acceptable approximation error).

In response to your second question, after you acquired an approximation such as $\hat{h}[n]$ for the impulse response, then you can apply it to any desired input by using convolution.


Scientific names to piano key numbers Is there a formula that converts scientific note names to piano key numbers? As an example, given A4 it would output 49.

The frequencies are on an exponential scale.  You can use any note as your base value.  On a Chromatic scale, each note is a twelfth step of an octave (doubling of frequency, same tone).  So, a very common formula is to use the 440Hz A note.  (Search on 432Hz to see a controversy about this.)

$$ f = 440 \cdot 2^{\frac{n}{12}} $$

Where $n$ is the semi-tone count from your base.

Here is a good reference:  Note names, MIDI numbers and frequencies

Upon rereading, I didn't really answer your question.  Because of the gaps in note names (black keys on the piano), this is more of a programming exercise than a simple formula.

In Python:

#=======================================================
def NoteConverter( ArgNoteName ):

        theLetter = ArgNoteName[0:1]
        theNumber = int( ArgNoteName[1:] )

        theOrdinal = "C D EF G A B".find( theLetter )

        theMidiNumber = theNumber * 12 + theOrdinal + 12

        return theMidiNumber

#=======================================================

This is for Midi numbering.  Subtract 20 to get your piano key number (on most pianos).


Removing Glare from Image Disclaimer: I am in NO way engaged in signal processing.  Just very curious...

Ny question is, very simply: is it possible to apply software processing to an image to remove or reduce headlight glare?

It would be a picture such as this:

If this IS possible, is there a software library available for this?  I don't care which language, I just want the capabilities.  Slash does anyone know if there is current research on this?

The problem you might have is that when bright lights are present, the camera will be using a fast shutter speed to keep the image from saturating too much.  All the detail you are interested in is in the bottom few bits of the data.

So even if you remove the glare and then try and pull the rest of the image up in level, it's very noisy.  For example, if you take your image and play with the levels to pull up the darker sections and just allow the glare to clip, this is the result:

If that's OK for onward processing (I don't know what you have in mind) then you could maybe remove the glare by looking at distributions of saturated pixels, and expand into the "glow" around them.

A better solution (if you have the option) can be to use a high-dynamic range camera which has an non-linear response in the pixels, and also often 10 or 12 bits per pixel of usable resolution, which means you can keep 7-8 bits for the darker sections.


Which information do we get from magnitude and phase spectrum? I am learning image processing. I want to ask very basic question related to FFT topic

  Which information do we actually get from "phase spectrum" and "magnitude spectrum" about an image?

You get the exact same information you will receive if you analyse a 1-D signal using the Fourier analysis tools, for example. To illustrate this consider the following examples

We perform Fourier Transform on it and obtain the following spectrum

As you can see there are two symmetric dots representing the frequency that is present in the image and the central, DC, component.
Next up

Again the Fourier Transform

It is obvious that the higher frequencies are further away from the DC, central, component. To understand the logic behind this you shall remember that we are actually working with complex numbers and their properties kick in.
To answer your second question: You are usually working with sparse data structures (vectors, arrays, tensors, etc.) in the frequency domain. Now you tell me which is easier: Working on the original data set consisting of 500 different values or the transformed one with half of them 0s ?
Now, the importance of phase - drum rolls - it is the same, meaning HUGE
Again, a sample image

Fourier Magnitude Spectrum

Fourier Phase Spectrum

Now, we inverse the transformation by using just the amplitude and then just the phase information
Inverted amplitude spectrum

Inverted phase spectrum

I think you get the idea :D
Now, here are some sample images and their respective amplitude and phase spectra, so you can practice a wee bit


Rate offset and jitter concealment for streaming audio Given streaming audio samples between two separated devices with separate non-synchronized but almost identical sample rates (say both 44100+-2ppm), what are some methods to best conceal any difference and/or jitter in rates from causing audible artifacts in the audio output, and with minimal required added latency for the concealment algorithm?

at +-2ppm, we'd win or lose a single sample any 45 seconds – my pragmatic approach here would be ignoring the problem, which, of course, has extremely low latency :)

But yes, especially if the producer is a bit slower than the consumer, you'll have a problem in the long run.

Now, there's no easy solution to this – first of all, your receiver needs to be able to see that problem at all.

But let's assume it does – simply by counting (recovered?) clock phases of the receive signal and comparing to its own clock cycles (this gets easier for serial audio data buses, e.g. I2S, because there's many bits per sample, so your bit clock is higher, so the counting offset gets visible much faster, but only if the transmitter continously clocks the bus).
If that can't work for technical reasons (non-continuous bursting of audio data, bus clock is decoupled from ADC clock, no access to recovered clock, etc), the only way I could think of would be the one of keeping track of buffer levels – I assume your receiver does have one, albeit small, buffer between its digital input and the ADC, because otherwise, the ADC would have to run directly off the bus clock, right?

You could observe the buffer "fillage" or clock offset, low pass filter and or caclulate a derivative of that info, define upper and lower boundaries, and simply drop a whole sample if you've got a next-to-critical source-faster-than-sink mismatch, or add a zero sample in the opposite case.

Of course, the "right" way to do this would be resampling – but from your question it seemed to me you're more interested in latency than in accuracy; you use the same frequency offset info as gotten from either bus/clock shift or buffer fillage to drive your own arbitrary resampler:

For example, what's pretty popular among these parts (PC-style computer based SDR), is a resampler based on MMSE-combined filters. They're not overly long.
The next more elegant thing is the Polyphase filterbank arbitrary resampler: again, filterbank with different filters, this time actually used just for a rational interpolator/decimator, but followed by logic that "combines" the output of the next lower rational resampler "branch" with the next higher linearly to get the wanted real resampling ratio. Advantage is that with a big filterbank, your error gets pretty small at reasonable CPU cost, and latency is bound by the length of the rational resampler bank filter - which typically won't be hundreds of samples long :)


Power and Energy of an Audio file I have a question about calculating the power and Energy of an audio recording.

I have imported a waveform of an audio recording in WAV format into Matlab and then calculated the energy and power in the time and frequency domain. Using Parseval's theorem, the energy in time and frequency is the same.

 clear all; close all;

 %Wav File einlesen und die Daten in data und Abtastrate in Fs
 [data,Fs]=audioread('120A_v20_05s.wav');  %Ermittlung der Arraygröße
 [nSamples,nChannels]=size(data); %Länge der Audiodatei ermitteln
 waveFileLength=nSamples/Fs;

 t=[0:length(data)-1] / Fs;

 %Signal darstellen figure
 plot(t,data)
 ylim([-0.3 0.3])
 title('120A v20')
 grid on xlabel('Zeit')
 label('X(t)')

 E1_timedomain=sum(abs(data.^2)); % Energy time domain

 L=length(data);
 Ptime=(norm(data)^2)/L;  % Power time domain

 y = fft(data);

 Pfrequ=sum(y.*conj(y))/(L^2); %Compute power with proper scaling.
 E1_frequdomain=sum(abs(y.^2))/nSamples; % Energy frequency domain

E1_timedomain = 4.003
Ptime = 1.8154e-04

My question is which units have the energy and power in this case? Joule and Watts or do I have to consider a factor?

Many Thanks!
Greetings Mathias

In principal, you get $\text{Ws}=\text{J}$ for energy and $\text{W}$ for power, but as you are calculating these measures from a wav-file, you do not get any unit. What you get is just a ratio to the maximum possible value of a wav-file. A square wave with amplitude one has a level (power) of $0\text{ dBFs}$, which means dB Fullscale. Your signal has roughly $-75\text{ dBFs}$, so that would be a valid result. If you want to get a result in actual physical units, you will need to know several things:

the sensitivity of the sensor (e.g. a microphone), that was used for audio recording. This is usually given in $\text{dBV/Pa}$
the calibration of the A/D converter used, usually given in form of the overload point in $\text{dBV}$, meaning that a sine wave with this level just drives the ADC to fullscale ($0\text{ dBFs}$), any increasing of the voltage leading to clipping.

With this information, you could calculate back the real power the actual sound event had at the point of the microphone.


Sensitiviy and Resolution in Windowing Functions When I read the  article Windowing Functions, there are two terms that are does not understanding them well, sensitiviy and resolution.

I looked for this two terms on the web, And I conclude that sensitivity are how much you can distinguish the components of the signal, and resolution how many components you can have for a signal. Is that correct?

If my above conclusion was correct I see that the more resolution the more components I have, yields less distinguishing of that components, which means the more resolution the less sensitivity, does that also is correct? Also, does the overlpping windows is a solution for the little resolution that some windows has, for example hanning window?

I object to the Wiki author's use of the words "resolution" and "sensitivity." Those words are vague and the author does not define or explain them in any meaningful way. If I was forced to use the author's poor terminology I'd say, "A fine-resolution window is a window having a narrow main lobe. (And windows whose main lobe widths are narrow have high-level sidelobes.) A window with good sensitivity is a window with low-level sidelobes. (And windows with low-level sidelobes have wider main lobes.)

I do not understand the 2nd and 3rd paragraphs of your question.

A window with a narrow main lobe improves your ability to distinguish (recognize, discriminate), in the frequency domain, two closely spaced spectral sinusoids.  A window with low sidelobe levels improves your ability to detect, in the frequency domain, a low-level sine wave that is close in frequency to a high-level sinusoid.


Location of Matched Filter In  a digital communication receiver, matched filter(SRRC filter) is used. In different literatures, the location of the matched filter is mentioned in two different places. One is just after the ADC followed by the resampler of a timing recovery system, and the other one is after the resampler of timing recovery process followed by a timing error detector (TED).

To my understanding, practically it should be after the resampler. Lets say, the interpolation factor is $4$ at the transmitter side, that means $4$ samples per symbol. So, if the Symbol rate is $2\textrm{ MHz}$, the sample rate is $8\textrm{ MHz}$.

Now, in the receiver side, let's say the ADC I am using is sampling at $40\textrm{ MSps}$ i.e the sample rate is $40\textrm{ MHz}$, which is $20$ times the symbol rate, or $5$ times the transmitter sample rate. It is only after resampling that we will get back our original transmitter sampling rate, so that SRRC filter (with oversampling rate $4$) can be employed. Otherwise, if we use before the resampler, we have to use a SRRC filter (with oversampling rate $20$), which will increase the number of filter coefficients and hence the hardware cost.

Is my understanding correct?

The time span of the zero crossings increases after the final RRC filtering (and the symbol sampling locations converge which is the goal for the benefit of zero ISI but the zero crossing increase in the process is to the detriment of timing recovery!). So if you are using a Gardner TED which is sensitive to this, it is better to have TED prior to RRC filtering, as the timing SNR will be higher.

However a synchronizer such as the Mueller and Mueller which operates on symbol decisions (1 sample per symbol) has better performance after the RRC filter.

Below are details related to the use of the Gardner TED, to show the considerations involved:

See my figures below from when I studied this more closely with the Gardner TED, showing also the trade-off to be made with a lower roll-off (alpha value) for better spectrum efficiency but lower timing SNR due to reduced discriminator slope and higher pattern noise. Shown is the "Mean(TED)" which is the timing discriminator for a Gardner TED (higher slope means higher loop gain/sensitivity), and the pattern noise of the detector versus timing offset. Notice that a higher timing SNR is achieved if you do your timing error detection prior to RRC filtering.  This is because the time span of the zero crossings increase after RRC filtering, while the symbol decision locations converge (see figure of waveform before and after RRC).

For comparison, if desired, I can also later include the M&M synchronizer since I had studied that in similar fashion as well, but for now I bottom lined the comparisons I derived at the end of this post.

Gardner TED for QPSK / QAM:

$TED = I_{2n}(I_{2n+1}-I_{2n-1}+Q_{2n}(Q_{2n+1}-Q_{2n-1}) = real[conj(y_n)(y_{2n+1}-y_{2n-1})]$

This shows the frequency spectrum characteristics of the pattern noise from the Garndner Timing Error Detector:

For comparison (and consideration), below is the Mueller & Mueller (M&M) Synchronizer:

Without elaborating on the Mueller & Mueller synchronizer (M&M) in similar detail, I will at least include below the key take-aways I made when comparing the Gardner TED vs the Mueller & Mueller synchronizer:


Erosion compotition rule $\left ( f \ominus g \right ) \ominus h = f \ominus \left (g \oplus h \right ) $ I want to prove that erosion of a signal by a function $g$ followed by erosion with another function $h$ is equivalent to erosion of the signal by the dilation of the two functions (Erosion property):

$$\left ( f \ominus g  \right ) \ominus h = f \ominus \left (g \oplus h  \right ) $$

My approach:

By definition of erosion we get that: $$\left ( f \ominus g\right )(x) = ⋀_{y} f(x+y)-g(y)$$

So, we get that $$\left ( f \ominus g  \right ) \ominus h = ⋀_{y}\left ( ⋀_{z} f(x+y+z)-g(y)-h(z)\right )$$

Is my approach correct so far? And if yes, how do I go on?

Morphological Erosion is not associative as dilation. Much like addition is associative but not subtraction. Consider the following basic definitions first:$$(A \ominus B)^C = A^C \oplus \widetilde{B} \tag{1}$$
$$(A \oplus B)^C = A^C \ominus \widetilde{B} \tag{2}$$

where $\widetilde{B}$ is the reflection of kernel $B$ and $A^C$ is the compliment of image $A$, inverse of image $A$. Now let us evaluate the expression in the question:
$$(A \ominus B) \ominus C = ((A\ominus B)^C \oplus \widetilde{C})^C \tag{using eq. 1}$$

$$(A \ominus B) \ominus C = ((A^C\oplus \widetilde{B}) \oplus \widetilde{C})^C \tag{using eq. 1} $$

$$(A \ominus B) \ominus C = (A^C \oplus (\widetilde{B} \oplus \widetilde{C}))^C \tag{since dilation is associative}$$

$$(A \ominus B) \ominus C = A \ominus \widetilde{(\widetilde{B} \oplus \widetilde{C})} \tag{using eq. 2}$$

$$(A \ominus B) \ominus C = A \ominus (B \oplus C) \tag{distributing the reflection}$$

Hope that answers the question


Changing Units of Axis I'm writing a Code for STFT in matlab and I am in the initial stage. Right now I have managed to do the FFT of a small chunk of a time signal. I used zero padding.

The problem I am facing is that after FFT there are 3500 bin frequencies. (In total there were 7000, I took half of them). My sampling frequency was 1000. So that equates 1 Bin = 10/35 Hz.

When I plot the data, I have the bin values in the X-Axis. In this was my Signal seems to have amplitudes in the 421st Bin or 60Hz.

How can I make my X-Axis such that it progresses by 10/35 and shows the spike in the correct place at 60 instead of 421. Thanks.

Thanks for the clarification.

If you're currently doing

plot(y)

then you should do

plot(x,y)

where

x = linspace(0,1000,length(y));


Likelihood of Unscented Kalman filter Short version:

How to calculate a likelihood of Unscented Kalman filter?

Long version:

Likelihood for linear Kalman filter (KF) is:

$$\mathcal{L} = \frac{1}{\sqrt{2\pi S}}\exp \left[-\frac{1}{2}\mathbf y^\mathsf T\mathbf S^{-1}\mathbf y\right]$$

Other words likelihood of KF depends on residual and uncertainty:

$$\mathcal{L} (\mathbf y, \mathbf S)$$

We can see that uncertainty depends on process noise $\mathbf Q$ which is function of a time $\mathbf t$:

$$\mathbf{S} (\mathbf P) \Rightarrow \mathbf{P} (\mathbf Q) \Rightarrow\mathbf{Q} (\mathbf t)$$

This allows conclude that likelihood of KF depends on time:

$$\mathcal{L} (\mathbf y, \mathbf t)$$

For Unscented Kalman filter (UKF) instead of uncertainty $\mathbf S$ is used covariance of sigma-points in measurement space $\mathbf P_z$, where:

$$\mathbf P_z = \sum w^c {\left(\mathcal Z - \mu_z\right)\left(\mathcal Z - \mu_z\right)^\mathsf{T} + \mathbf R}$$

In this case the likelihood can be found as:

$$\mathcal{L} = \frac{1}{\sqrt{2\pi P_z}}\exp \left[-\frac{1}{2}\mathbf y^\mathsf T\mathbf P_z^{-1}\mathbf y\right]$$

But, neither $\mathcal Z$ no $\mathbf R$ depends on process noise $\mathbf Q$, so likelihood of UKF doesn't depend on time.

In practice it means the following: if I try to find likelihood of KF consequentially increasing time which passed after last measurement then likelihood will be reduce, but for UKF likelihood will not be changed.

Where is my mistake?

Your mistake is that $\mathcal Z$ does depend on process noise. It's just a bit more obscure than the linear filter. $\mathcal Z$ is the projection of the sigma points through the observation function. The sigma points are a $\mathcal f(P)$. Therefore $\mathcal Z$ is a $\mathbf f(P)$, and therefore related to $\mathbf Q$.


Can the (point-wise) multiplication of two discrete band-limited functions be inherently aliased? I'm not fully used to the language of DSP, so forgive me if I confuse some of the notation, or if the below is trivial.

I'm interested in the case where we want to do the point-wise multiplication of two discrete signals, $f[x_n]=g[x_n]h[x_n]$ for $n=0,\dots,M-1$. This is easy enough to do in the computer. However, I'll take a detour and formulate the problem in the (Discrete) Fourier domain, where the multiplication becomes a convolution:
$$
\widehat{f}[k_n] = \sum_{j=-M}^M \widehat{g}[k_j]\widehat{h}[k_n-k_j].
$$
Interesting is that (just like the MATLAB conv command) we get $2M-1$ outputs; i.e., the convolved signal is inherently longer than the original signal. If we take the result back to the sample domain, we have to throw away $M$ samples to get back to the original domain (or get back to sample domain and pick every second entry).

I ran the following test in MATLAB that seems to confirm my thought,

M=100;
x=[0:M-1];
g=sin(x/M*6*pi);
h=cos(x/M*6*pi);
f=g.*h;

ff = fft(f);gf = fft(g);hf = fft(h);
v  = ifftshift( (-1/2:1/M:1/2-1/M) + mod(M,2)/(2*M) ); % Normalized frequencies

figure(1); plot( x, g, x, h, x, f )
title('Discrete multiplication')
legend('g','h','f')

figure(2); plot( v, abs(gf), v, abs(hf), v, abs(ff) )
title('Spectrum')

This generates the following two figures:

and

that appears to confirm my thoughts. The spectrum of $f[x_n]$ (in yellow) is broader than the spectrum of either input signals. Then to get to the question in the title: (1) is it possible that the point-wise multiplication of two discrete, band-limited functions is aliased itself? And (2) although I was able to recreate the $f$ also using the convolution theorem using

ghf = conv( fftshift(gf), fftshift(hf)./M, 'same' );
gh  = ifft(ifftshift(ghf),'symmetric'); % = identical to the vector 'f'

would it be possible to circumvent the aliasing by outputing the 'full' convolution, thereby getting $2M-1$ samples back, rather than $M$ samples?

(1) is it possible that the point-wise multiplication of two discrete, band-limited functions is aliased itself?

Yes.  Or at least it's possible that the result won't fit in a tidy manner into your sampling rate.  Because multiplying two signals together creates energy at different frequencies than the frequencies of the original (this is easiest to do by a single-frequency example: $2 \cos \omega_1 t \cdot \cos \omega_2 t = \cos (\omega_1 - \omega_2)t + \cos(\omega_1 + \omega_2)t$.  If that $\omega_1 + \omega_2$ is to great to fit into your sampling rate, then it'll alias.

  (2) would it be possible to circumvent the aliasing by outputing the 'full' convolution, thereby getting $2−1$ samples back?

Yes and no.  Yes, you could, but to be consistent the result would have to be at a sample rate $\frac{2M - 1}{M}$ times the original.

An equivalent operation would be to upsample both signals by a factor of two and then multiply -- if the upsampling were done in a way that results in correctly band limited signals, once you multiplied them they would fit into the new, bigger bandwidth made available to them by upsampling.


Can a unit delay in discrete time be represented by exponential functions? If we have a signal $y[n]$ and its unit delayed version $y[n-1]$, can we write $y[n-1]$ in terms of $y[n]$ times some exponential?

The reason I want to do this is to then take $y[n]$ common and simplify my analysis that I am doing for my research. The above is a generic form.

In continuous time, we can, though it is a mathematical curiosity of little practical use. Suppose we have a smooth continuous-time signal $x$. Using the Taylor expansion,

$$x (t-1) = x(t) - \dot x (t) + \frac 12 \ddot x (t) - \cdots = \left( 1 - \frac{\mathrm d}{\mathrm d t} + \frac 12 \left(\frac{\mathrm d}{\mathrm d t}\right)^2 - \cdots \right) x(t) = \color{blue}{\exp \left( - \frac{\mathrm d}{\mathrm d t} \right) \, x (t)}$$


Minimum sampling frequency, quantization, and bitrate calculation An analogue sensor has a bandwidth which extends from very low frequencies up to a maximum of 14.5 kHz. Using the Sampling Theorem what is the minimum sampling rate (number of samples per second) required to convert the sensor signal into a digital representation?
If each sample is now quantized into $2048$ levels, what will be the resulting transmitted bitrate in kbps?
EDIT
The answer I get is 11 bits per sample, $2048 = 2^{11}$, is this calculation correct?

It is not clear in which part of the problem statement you experience issues.
However, a few hints may help you get the gist of such problems

The sampling frequency, it is intuitive that in order for a system to monitor the state of another system it must poll it at higher rates than the highest rate in which the monitored signal is able to change its state. But what is the minimum sampling frequency that gives our sampling system the ability to correctly approximate the continuous signal in the discrete time domain?

This frequency threshold is at two times the highest frequency of the continuous signal:
$$
f_{sampling} \gt 2\cdot f_{bandwidth}
$$
The theorem that the problem statement suggests is the Shannon-Nyquist sampling theorem, the theorem that the above inequality comes from.
In real sampling systems, the sampling frequency is never taken exactly on this threshold, usually the threshold is taken a little bit higher than $2\cdot f_{bandwidth}$, for example at $2.2\cdot f_{bandwidth}$.
However, in theory it is safe to assume a sampling frequency at the Nyquist rate (i.e. the aforementioned frequency threshold).
Your analogue signal's highest frequency is $14.5\, Hz$, ranging from very low frequencies so it is safe to assume, this actually the bandwidth frequency.
All the above summarized in the following equation:
$$ f_{sampling} = 2\cdot f_{bandwidth}$$

Bitrate, in order to calculate the bitrate of your sampling system it is necessary to know the sampling frequency. Fortunately, this was the first the step of the exercise. The second component for the bitrate calculation is the actual length of each sample in bits. The minimum number of bits that are necessary in order to represent a quantized value in x distinct levels is given by:

$$
2^{bitlength}\, =\, x
$$
or
$$
bitlength = \lceil\log_2{x}\rceil
$$
You could actually deduce the bitrate equation by performing a dimensional analysis, because it is stated that it is measured at kbps(kilobits per second). For instance, we know that $Hz$ is measured in $\frac{samples}{sec}$ and kb stands for kilobits meaning $1000\, \frac{bits}{sample}$, so:
$$
bitrate = bitlength\cdot f_{sampling} \left[\dfrac{bits}{sample}\dfrac{samples}{sec} \right]
$$
or in kbps
$$
bitrate = bitlength\cdot f_{sampling}\cdot 10^{-3} \left[kbps\right]
$$


Is spectral density sometimes normalized by sampling rate rather than bin size? I'm a scientist conducting an experiment that requires some signal processing. My expertise is not in signal processing, thus here I am. We've basically re-created an experiment conducted by other scientists, attempting to check their results. Here is a link to their paper: Ultrasensitive Inverse Weak-Value Tilt Meter
In short, a laser bounces off of some mirrors, one of which is oscillating at a controlled sinusoidal frequency, onto a quadrant detector, which outputs an electrical signal to an oscilloscope where we record it. So, you end up with a noisy record that has a tiny, known sine wave hiding in it.
Everything I've read indicates that to calculate spectral density, you must:

Get the spectrum by performing an FFT* on the record
Normalize the spectrum by the bin size, which is sampling rate divided by number of samples (Fs/N)

*For clarification, when I refer to an FFT, I'm referring to the single sided, absolute value, of the FFT, normalized by the number of sample points, N. So, we took the FFT of the signal, threw away the negative frequencies, doubled the positive frequency values (except DC and Nyquist), and divided by N. I checked this method by feeding signals directly from a function generator to the oscilloscope and verifying that the resulting peaks matched the frequency and amplitude of the inputs.
But, in the paper linked above, they seem to have normalized their spectrum by sampling rate only. I say this because at the top of the first column on page 3, they point out that the sampling rate is 1 kHz, and in the footnote on page 3, they point out that the peak in their spectral density plot (Figure 4) is 1.6 nrad / sqrt(1kHz). They make no mention of bin size or number of samples (N). Since I'm trying to directly compare my numbers to theirs, I need to know definitively what is going on here. Are there two definitions for spectral density? Thanks in advance.

The use of $rad/\sqrt{\text{Hz}}$ suggests that this is phase noise specifically (a spectral density due to phase fluctuations), and typically in my use this has been described as a power spectral density (units of $rad^2/\text{Hz}$), so this is just the square root of that quantity.
The reason the DFT (of which the FFT computes) is divided by $N$ is to normalize the FFT to be the same units of the time domain signal, specifically using the following normalized form of the DFT:
$$X_1(k) = \frac{1}{N}\sum_{n=0}^{N-1}x[n]W_N^{nk}$$
Versus the typically version that is not normalized which the FFT returns:
$$X(k) = \sum_{n=0}^{N-1}x[n]W_N^{nk}$$
With such a normalization, the magnitude of $x[n]$ at any specific frequency will match the magnitude of $X(k)$ for that frequency. For example, if we had a time domain waveform of a sinusoidal phase error versus time given as:
$$\phi[n] = A\cos(\omega n) = \frac{A}{2}e^{j\omega n} + \frac{A}{2}e^{-j\omega n} \space \text{rad}$$
Then assuming $\pm\omega$ were exactly on a bin centers (for the DFT due to its circular nature $-\omega = N-\omega$), the resulting two bins in $X_1(k)$ would have a magnitude of $\frac{A}{2}$, matching the magnitudes of the time domain waveform.
As a power spectral density (meaning we are interested in the power over a given frequency range) the normalized power of each frequency index in the DFT (aka bin) is then:
$$|X_1(k)|^2 = \frac{|X(k)|^2}{N^2} \space \frac{\text{rad}^2}{\text{bin}}$$
(Where the units of $\text{rad}^2$ for the power quantity $|X_1(k)|^2$ only make sense if x[n] was the phase noise in units of radians).
$\frac{\text{rad}^2}{\text{bin}}$ is a power quantity per bin. To make this the recognized form of the power spectral density in power/Hz we recognize that $Nd = f_s$ where $N$ is the number of samples in the DFT, $f_s$ is the sampling rate, and $d$ is the spacing of each frequency index (bin as the OP used) in Hz resulting in the spectral width of each bin in Hz:
$$d = \frac{f_s}{N} \space \frac{\text{Hz}}{\text{bin}}$$
Thus
$$ \frac{|X(k)|^2}{N^2} \frac{\text{rad}^2}{\text{bin}} \times d^{-1} \frac{\text{bin}}{\text{Hz}} = \frac{|X(k)|^2}{N^2}\frac{N}{f_s} \frac{\text{rad}^2}{\text{Hz}} = \frac{|X(k)|^2}{N f_s} \frac{\text{rad}^2}{\text{Hz}}$$
This result would specifically be what we typically notate as $\scr{L}_{\phi}(f)$ as the two-sided power spectral density due to phase fluctuations (since the DFT contains both sides of the spectrum, in contrast to the one-sided PSD which is $S_\phi(f) = 2\scr{L}_{\phi}(f)$.).
Note we say "due to phase fluctuations" since the units here were phase. It is also interesting how the phase unit in radians when squared is the power unit relative to the carrier (often expressed as dBc/Hz). This is clear for small angles given the small angle approximation $sin(\theta) \approx \theta$, or geometrically the quadrature component being the noise as phase noise relative to the in-phase component being the carrier that has been rotated due to that phase, such that the ratio of the two is the phase unit in radians, for small angles!) This is why when phase noise is dominant this computation will match the actual power measurement we see under test with a spectrum analyzer.

Further update:
The OP clarified in his comments that his question is specific to the peak at 30 Hz offset as shown in this plot:

It isn't specified but assuming this is a two-sided spectral density, the peak of a single tone would have a total power independent of density, so we would typically report its result as $\text{rad}^2$ and not $\text{rad}/\text{Hz}$ (or the magnitude quantity as the square root $\text{rad}$ as used in this plot, meaning this plot is $\sqrt{\scr{L}_{\phi}(f)}$). The paper also incorporates a moving average of 5 and suggests in a foot note that the peak would be $\approx 1.6 \text{nrad}/\sqrt{1\text{kHz}}/5$, and the plot was scaled (moved up or down) such that the level of the tone landed on this expectation.
I suggest that the peak would be at either $\approx 1.6 \text{nrad}/20$ or $\approx 1.6 \text{nrad} \sqrt{2}/20$ depending on if the spectrum is intended to be double-sided or single-sided which should be specified. The sampling rate does not change the value of the tone on the spectral density when the units are already in nrad, so there should also be no $\sqrt{1\text{kHz}}$ in that answer - The sine wave theoretically occupies zero bandwidth, or for practical reasons we can assume we integrated that power over a small bandwidth in order to measure the peak we see. Either way the density becomes a single figure for the tone independent of bandwidth. Any windowing applied in the time domain prior to the FFT (other than the rectangular window) will also shift the value of the tone differently from the values for the noise. Further details below.
To confirm that assumption, here is my prediction of where such a tone would be:
The 1.6 nrad oscillation is specified as the peak to peak value and thus is of the form:
$$\phi(t) = \frac{1.6}{2} \cos(2\pi f t) \space\space \text{nrad}$$
with $f=30e3$
If the spectrum is two-sided (as $\sqrt{\scr{L}_\phi(f)}$ rather than one-sided as $\sqrt{S_{\phi}(f)}$), then the spectrum is only showing the upper half of this two-sided spectrum, with both sides given by:
$$\phi(t) = \frac{1.6}{2} \cos(2\pi f t) = \frac{1.6}{4}e^{j 2\pi f t} + \frac{1.6}{4}e^{-j 2\pi f t} \space\space\text{nrad}$$
Thus prior to the effect of the moving average filter (MAF), I would predict the tone shown on a double-sided spectrum to be at:
$$\frac{(1.6e-9)}{4} = (4e-10) \space \text{rad}$$
Notice the units are $\text{rad}$ and not $\text{rad}/\sqrt{\text{Hz}}$ as the standard deviation of the tone itself is not a density spread across frequency, unlike that of the noise.
I assume the moving average filter that is mentioned was done on the frequency domain samples. If in the time domain there would be an additional loss of 0.963 but I don't see evidence of such a moving average response in the plot, in which case with a moving average of frequency samples, the tone is reduced by a factor of 5 as the author had done, resulting  in $(4e-10)/5 = (8e-11)$.
If the plot was supposed to be a single-sided spectrum $\sqrt{S_{\phi}(f)}$, then the result would be $\sqrt{2}$ larger or $1.13e-10$, which is consistent with the standard deviation of $\phi(t)$ reduced by the MAF.
Neither of these results match the plot, but this is where I would expect a 30 Hz tone after a moving average of 5 samples when sampled at 1 KHz if the units of the spectral density are $\text{nrad}/\sqrt{\text{Hz}}$, for either case of a single-sided or double-sided spectral density. Also note that my computation was independent of the bin size or number of samples since as the author of the paper was intending to do (and perhaps did if I made an error in my prediction) was to predict the expected value of that tone and then scale the plot accordingly. My earlier answer shows how I would scale the result from the DFT directly in which case the bin size and number of samples would be involved.
As a further note since these spectrums are being derived from FFT's and since the OP is interested ultimately in assessing noise: We must also be careful to account for the equivalent noise bandwidth due the effect of windowing especially if we are normalizing the plot based on the power of a tone. (and other effects such as scalloping loss etc which have been minimized by choosing a tone at or near a bin center as was done). Any windowing done on the time domain signal other than the rectangular window will widen the bandwidth of each bin beyond the single bin as given by the rectangular window, which means that the noise measured will be larger than the actual noise! Further the window has a loss reducing the signal from the tone and the noise, but because of the effectively wider noise bandwidth of each bin the noise will go down less than the tone (the tone only occupies one bin)! The effect of the moving average in frequency on SNR is also affected by the window since the adjacent noise bins are no longer uncorrelated. I detail this further in this post: Find the Equivalent Noise Bandwidth


normalizing STFT output and calculating energy of signal I have two questions regarding power spectra and energy of a signal:

For power spectra calculation I used STFT, and I saw that the value of STFT should be normalized before applying logarithmic function. The normalization is done by multiplying them by $2/\mbox{window length}$. What is the reason?
I also want to calculate the energy of the signal according to Parseval theorem and using STFT. Should I use the average of STFTs or the sum of them?

The STFT yields complex values. Energy conservation (up to a scaling factor) should apply to the square of its magnitude, also called spectrogram. Depending on

the initial normalization of the FFT,
the window shape,
the lag (do you compute the FFT every $h$ sample, $h$ is called the hop, sometimes),

the scaling factor will differ, but will be the same for all signals.

$2$/window length is specific to a choice of the three parameters above. A scaling factor will only result in a $\log$-magnitude shift, which does not really matter for display or relative energy.
If you can compute the energy normalization based on the three parameters above, good. If not, take a randomly-picked signal (non zero), compute the ratio between the energy of the signal and that of the spectrogram, and you have your factor. Beware of the border effect, better take a centered signal with enough zeros (above the window support) on the right and the left.


Convolution that results in an all-zero sequence I am asked to find a pair of sequences, each one of which contains three distinct values and where the convolution is an all-zero sequence.

I've came to the conclusion that one of the sequences must be infinite but just can't think of any examples...

Expanding on my comment when this question was closed, suppose that
$$x[n]=e^{j2\pi n/3},~\forall n,\quad y[n]=\begin{cases}e^{j4\pi n/3}& n=0,1,2,\\0, &\text{otherwise.}\end{cases}$$

Then,
$\displaystyle \sum_{n=-\infty}^\infty x[k−n]y[n]= \sum_{n=0}^2 e^{j2\pi (k−n)/3}e^{j4\pi n/3}
= e^{j2\pi k/3}\sum_{n=0}^2 e^{j2\pi n/3}=0.$

More generally,   $\displaystyle x[n]=e^{j2\pi n/N},~\forall n,\quad y[n]=\begin{cases}e^{j2\pi (N-n)/N}& n=0,1,2,\ldots, N-1\\0, &\text{otherwise,}\end{cases}$
works as a solution for sequences taking on $N$ values.
Of course, nitpickers might argue that y takes on $N+1$ values, not $N$
as required by the OP, but in the case $N = 3$, there is an easy workaround.
Define
$$z[n] = \begin{cases}1, & n = 0,\\
-e^{j2\pi/3}, & n = 1,\\0, &\text{otherwise,}\end{cases}$$
which is clearly a three-valued sequence, and note that
$$\displaystyle \sum_{n=-\infty}^\infty x[k−n]z[n]=
\sum_{n=0}^1 e^{j2\pi (k-n)/3}z(n) = e^{j2\pi k/3}\cdot 1
+ e^{j2\pi (k-1)/3}\cdot (-e^{j2\pi/3})
= e^{j2\pi k/3} - e^{j2\pi k/3} = 0.$$ The real pezzonovante nitpickers
who wish to complain that the three values that $z$ takes on are
not the same three values that $x$ takes on are invited to work out
a more satisfactory solution for themselves.

It is tempting to extend $y$ periodically
to be nonzero for all $n$ so as to totally avoid the nitpick, but unfortunately,
the infinite sum $\sum_{n=-\infty}^\infty e^{j2\pi n/3}$ that is the
result of the convolution is not absolutely
convergent.  Of course, since both sequences are of period $3$, their
 convolution could also be defined as a periodic convolution,
which is the sum of just $3$ terms, and we are back in business.
More generally, note that the periodic
convolution of any two rows of the $N\times N$ DFT matrix is $0$,
and so we can extend this notion to sequences that take on $N$
distinct values for any $N \geq 3$.  For another viewpoint on this
approach, see the answer by @Hilmar.


Why Convolution between Signal and a Channel Is Not Done Directly in SIMO Channel In OFDM system, I need to transmit a signal $x$ in SIMO channel $H$ with one $Tx$ antenna and 4 $Rx$ receiver's antennas. The initial equation for that convolution is $r=Hx+n$ where $n$ is the noise. $r$ is the received signal, $x$ is the transmitted signal, and $H$ is a toeplitz matrix as built in below code.

The question, When I created the toeplitz matrix $H$, the dimensions of  transmitted signal $x$ and toeplitz matrix $H$ are not equal,  so we are unable to do the multiplication!. I checked online, I got, before processing that multiplication, I should do two steps. first get $x1$ and then $x2$ as in below code.  the $x2$ is the flip of $x1$, is theoretically known, but what's the benefit of getting the $x1$ ?

EDIT:
Regaring the dimension of toeplitz matrix $H$, is as below:
if we have the length of signal $Q$, the length of channel $M$, the number of received antenna $p$. So, the dimension of toeplitz matrix will be (Q*p x Q+M)

thnx

I've said, in another comment, convolution using the function conv "i.e in MATLAB" and convolution using the Toeplitx matrix must give the same results. That's ok.

Now, according to your code, the received signal $r$ is the results of convolution between channel $h$ and emitted signal $x$, which means $r = h*x + n$
* indicate to the convolution (which is circular convolution in case of using $CP$ with OFDM or with any other system).

So, in that conventional known case, we are using SISO system, where 1 antenna is used as transmitter and 1 as receiver, The length of our parameters should be:

$x = N$ x 1 ; $h = L$ x 1;

$N$ is the length of our signal, in your case $Q$

$L$ the length of channel or we sometimes call it IR, in your case $M$

Till here, it's clear, it's the normal process which can be read anywhere. Now suppose that you are using $SIMO$ system with $P$ antennas at reciever instead of 1 (by the way fractional sampling is also equivalent to SIMO system).  in that case you suppose to have $P$ copies for your signal instead of one. It's like you are doing $P$ times convolutions for your emitted signal with $P$ different channels. let's say $P$ = 4; means you have 4 receiver's antennas equivalent in our case into 4 different channels.

As mentioned, you suppose to have 4 copies for your signal, and 4 different channels compared with conventional case, so let's say parameter $H$ is toeplitz matrix which will represent those four channels in SIMO system. So we will have

$R = HX + N$

$R$ is received signal in SIMO, $H$ is the toeplitx matrix, $X$ emitted signal in SIMO system too and also $N$ represent the noise.

Now, your question how to build $H$ and what's dimension of $H$ when using SIMO system. (by they way, your code is correct)

The emitted signal $X$, represents 4 copies of $x$, so its dimension is $(N +L)$ x 1;

The toeplitz matrix $H$ has a dimension $NP$ x $(N+L)$ where $P$ = 4 in your case.  And the noise $N$ should be $NP$ x 1. (N as noise is different about N of emitted signal in this case).

So, Now, you can conclude the dimension of received signal $r$ and $R$ in every case easily.

and regarding your question, why using $x1$, it's to have dimension of emitted signal $X$ of $(N +L)$ x 1 instead of dimension $N$ in conventional case

thnx


Understanding the Parameters of the Bilateral Filter Allow me to show you a snip from Szeliski's Computer Vision book on bilateral filter. I don't know how to interpret these kernel definitions:

In the actual summation what values do the $k$ and $l$ take? E.g For a kernel of size 5x5 does $k$ take values between negative and positive inf? or between 1 to 5 (e.g. in matlab programming)? Or does it take values from -2 to 2?

Also what do those double bars mean in the definition of range kernel?

A small discussion of the details you've asked for:

The bilateral ﬁlter is controlled by two parameters: $\sigma$d and $\sigma$r.

• As the range parameter $\sigma$r increases, the bilateral ﬁlter becomes closer to Gaussian blur because the range Gaussian is ﬂatter i.e., almost a constant over the intensity interval covered by the image.

• Increasing the spatial parameter $\sigma$d smooths larger features.

An important characteristic of bilateral ﬁltering is that the weights are multiplied, which
implies that as soon as one of the weight is close to 0, no smoothing occurs. As an example, a large spatial Gaussian coupled with narrow range Gaussian achieves a limited smoothing although the ﬁlter has large spatial extent. The range weight enforces a strict preservation of the contours.

"The Norm or the double bars indicate the gaussian distance in the equation.This distance is deﬁned by Gσ(||p − q||), where σ is a parameter deﬁning the extension of the neighborhood."

Since (i,j) and (k,l) are simply spatial points on an image, they will vary from the start of the image to the end of the image. Typically this would be dependent on your indexing measure, for example in a matrix, you've have to range from 0 to 5 for a 5x5 image. The actual values of k,l matter little. Its the relative value which matters far more.


Normalizing a spectrogram or a pitch class profile I am developing a chord recognition program using neural networks.I need to normalize the spectrogram or my pitch class profile to provide inputs for the neural network. The pitch class profile classifies the frequency into 12 class semitones (From note A to G).

The problem here is it produces really high values and low values: they can be on the order of ten, one million, or even a trillion. Meanwhile, the input values for the neural network should be between 1 and -1.

What is a normalization formula that I can use to transform the values to be between 1 and -1?

Audio volume is usually measured in decibels (a logarithmic scale) precisely because there's such a wide range of input energies. But in your case, it's unlikely that you need volume at all. What you need is relative volume: what fraction of the total energy in the signal is in each FFT bin? Since the energy in each FFT bin is positive, it follows that the relative contribution of each FFT bin is strictly in the range [0,1]. Multiply by 2, subtract 1, and you're in the correct range. And as a bonus, the extremes of the range are quite unlikely in practice.


Head detection using HOG and SVM I am currently trying to detect heads in a sequence of real-footage images and am using HOG feature descriptor and SVM as classifier.

Currently I am using Dalal's HOG implementation code in MATLAB found in this link:
http://www.mathworks.com/matlabcentral/fileexchange/46408-histogram-of-oriented-gradients--hog--code-using-matlab

Currently I am using libSVM MATLAB version found in this link:
http://www.csie.ntu.edu.tw/~cjlin/libsvm/

I prepared 350 positive training images and 1243 negative training images.

The hog feature vectors extracted from the training images are converted to libsvmFormat and inputted in the libsvm training method, to obtain a model. The hog vector length is that of 1764.

Regarding the libSVM I chose these as parameters:

-s 0  (i.e. C-SVC)
-c 1  (i.e. cost = 1)
-t 2  (i.e. kernel = RBF)
-g 3  (i.e. gamma = 3 (this is for kernel))

Regarding the HOG version, I left the cell, bin & block parameters as they were in the implementation shown in the link above.

I am using a scanning window of size 128x128 and 256x256 to scan through the whole image to detect possible heads.  At each window, the hog feature vector is extracted for each image and is inputted in libsvm predict, to test whether it should be classified as a head or not.

However, after doing all the above, I have a numerous amount of false negatives and can't figure out what I am doing wrong.

Can someone experience please offer some advice on what is possibly wrong?  I really need to figure this out please.  Much appreciated!

The thing to keep in mind here is that HOG is not invariant to in-plane rotation. A change in orientation of more than 10-15 degrees will probably throw it off.  So if your head can have different orientations in the image, you would either need to train multiple detectors or use something other than HOG.

By the way, there is a function extractHOGFeatures in the Computer Vision System Toolbox, and SVM is available in the Statistics Toolbox. See this example.

Edit: from your sample images it looks like in most cases you need to detect profile faces. vision.CascadeObjectDetector in the Computer Vision System Toolbox can do that for you. For greater robustness, I would run two detectors in succession: one for profile faces and one for frontal faces.  And also, HOG may not be your best choice here. For face detection Haar or LBP features generally work better.

No matter what you do, you will have some false negatives.  So, if you are using a video camera, I would also recommend tracking the faces.  See this example. This will give you an opportunity to correct mistakes made on individual frames.

Edit2: You can also train your own detector with trainCascadeObjectDetector, and you can use trainingImageLabeler app to label the positive samples.


determine the impulse response of given second order difference equation $y[n]-3y[n-1]-4y[n-2]=x[n]$

i solved it as follows :

1.first calculated natural-response(ZIR i.e,zero input response)

by using characteristic equation it came out as

  $y_n[n]=[c_1(-1)^{n}+c_2(4)^{n}]$ here i couldn't be able to calculate $c_1$ and $c_2$ because there were no initial conditions given in the question.

2.And,Then calculated forced-response(ZSR i.e,with all initial conditions zero) and ,since we know it consists of two parts

homogeneous solution and particular solution i.e,($ y_h[n] +y_p[n]$)
i calculated both of parts separately as follows:
first $ y_p[n]=0 $ because input is $\delta [n] $
and secondly $y_h[n] $ =$[C_1(-1)^{n}+C_2(4)^{n}$ (here i was able to calculate values of $C_1$ and $C_2$ as $\dfrac{-1}{5}$ and $\dfrac{6}{5}$ respectively by using zero initial conditions i.e, $y[-1]=y[-2]=0$,$x(0)=1$ ,as per definition of forced response)

3.and finally i added all solutions to get total solution i.e, $y[n]= y_n[n]+y_h[n]+y_p[n]$ and got answer as $ y_n[n]=[c_1(-1)^{n}+c_2(4)^{n}-0.2(-1)^{n}+1.2(4)^{n}+0] $
but answer in the book is given as $y [n]=[-0.2(-1)^{n}+1.2(4)^{n}]$ so my problem is why didn't they incorporated natural response term in final answer  .
thus, my questions are

  1.does "finding impulse response in numerical" means finding only "forced response" and leaving natural response.terms.unlike we do in other input cases(like $x_\left[n\right]=\left[\dfrac{1}{2}\right]^n$) where we also calculate natural responses(where of course initial conditions are provided so that we can calculate $c_1$ and $c_2$ as well)...... is impulse input an exceptional case of inputs?

2.same they have done for questions regarding 'step response'

i haven't read z-transform and generating function method to solve difference equations so please do not give answer by using those techniques.

The impulse response of a system is its zero-state response to an impulse at the input. If the system is linear and time-invariant (LTI), then the system's response to any input signal can be described in terms of the impulse response. The impulse response only depends on the structure of the system, and it cannot depend on any initial conditions. In fact, a system with non-zero initial conditions is not linear because a part of the output (the zero-input response) does not depend on the input signal but solely on the initial conditions and the system itself.


Create an aliased image in FD Is there an accurate way to create an aliased image from the Fourier transform of the original image?

in other words, i have the Fourier coefficients of an image and i want to make down-sample in frequency domain.

I am trying to construct an aliasing image.

I=rgb2gray(imread('Liver.png')); % a 512*512 image
subplot(1,2,1),imshow(I),title('Original image')  % shown on the left side of figure below

when you implement F=fftshift(fft2(I));, the corresponding axis represent the frequency range of [-255:256]/512; and the cut-off frequency U = V = 256/512 = 0.5;

based on Nyquist Sampling Theorem, the sampling periods in spatial domain should satisfy

delta_x<=1/(2*U)=1;
delta_y<=1/(2*V)=1;

to avoid the aliasing.

So if your bandwidth of recorder is 64 * 64, the corresponding point interval in spatial domain is 1 / (64/256) = 4. Thus,

Ia=imresize(I(1:4:end,1:4:end),[512 512]);
subplot(1,2,2),imshow(Ia),title('Aliased image') % shown on the right side of the figure below

You can observe the aliasing effect at the edge.


Interpretation of sample rate I'm a little confused about the sample rate:

If I have an audio file that was recorded with a sample rate of 48kHz, this means that there are 48,000 values per second in this file, right?

So when I read this file with a smaller sample rate of lets say 44.1 kHz I am summarizing every 44,100 values to one second so effectively I make the audio slower and thus longer than it actually is. Am I getting this correctly or do I have a total misunderstanding?

Yes that is correct-- if you sampled the audio recording at 48 KHz and then played those same samples back at 44.1 KHz you would be "time-stretching" the recording and it would be slower. However you can also resample the samples by the ratio of 147/160 to then have the 44.1KHz samples play back at the same speed.