What DSP method are consumer products like the game Rocksmith and guitar polyphonic tuners likely to use for multiple f0 detection? The game Rocksmith and some polyphonic guitar tuners (like the PolyTune) are able to recognize, at least to some extent, multiple fundamental frequencies in single-source audio.

I know there are some methods for multiple fundamental frequency detection, like MUSIC and ESPRIT, but is it likely that these products implement them or leverage several assumptions (such as: audio is from a guitar, audio is from a single source) about the application domain to make a simpler method work?

In the second case, which simpler methods are likely to be used? Simple peak-picking on a spectrum, with domain-specific heuristics? Or something more complex?

Check out chapter 1.3 of this IRCAM paper on multi-F0 estimation. It discusses the difficulties in extracting multiple F0s from a recording, including the handling of overlapping partials, transients, and reverberation, as well as the modeling of domain-specific sources with varied spectral properties.


Pole zero plot, normalizing frequency response plot? I'm asked to plot the frequency response (amplitude) given a specific pole-zero diagram.

$$ H(z) = H_0 \frac{\prod\limits_{m=1}^{M} (z - q_m)}{\prod\limits_{m=1}^{M} (z - p_m)}$$

$$ H(e^{i\omega}) = H_0 \frac{\prod\limits_{m=1}^{M}(e^{i\omega} - q_m)}{\prod\limits_{m=1}^{M}(e^{i\omega} - p_m)}$$

If I understood it correctly, the amplitude at frequency $\omega$ is (the magnitude of the distance from $e^{i\omega}$ to all the zeroes) divided by (the magnitude of the distance from $e^{i\omega}$ to all the poles), i.e:

$$ \Big|H(e^{i\omega})\Big| = \Big|H_0\Big| \frac{\prod\limits_{m=1}^{M}|e^{i\omega} - q_m|}{\prod\limits_{m=1}^{M}|e^{i\omega} - p_m|}$$
where $q_m$ are the zeroes and $p_m$ are the poles and $H_0$ is the constant gain factor.

The problem is that in the graph where I need to draw the frequency response, the frequency and amplitude range $0\to1$ like so:

After I get a value by calculating the poles and zeroes I almost always get a value above $1$. What do I need to do with the value to fit it into the graph? How do I normalize(?) the value?

For the frequency axis, you're going to have to normalize by $\pi$ to plot digital frequencies in $[0,\pi]$.

For the magnitude axis, normalize by $\mathrm{max}\left(|H\left(e^{i\omega}\right)|\right)$.

I'm surprised you weren't asked to plot magnitude response in dB (a logarithmic axis).


bit reversed order of FFT matrix for channel estimation Although the question seems similar to this one HERE, but  what I want to ask about is little bit different.
Assume we have a signal $X$ whose length is $N$x$1$, we convert that signal into time domain using bit reversed order FFT as: $x = F^{'b}X = p*F'*X$ , where $F^{'b}$ is the inversed bit-reversed order FFT matrix ($F^{b}$). $p$ is row permucation matrix and $F'$ is the inversed DFT matrix.
A cyclic prefix is added into the signal $x$ having $x_{cp} = [cp; x]$. Then that signal is transmitted over a frequency selective channel $h$ resulting $y_{cp} = h*x_{cp}$  where $*$ denotes to the circular convolution operation.
The CP is first removed from the signal $y_{cp}$ resulting the signal $y$.
Based on my analysis, when transforming the signal $y$ into frequency domain using $F^{b}$, as following $Y = F^{b}y = p'*F*y$, the channel is diagonalized where I can get the channel $h = F^{'b}(Y./X)$ it means in that case I use the whole signal $X$ as pilots.
The problem I am facing, when using such pilots from $X$, let's be $X_{1:4:end}$,  and then perform the same steps, the results of $h = F_2^{'b}(Y_{1:4:end}  ./ X_{1:4:end})$ is not equivalent to the channel !! where $F_2^{'b}$ is similar to $F^{'b}$ but with size similar to $(Y_{1:4:end}/X_{1:4:end})$ and $./$ is element wise division.
$NP$: When performing these steps with the normal $FFT$ matrix, everything is working well and I can recover the exact channel using the pilots data.

So, writing your $F^b$ and $F^{'b}$ as permutations $p, p^{-1}$ of the DFT matrix $F$ , IDFT $F'$, and the cyclic prefixing and removal also as Matrix operation $P$, $R$ (which boils down to cyclically extended or truncated identity matrices), and the circular convolution with the channel impulse response as circulant matrix $C$,
\begin{align}
x&=F^{'b}X=p F' X\\
x_{cp} &= PpF' X\\
y_{cp} &= y*x_{cp} \\
y&=R(y*x_{cp})= CpF'X\\
Y&=F^By \\
&= p'FCpF'X & F, F' \text{ are symmetric matrices!}\\
&=p'FCF'pX\\
&=p'\underbrace{FCF'}_{\text{IDFT-Channel-DFT}}\underbrace{pX}_{\text{permuted symbols}}
\end{align}
So, all that's really happening here mathematically is that you permute the symbols before transmission, then de-permute them at the receiver. In between, you get normal OFDM. There's no advantage or disadvantage to the permutation w.r.t. to channel diagonalizability, but we've arrived here in your previous question already.
Now, for channel estimation, you'll want to insert pilots such that they make sense for the central $FCF'$ part – in other words, you need to think about where you'd want your normal OFDM system to have these pilots, and then shift these positions according to your $p$.
Of course, the channel you estimate left of your final reordering isn't the same as the channel you estimate with OFDM, it's a permutation. As shown above, it's however "just" a permutation, there's nothing in the maths here that changes anything about the channels properties at all, so if your results aren't the same as for straightforward OFDM, then you've got an implementation bug somewhere.


How does it go about implementing pitch shifting without affecting the duration or speed on matlab? I am creating a project that will record a sound sample from the world, for example, a bird's chirping. And then through pitch scaling, I could create different notes.

I first tried to change the sampling rate fs below in the Matlab code, and I noticed that it will change both the duration and the pitch. After some research, i found out this method is called sample rate conversion. but I still do not quite understand why with a higher sampling rate, the octave created will sound slower and lower, and vice versa.

It would be really helpful if someone could provide some insight into this. Thank you in advance!

function [x1, x2, x3, x4, x5, x6, x7] = signal_generator(n1,n2,n3,n4,n5,n6,n7,fs)
% n is the signal frequency, fs is the sampling frequency

m = fs; % how many samples per second
t = 0:(m-1);
x1 = sin(2*pi*n1/fs*t);
x2 = sin(2*pi*n2/fs*t);
x3 = sin(2*pi*n3/fs*t);
x4 = sin(2*pi*n4/fs*t);
x5 = sin(2*pi*n5/fs*t);
x6 = sin(2*pi*n6/fs*t);
x7 = sin(2*pi*n7/fs*t);

end

For the speed change: Sample rate is the rate at which samples are read right? So, as per that definition, since you have a higher sample rate, you would be reading more samples in one second. The length of the file remains the same but now since the samples are read faster, the overall length of the file would become slower (with an increase in sampling frequency).

For the pitch change: Now, let us consider the case of increasing the sampling frequency. I'm sorry for giving a rough explanation for this but imagine increasing the sampling frequency would squish the signal in the time axis since now lesser samples represent the same length in samples. This would mean an increase in frequency, thus resulting in a pitch shift.


ARMA Filter Output Stationary and set up? I have questions regarding ARMA Filters.
Is the output of a ARMA Filter stationary or just wide sense stationary?
I do know that you can obtain an ARMA filter by connecting an MA filter with an AR filter (MA,AR). Does that work the other way around? Can i obtain an ARMA filter by switching the order and connecting an AR filter with an MA filter? (AR,MA)

Hi: The use of the term "stationary" is vague because there are a few different kinds. But usually, the term "stationary" is meant to mean wide-sense stationary. But, since you've also used the term wide sense stationary in your question, I'm thinking that you're use of the term stationary might just be referring to the output ( by output, I'm assuming you mean the estimated residual noise term )  having a constant (zero) mean and constant variance.
Wide sense stationary refers to constant mean, constant variance and $corr(X_{t1}, X_{t2})$ being only a function of, the time difference, $t2 - t1$.
But, the general process in ARMA modelling ( in statistical time series that is. dsp could be different), is to take the estimated residuals and check that they have a  constant mean of zero, constant variance and that they are uncorrelated. This is because the ARMA framework assumes that the error term is white noise so the estimated residuals should be uncorrelated and have zero mean and constant variance.
So, the test that is done after the ARMA model is identified and estimated
is testing that the estimated residuals are white noise. But, white noise is a specific case of stationarity so one could argue that you're kind of testing for wide-sense stationarity ( except that corr should be zero so not even a function of the time difference ).
Note that there is another term called "strict sense stationarity" which is a stronger form of stationarity where one requires that the joint distribution of  $X_{1}, \ldots, X_{n}$ is the same as the joint distribution of  $X_{t+1}, \ldots, X_{t+n}$.


Is discrete cosine transform differentiable? Is DCT differentiable? I.e. is $\frac{\mathrm{d}\mathrm{DCT}(x)}{\mathrm{d} x}$ defined?

I've seen a few implementations in automatic differentiation software packages, so I suppose the answer is positive, but I'm still missing a mathematical proof.

The Discrete Cosine Transform of a discrete sequence $x$ can be defined as:

$$
{\tt DCT(}x{\tt )} = X_k = \frac{1}{2} (x_0 + (-1)^k x_{N-1})
 + \sum_{n=1}^{N-2} x_n \cos \left[\frac{\pi}{N-1} n k \right] \quad \quad k = 0, \dots, N-1.
$$

If, by differentiable, you mean:
$$
\frac{\partial {\tt DCT(}x{\tt )}}{\partial x}
$$
then isn't it just
$$
\frac{\partial {\tt DCT(}x{\tt )}}{\partial x} = \frac{\partial X_k}{\partial x} =\left [
\begin{array}{ccccc}
\frac{1}{2} & \cdots & \cos \left[\frac{\pi}{N-1} n k \right] & \cdots & (-1)^k\frac{1}{2}
\end{array}
\right ]
$$
where $x = [ x_0, x_1, \ldots, x_{N-2}, x_{N-1} ]^T$ ?

What do you mean by a "mathematical proof"?


Can you count signal energy equivalent in frequency domain? I read about method of analyzing audio signals where FFT result is bandpass-filtered and the energy of this band is counted.
Can you count the momentary energy in frequency domain?
Does this mean summing magnitude value from all frequency bins (after filtering)?

(I don't care whether this will yield the value similar to the one counted in time domain but just if this can be of some benefit - reflect the energy present within this band.)

Parseval's theorem applied to the signal itself:

$$\sum_{n=0}^{N-1}|x[n]|^2 = \frac{1}{N}\sum_{k=0}^{N-1}|X[k]|^2$$

So, if you sum the magnitude squared of the frequency coefficients of the DFT, you get the energy of your time-domain vector.

If you want power (i.e. energy per sample), divide both sides by $N$.


Why is dirac delta used in continuous signal sampling? Why this concept is the most widely accepted model of signal sampling? By multiplying the continuous signal value with the dirac delta we get an infinite value. However if we perform convolution of our signal and the dirac comb, we get our signal again - how is it useful?

The answer to a similar question has been given here, but it's still not 100% clear to me.

Sampling of a continuous function: Kronecker's or Dirac's delta?

The best explanation of this that I've seen is from the "Digital Signal Processing Handbook" by Madisetti. Essentially the multiplication by the delta function is equivalent to sampling because the Fourier transforms are the same. So although the result of $s(t)\delta(t-nT)$ may not make much sense it's Fourier transform does exist because we are taking the integral over the delta function.

I'm cutting and pasting his text and redoing the equations. Note CT=continuous time, DT= discrete time. I have not included figures 1.1 and 1.2.

The relationship between the CT and the DT domains is characterized by the operations of sampling and reconstruction. If $s_a(t)$ denotes a signal $s(t)$ that has been uniformly sampled every T seconds, then the mathematical representation of $s_a(t)$ is given by

$s_a(t) =\sum_{n=- \infty}^{\infty} s(t)\delta(t-nT) \qquad$   (1.1)

where $\delta(t)$ is a CT impulse function defined to be zero for all $t\not=0$, undefined at t=0, and has
unit area when integrated fromt $t=-\infty$ to $t=\infty$. Because the only places at which the product $s(t)\delta(t−nT)$ is not identically equal to zero are at the sampling instances, $s(t)$ in (1.1) can be replaced with $s(nT)$ without changing the overall meaning of the expression. Hence, an alternate expression for $s_a(t)$ that is often useful in Fourier analysis is given by

$s_a(t) =\sum_{n=- \infty}^{\infty} s(nT)\delta(t-nT) \qquad$       (1.2)

The CT sampling model $s_a(t)$ consists of a sequence of CT impulse functions uniformly spaced at intervals of $T$ seconds and weighted by the values of the signal $s(t)$ at the sampling instants, as depicted in Fig.1.1. Note that $s_a(t)$ is not defined at the sampling instants because the CT impulse function itself is not defined at $t$. However, the values of $s(t)$ at the sampling instants are imbedded as “area under the curve” of $s_a(t)$, and as such represent a useful mathematical model of the sampling process. In the DT domain the sampling model is simply the sequence defined by taking the values
of $s(t)$  at the sampling instants, i.e.,

$s[n]=s(t)|_{t=nT}  \qquad $(1.3)

In contrast to $s_a(t)$, which is not defined at the sampling instants, $s[n]$ is well defined at the sampling instants, as illustrated in Fig.1.2. Thus, it is now clear that $s_a(t)$ and $s[n]$ are different but equivalent models of the sampling process in the CT and DT domains, respectively. They are both useful for signal analysis in their corresponding domains. Their equivalence is established by the fact that they
have equal spectra in the Fourier domain, and that the underlying CT signal from which $s_a(t)$ and $s[n]$ are derived can be recovered from either sampling representation, provided a sufficiently large sampling rate is used in the sampling operation (see below).

Cheers,

David


Effect of zero padding an odd symmetric FIR filter in the time domain I have a symmetric lowpass FIR filter with 1149 time domain taps (all real  coefficients). For implementation purposes, it would be easier if the filter had 1200 taps.

Since it has an odd number taps, can I zero pad it to 1200 time taps and keep the phase and magnitude response the same?

It seems like these are my options:

[25 zeros, 1149 taps, 26 zeros] which delays the filter but doesn't change the frequency response
[51 zeros, 1149 taps]
[1149 taps, 51 zeros]

The three options you mention all result in the same magnitude response. The ones with zeros before (to the left) of the actual impulse response will just add delay (as many samples as there are zeros). The last option has exactly the same frequency response as the original filter, without any extra delay.

I think it's pointless to add extra delay by prepending zeros. Zero padding at the end is OK, if you have a good reason for it.
If you really want a filter with $1200$ taps, why wouldn't you use the extra taps to get a better filter?

If you need a linear phase filter with an integer delay (in samples) then you could design a filter with $1199$ taps and append one zero tap at the end (if for whatever reason you need exactly $1200$ taps). This is because the delay in samples of a length $N$ linear phase filter equals $(N-1)/2$, which is only integer if $N$ is odd.


Why can't realisable IIR filters have linear phase? I am studying IIR Filter Design and came across this arbitrary statement in my textbook which says that 'physically realisable and stable IIR filters can not have linear phase'. Would really appreciate if anyone can elaborate on this.

The frequency response of a real-valued discrete-time system with linear phase has the form

$$H(e^{j\omega})=A(\omega)e^{-j\omega\tau},\qquad\omega\in [-\pi,\pi]\tag{1}$$

where $A(\omega)$ is either a real-valued even function or a purely imaginary odd function, and $\tau$ is some real-valued parameter (the delay). If $A(\omega)$ is purely imaginary, then there is a phase jump of $\pi$ at $\omega=0$, but otherwise the phase is linear. This latter type of phase is usually referred to as "generalized linear phase", and it is the type of phase exhibited by filters with odd symmetry in the time domain.

The impulse response of the corresponding system is given by

$$h[n]=(a\star d)[n]\tag{2}$$

where $a[n]$ is the inverse discrete-time Fourier transform (IDTFT) of $A(\omega)$, $d[n]$ is the IDTFT of $e^{-j\omega\tau}$, and $\star$ denotes convolution.

The sequence $a[n]$ is either symmetrical or anti-symmetrical due to the properties of $A(\omega)$ (real or purely imaginary), and hence non-causal. The sequence $d[n]$ is given by

$$d[n]=\frac{\sin[\pi(n-\tau)]}{\pi(n-\tau)}\tag{3}$$

which simplifies to $\delta[n-\tau]$ if $\tau$ is an integer.

Consequently, if $a[n]$ has infinite length (as is the case for IIR filters), $h[n]$ cannot be causal and have linear phase at the same time. Only if $a[n]$ has finite length (thus FIR) can it be made causal by shifting to the right. That's why FIR filters can have linear phase whereas causal IIR filters cannot.


Suggestion for a DSP board to implement a Kalman filter in real-time at a rate of around 4MHz or higher I have been attempting to implement a Kalman filter on an FPGA in VHDL and have got to the point that it works. However I have the issue that each variable requires 64 bits to be represented to get the expected behaviour in simulation. With my current implementation about 95% usage of the FPGA is required for a 2D Kalman filter, when I really require a minimum of a 3D filter for my application. It also only runs at a rate of ~400KHz when I require 4MHz (ideally more) for my application.

I require an operating rate of 4MHz, i.e. it reads from the ADC, performs an iteration of the Kalman filter (a 3D Kalman filter) and writes to the DAC in 250ns. Ideally I would have a 16 bit ADC and DAC (could maybe be slightly lower) between $-1V$ to $+1V$.

A colleague proposed I use a DSP instead as the floating point implementation the programming a lot easier. In addition, as I have much more experience in writing efficient C code I am more inclined to use a DSP to solve my problem, in fact I already have a working C implementation that I prototyped before implemented the VHDL implementation. The C code I completed in a couple of days, whereas the VHDL implementation has so far taken me a couple of months of work (albeit not continuous work).

In my research I came across the TMS320C6678 Lite Evaluation Module which, with it's high GFlop rate of around 20GFlops (with a single core, 160GFlops if I can make use of all 8), should allow me to have ~5000 Flops with which to implement a single Kalman filter operation, which should be more than enough. Another StackOverFlow user, however informs me that connecting an ADC/DAC with 16 bit precision and a sample rate of around 4MS/s would be very difficult due to having to use a PCIexpress bus, is this correct? I saw some people online discussing using the SPI connection to connect the specified board to an ADC/DAC.

(100MS/s is what I can currently achieve with a NI 5781 ADC/DAC and a NI PXIe-7961R FPGA, although my current Kalman filter implementation cannot operate at these sorts of speeds) Is there a DSP which is more appropriate for my application and will allow me to achieve the sorts of speed I require and be able to interface with an ADC/DAC at the sorts of speeds I require? (Using either fixed point or floating point arithmetic)

If this is not the case, would it perhaps be worth learning to use Vivado High-Level Synthesis (HLS) so that my C code can be turned into VHDL, would this perhaps be more efficient than I can achieve with my current knowledge of/skill in VHDL?

Are there any other device options I have missed? I remember someone suggesting a software define radio at some point for my application?

If none of these is a reasonable option, could someone suggest some online/book resources for learning techniques to make my VHDL code more efficient?

The covariance matrix is a 2x2 marrix, each value is stored using 64 bits. My clock is a 40MHz clock downsampled by 9 to be 3.0769 (as my synthesiser reports it must be about 3MHz or less). Observations are 16 bit samples at 40/13/7 MHz ~400KHz (as each kalman filter iteration takes 7 clock cycles). Also my current implementation is all fixed point, done with the signed data type. For most values I multiply all the values by 2*N and truncate them to integers before storing them as the signed data type.

And for the elements of the covariance matrix, which are smaller, I use a different value, 2**M offset. I once tried to change it so they each had their own fixed point location/multiplicative offset but each time I operate 2 numbers together I need to bit-shift them so they have the same offset/fixed point location otherwise I couldn't correct the answer back at the end by bit-shifting. This lead to increased resource usage and therefore I stopped attempting this.

The conclusion of my discussion with hops and Marcus was that getting a DSP to operate with an ADC and DAC at the latency required for my application (<=250ns for 4MHz) would be difficult and require writing some custom interface. We concluded that the best option is to continue to work with the FPGA and improve my VHDL implementation using pipelining, bit widening to reduce my usage of large numbers of bits and turning my functions into architecture blocks. An alternative solution was to attempt to use Vivado HLS to generate VHDL code from C code.


Does a collection of Gaussian random variables necessarily constitute a Gaussian Process? If $\{X(t)\}$ is a Gaussian Process then the random variables $X(t_k)$ where $k = 1,2,3...n$, are jointly Gaussian.

If each random variable $X(t)$ is a Gaussian variable, then will the random process $\{X(t)\}$ be a Gaussian Process?

I am asking because while studying the Statistical proof of why we model noise in communication as Gaussian Process $N(t)$, I came across the fact that noise $N(t)$ at time $t$ is basically a sum of micro phenomenon which makes the random variable $N(t)$ a Gaussian random variable. Now that random variable $N(t)$ is Gaussian how do we prove that the random process formed by these random variables is a Gaussian Process?

As you correctly say, if $\{X(t)\colon t \in \mathbb T\}$ is a Gaussian process, then for every possible choice of $n$ random variables $X(t_1), X(t_2), \ldots, X(t_n),$ where $t_k \in \mathbb T$ for $ k = 1,2,\ldots, n$, the $n$ variables are jointly Gaussian. Indeed, that is the definition of a Gaussian random process: a collection of random variables such that every finite subset of random variables is jointly Gaussian, that is, has a multivariate normal density.  Note that the individual random variables are also themselves Gaussian random variables.  Note also that the definition does not require that the index set $\mathbb T$ be a countable set such as $\mathbb Z$, the set to all integers (which would make the process a discrete-time process); having $\mathbb T = \mathbb R$, the set of real numbers, is fine (this makes the process a continuous-time process): but in either case,the requirement is still the same that every finite subset of random variables that you can choose from the process has a jointly Gaussian distribution. There is no such thing as an infinite-dimensional joint distribution function (joint Gaussian or otherwise), and so the definition is silent about what can be said about the joint distribution of infinite subsets of the random variables.

You are asking whether it is true that if we have

  a collection of Gaussian random variables, specifically the set $\{X(t)\colon t \in \mathbb T\}$, such that each member $X(t)$ of the set is a Gaussian random variable,

then

  is it always the case that the set constitutes a Gaussian process?

The answer is No if you accept the definition above because Gaussian random variables don't have to be jointly Gaussian, and is so there is no guarantee that the chosen collection of Gaussian random variables is also a collection of jointly Gaussian random variables, as required in the definition of Gaussian random process.  An incredibly diverse collection of joint densities of two Gaussian random variables that are not jointly Gaussian can be found in this answer on stats.SE.

Turning to your "micro phenomenon" comment, it seems that you are (or the author of your textbook is) thinking that the assumption of individual Gaussianity is justified from a central limit theorem type view of the distribution of each individual random variable.  This is not quite correct. The Gaussian random variable model usually comes from an assumption that one  is dealing with thermal noise which is caused by the random motion of electrons inside a conductor. The (small) charge on each electron induces a (small) voltage across the terminals of the conductor depending on the position of the electron with respect to the terminals, and all the zillions of such voltages add up (or cancel out) to produce a net voltage  which can be modeled a Gaussian random variable. This net voltage varies with time as the electrons jostle about, creating a random process that is generally modeled as a Gaussian process, that is, we make a leap of faith that the random variables are not just individually Gaussian (which we can justify as above in a hand-waving sort of way) but also jointly Gaussian. One justification for this leap of faith is that if we imagine a linear system with white Gaussian noise as input, then the output process is a Gaussian process. More important, the linear system with white Gaussian noise input model has considerable support from experimental evidence it matches up well with physical reality.

@AlexTP Uncorrelated Gaussian random variables are not necessarily jointly Gaussian random variables. The classic example is $X \sim N(0,1), Y = ZX$ where $Z$ takes on values $\pm 1$ with equal probability $\frac 12$. Then, $Y \sim N(0,1)$ also. Furthermore, $X$ and $Y$ are uncorrelated random variables, but $X$ and $Y$ are not jointly Gaussian random variables. In fact, all the probability mass lies along two straight lines intersecting at right angles instead of being distributed with circular symmetry in the plane as happens with jointly Gaussian random variables with equal variances.


Why is the range of frequency for discrete time Fourier transform $-\pi<\omega<\pi$? In my class we are taught that the range for the frequency is $-\pi<\omega<\pi$ for discrete time Fourier transform, however for continuous time the limit is $-\infty<\omega<\infty$
why is this the case?

Short intuitive explaination: In a time discrete system, the spectrum repeats itself at multiples of the sampling frequency (normalized sampling frequency is $2\pi$). Because of that, the limits of the fourier integral can be reduced from $\pm \infty$ to $\pm \pi$.


Gaussian noise generation for a given SNR ? I am trying to add a Gaussian noise, normal distributed to a signal I have simulated (sig_noiseFree), to get a noisy signal (sig_noisy). lease see the code.

I wrote the function my self according to what I understood. Can please anyone confirm that I am doing this right ? especially the noise addition and SNR calculation ?

I have a second question : I understand that the variable "k" controls the noise level, how can I set the SNR and generate the suitable k so I get the noisy signal with the SNR wanted ?
Thank you in advance.

def SNR(sig, noise,dt):
    Signal = np.sum(np.abs(np.fft.fft(sig)*dt)**2)/len(np.fft.fft(sig))
    Noise = np.sum(np.abs(np.fft.fft(noise)*dt)**2)/len(np.fft.fft(noise))
    return (10 * np.log10(Signal/Noise))

def Gauss_noise(k,sig):
     return np.random.normal(scale=k*np.max(sig), size=len(sig))

def noisy_sig(sig,k, dt):
    return np.fft.ifft(np.fft.fft(sig)*dt+np.fft.fft(np.random.normal(scale= k*np.max(sig), size=len(sig)))*dt)/dt

def main(argv):
    k = 2e-2
    noise_g = Gauss_noise(k,sig_noiseFree)
    sig_noisy = noisy_sig(sig_noiseFree,k, dt)

if __name__ == "__main__":
main(sys.argv)

Since nobody answered this question yet, i'm going to do my best to do it. Be aware i'm no expert on DSP.

I believe your first function is calculating the signal and noise average power based on each one's Fourier Transforms. So you're returning the SNR in decibels at the end.

\begin{equation}
SNR_{db} = 10 * \log_{10}{\frac{P_{signal}}{P_{noise}}}
\end{equation}

Your gaussian noise function generates the noise based on a scaling factor k of the signal max amplitude. Since you want to scale the amplitude of the noise based on your signal, i believe you want a relationship of:

\begin{equation}
k=\frac{A_{noise}}{A_{signal}}
\end{equation}

With each A meaning RMS amplitude. That can be generated used the following equation:

\begin{equation}
k = \sqrt{\frac{1}{SNR_{linear}}}
\end{equation}

or, to be more clear:

\begin{equation}
k = \sqrt{\frac{P_{noise}}{P_{signal}}} =\sqrt{\frac{A_{noise}^2}{A_{signal}^2}}
\end{equation}

If you were given the SNR in decibels and was asked to generate a noise based on it, you can use the following equation:

\begin{equation}
k = \frac{1}{10^{\frac{SNR_{db}}{10}}}
\end{equation}

In the third function you're generating the output signal by adding the frequency components of each signal, but if it's just an additive gaussian noise, you could just add the noise to the signal. (I'm not exactly sure on this).

You could also generate the linear SNR from your SNR in decibels, I've used this function in one of my projects once:

def linear_snr_from_db(snrdb):
    return 10.0 ** (snrdb/10.0)

I got all of the equations based on the wikipedia page for SNR and my little bit of experience.


Why Frequency Domain conversion is Important in Digital Image Processing? Any Image enhancement technique can be easily applied to the normal images. Still, most of the folks say that frequency domain conversion of the images can lead to better enhancement. Why do we need to view images or apply enhancements in frequency domain when they can be easily used in spatial-domain? Can someone please explain this in plain and simple English rather then mathematical equations for the sake of simply understanding the concept.

Not sure if frequency domain is necessarily better for image processing but frequency domain does give you an advantage for certain types of problems because it is more quantitive.

First you need to understand what frequency domain of an image means as it is intuitively difficult to look at frequency domain of an image and make sense out of it. Topic is beautifully handled here : What does frequency domain denote in case of images?

Basically frequency domain represents the rate of change of spatial pixels and hence gives an advantage when the problem you are dealing with relates to the rate of change of pixels which is very important in image processing. For example :high frequency in the frequency domain represents rapidly or sharply changing pixels such as boundaries or edges in an image. A high pass filter can be extremely helpful in identifying or removing these edges easily but the same problem is much more difficult in spatial domain (x-y domain).
Similarly a simple low pass filter can be used to get a smoother image. As you can guess images can be easily made sharper or smoother by manipulating different parts of frequency domain.

Some points on why it is so powerful :

Frequency domain gives you control over the whole images, where you can enhance(eg edges) and suppress (eg smooth shadow) different characteristics of the image very easily.
Frequency domain has a established suit of processes and tools that be borrowed directly from signal processing in other domains.
Some tools used for even image recognition such as correlation , convolution etc are much simpler and computationally cheaper in frequency domain.


Detect song changes in audio stream I have a stream of music that I'm reading(from a microphone, from an online streaming source or whatever). And I want to have a system that will detect when a song changes.

The question is, is this possible with any kind of accuracy?(even something 75% rate of detection would be helpful).

Alternatively, would it be possible to detect large enough changes in the music(so it might be a song change but it could just be a switch from a verse to chorus or chorus to bridge, etc...).

I am not asking for an exact algorithm -- though obviously if you know one then I obviously wouldn't mind -- but even a general approach on how to tackle such a problem would be helpful.

All pointers welcome, scientific papers, open source projects, some random guy's blog, whatever.

The very first thing to try would be to implement a rough speech/music classifier to detect a radio host/DJ; and to extract the amplitude envelope to detect fade outs and fade ins.

Once you have implemented that and/or if this doesn't solve your problem because your stream is continuous (DJ-style mix), you'll certainly run into the song boundary / song section boundary ambiguity - but if this is not a problem to you then there are solutions. What you are looking for is known in the academic community as ``unsupervised audio stream segmentation" - with a variety of techniques proposed over the years. You have to decide on two ingredients: which audio features to extract from the audio stream; and which statistical analysis method to use to detect changes/structure.

Regarding the features, the consensus is on a combinations of the following:

MFCC (which carry timbre/instrumentation information)
Chroma (which carry key/chord information)
Loudness
Tempo or other measures of rhythm (beat histogram, Pampalk's fluctuation patterns, etc.)

Very frequently, an algorithm can be repurposed just by a change of features... For example, there are papers on musicological analysis of piano music which will not use MFCC (same timbre throughout the piece), or papers on pop song structure analysis which will not use tempo (assumed to be constant throughout the piece, so not informative) - but the methods they describe could very well work for you if you decide to include other features - so don't be narrow when looking for literature.

Regarding the statistical analysis methods, the usual suspects are:

Computation of a self-similarity matrix, and "checkerboard" filter to detect blocks. The classic paper is by Foote, "Automatic Audio Segmentation Using A Measure of Audio Novelty". Foote's novelty detection function can be computed online and incrementally - over a buffer of past audio frames.
Hidden Markov models, used in an unsupervised way. You hope that the HMM states will map well into segments (Peeters and Rodet, "Music Structure Discovering Using Dynamic Audio Features for Audio Summary Generation: Sequence and State Approach"), or you overestimate the number of HMM states and look for changes in the HMM state distribution within small windows (Levy and Sandler, "Structural Segmentation of Musical Audio by Constrained Clustering"). These methods are more geared towards offline processing, and situations in which you have an idea in advance of how many segments you want to identify.
"Streaming" change detection methods. The principle is to split the past $N$ frames of audio into two $N/2$ frames blocks $S_1$ and $S_2$, and to perform a statistical test to decide which of the hypotheses "Frames from $S_1$ and $S_2$ are drawn from the same distribution $P$" and "Frames from $S_1$ and $S_2$ are drawn from a different distribution $P_1$ and $P_2$" is the most likely. Big peaks or dips in the test outcome indicate that the middle of the block is a segment change. There are a number of information criteria used to weigh the two hypotheses, such as the BIC (Classic paper: "Unsupervised Audio Stream Segmentation And Clustering Via The Bayesian Information Criterion" by Zhou). The limitation here is that you usually have little data to perform the test, and the data has weird-looking distributions. See the section 2.2 of Gillet et al's "Comparing audio and video segmentations for music videos indexing" for a few pointers on kernel methods that overcome these limitations.

For a real-life application (indexing of radio streams), see for example Ramona's "Combined supervised and unsupervised approaches for automatic segmentation of radiophonic audio streams" - this work combined the kind of unsupervised section detectors discussed here to a supervised speech/music classifier (to tackle the most obvious part of the problem).


How can I perform object tracking using noisy images? The problem

Take 50+ of Lena images, add gaussian noise to them to the point that you couldn't say it is Lena. I insist, gaussian noise, not gaussian blurring.

Now imagine a cropping window randomly moving around Lena's face (like a video camera focused on this part).

Now you have a handful of noisy pictures containing what used to be Lena's face somewhere in the image.

The question

Now I want to get back the face without knowing it in advance :)

First I have to identify overlapping parts (i.e. noisy faces here) from the set of images.
Then I must superimpose the faces.
And finally, add them to get back her beautiful face.

The problem is that I have no clue of how to do that (cross correlation, perhaps to identify the faces and the shift between them).
I've never had courses in image processing. I just need a reference or an algorithm name.

If the noise levels added to each pixel are independent, perhaps blurring the images first and then applying a pairwise cross-correlation could do the job.


Book recommendations on Gabor filter and Fourier series I'm starting to learn about Gabor filters and Fourier series. I need to make a presentation on Gabor filters in a few months, so I need quality references for the presentation.

Does anyone have any suggestions on introductory books, useful articles or other similar literature which I may use as references for theory on Gabor filters and Fourier transform?

This is a book that I would recommend as a starting point: Richard Szeliski: Computer Vision: Algorithms and Applications


Conditions for precoding matrix to preserve complex conjugate symmetry on DFT vector Suppose there is a DFT vector $\mathbf{X}$ with length N, which presents complex conjugate symmetry around its middle point, i.e., $X(1) = X(N-1)^*$, $X(2) = X(N - 2)^*$ and so forth.  $X(0)$ and $X(N/2)$ are the DC and Nyquist frequency respectively, therefore are real numbers. The remaining elements are complex.

Now, suppose there is a matrix $\mathbf{T}$, with size $N \times N$, which multiplies vector X.

\begin{align}
\mathbf{Y} = \mathbf{T}\mathbf{X}
\end{align}

The question is:

In what conditions, for matrix $\mathbf{T}$, the complex conjugte symmetry around the middle point of the resulting vector $\mathbf{Y}$ is preserved?

The motivation for this question is trying to come up with a precoder matrix $\mathbf{T}$ that results in a precoded (pre-equalized) symbol $\mathbf{Y}$ whose IFFT is real.

EDIT:

Thanks @MattL. and @niaren.
The difficulty about this question is to find necessary conditions. Matt's answer is indeed sufficient. It is also sufficient to make the following modifications:

The first row and the first column do not need to be zero. Instead, they could be non-zero, as long as its values present a complex conjugate symmetry around the middle point, its first value is real and its $(N/2+1)$-th value is real, just like the symbol. The same could be stated for the $(N/2+1)$-th column, the $(N/2+1)$-th row, and the main diagonal.

Secondly, the same correspondence between the matrix in the upper left corner and lower right corner could be made between the upper right corner and the lower left corner, that is, chose an $(N/2 -1)\times(N/2-1)$ matrix starting from $t_{2,N/2 + 2}$ to $t_{N/2,N}$, flip from left to right, flip upside down and take the conjugate, then put at the lower left corner. On MATLAB, that would be:

T(N/2+2:N,2:N/2) = conj(fliplr(flipud(Tisi(2:(N/2),N/2+2:N))))

This structure is similar to the structure of the DFT matrix. Would that be a necessary condition?

EDIT(2):

The following code implements such a valid operator for any real-valued $N \times N$ matrix $\mathbf{A}$:

N = 8;
A = rand(N,N); %must be real-valued
w = exp(-1j*2*pi/N); % twiddle factor
W = w.^(repmat(0:N-1,N,1).*repmat(0:N-1,N,1).'); % DFT matrix
T = W*A*W'

EDIT(3):

It is also interesting to note the $\mathbf{T}^{-1}$ presents the sufficient condition too. This comes from the fact that:

\begin{align}
\mathbf{T}^{-1} &= \mathbf{\left(W A W^{H}\right)}^{-1}\\
 &= \mathbf{\left(W^{H}\right)}^{-1} \mathbf{A}^{-1} \mathbf{W}^{-1}
\end{align}
where $\mathbf{W}$ is the DFT matrix.

Since $\mathbf{W^{H}} = N\mathbf{W}^{-1}$. This equation becomes:

\begin{align}
\mathbf{T}^{-1} &= \left(N\mathbf{W}^{-1}\right)^{-1} \mathbf{A}^{-1} \frac{1}{N}\mathbf{W^{H}}\\
&= \mathbf{W} \mathbf{A}^{-1} \mathbf{W^{H}}
\end{align}

Finally, since $\mathbf{A}^{-1}$ is real-valued, provided that $\mathbf{A}$ is full rank,$\mathbf{T}^{-1}$ is sufficient.

I think that the entries in your matrix $\bf{T}$ must obey $a_{N-n+1,N-m+1} = a^*_{n,m}$.
This is saying that the entries in row $N-n+1$ are the same as the coefficients in row n but where the coefficients are conjugated and reversed. The pattern in $\bf{T}$ for $N=4$ is

$T_4 = \left[ \begin{smallmatrix} a_{11}&a_{12}&a_{13}&a_{14}\\ a_{21}&a_{22}&a_{23}&a_{24} \\ a^*_{24}&a^*_{23}&a^*_{22}&a^*_{21} \\ a^*_{14}&a^*_{13}&a^*_{12}&a^*_{11}\end{smallmatrix} \right]$

I'm sure someone will come up with a better and more precise answer.


Help in analytical expression of BER for M-QAM I am new in learning deriving BER expression. I am familiar for the case of BPSK in AWGN but cannot understand how the BER expression in Eq(12) in this paper, An Enhanced Spectral Efficiency Chaos-Based Symbolic Dynamics Transceiver Design
download link

comes. http://www.raymaps.com/index.php/theoretical-ber-of-m-qam-in-rayleigh-fading/

tutorial presents the BER expression for M-QAM in AWGN.

Do you have any reference on how to analytically derive the BER for QAM and Rayleigh in AWGN for this case? How can I derive?

You need to read a good book to form a solid foundation first and then read research papers to understand advanced concepts. Eq (12) comes through a standard procedure for BER calculations. See Digital Communications by John Proakis, for example.


Why does the sound play faster I am trying to duplicate a stereo channel into four channels.
When I do that the sound plays twice as fast. it's PCM raw data
Basically trying to copy the stereo bytes into four channel buffer.
In Audacity
16kHz plays faster.
8KHz can hear it perfectly.

Let Sic be the i-th sample of the c-th channel, you started with stereo (2-channel) 16kHz. Your audio was like this
 s11, s12, s21, s22, s31, s32, s41, s42, s51, s52
When you made it a 4-channel audio you reassigned the same infromation to different channels like this
 s11, s12, s13, s14, s21, s22, s23, s24
In the stereo configuration you were using 2 samples every sampling period, so at 16kHz you play 32000 samples per second. With four channel you use 4 samples every sampling period, if you want to play 32000 samples per second your sampling rate needs to be 32000/4, 8kHz.


A few overview questions about FFT. Relationships between input, buffer size, and expected output? I have a FFT set up and I've found that it responds very differently based on the buffer size and input. The output varies wildly, with magnitude of the peak frequency ranging from 500 to <1.

I have two tests in particular I'm running, the first being a sine wave at amplitude 1.0 with an LFO sweeping it's frequency at 0.1Hz. I've found that smaller buffers are more accurate in their ability to track the frequency. The magnitude of the frequency playing reaches 300 pretty easily.

The next test is with a mic, and it responds incredibly differently. If I whistle or hum, it will pick up that frequency. If I play music, rather than getting an increase in magnitudes of all frequencies, everything becomes weaker, including any whistling or humming. Changing the buffer size to be smaller helps, but at most iterations of the FFT on a buffer it results in very tiny magnitudes across the board.

In general, the things I'm finding are that smaller buffers are more likely to get a very strong magnitude of a frequency, but the more frequencies there are with a (magnitude > 1), the less likely that any of them are very strong (> 200).  But none of this really means anything if it's a product of a broken implementation, and I can't seem to find any specific numeric information on what my expectations should be.

So, given an audio input and an output of frequencies and magnitudes (calculated by $\sqrt{\frac{\text{real}^2 + \text{imag}^2}{\frac{\text{numSamples}}{2}}}$), my questions are as follows:

Given a sole, unmoving sine wave, with a relatively small buffer size what kind of magnitudes are to be expected on what frequencies?
Given a sole, unmoving sine wave, with a relatively large buffer size what kind of magnitudes are to be expected on what frequencies?
Given a sole, sweeping frequency sine wave, with a relatively small buffer size what kind of magnitudes are to be expected on what frequencies?
Given a sole, sweeping frequency sine wave, with a relatively large buffer size what kind of magnitudes are to be expected on what frequencies?
Given static, with a relatively small buffer, what kind of magnitudes are to be expected?
Given static, with a relatively large buffer, what kind of magnitudes are to be expected?

I know this is a lot of questions, but the results I'm getting from these tests are not what I was expecting, and I want to start by making sure it's not my expectations that were wrong.

So, I'm not an expert at all, but I feel that I've perhaps learned enough to be somewhat helpful to anyone who ends up here with little introduction to FFT.  If any of this is wrong, see the "P.S" at the bottom.  tl;dr I'd rather share than leave this question unanswered entirely.

A few things are important to note: Parseval's theorem, bin frequencies, and normalization:

Parseval's theorem is something to the ends that the energy of the signal will be the energy of FFT bins (bins being each frequency/magnitude pair).  So, what this means in practice (I think) is that a small buffer size will make the bin magnitude of a sine wave frequency vary based off the amount of high amplitude wave is in the buffer.  It won't be consistent.

Bin Frequencies are the frequency of each bin from your FFT.  They aren't going to be the exact numbers you are looking for without a bit of work, but they should be close enough to give you some info.

Normalization is something that helps with all this.  For my project, I found the max bin magnitude, averaged it with the max from the last cycle, and then divided everything by it.  The result is that my visualization is showing something much more like what my ears feel they would see.

Now for the Q's (assuming all signals are normalized to -1 < x < 1):

Magnitudes of ~200 +-100, in the bins near the signal wave freq
Same, but magnitudes of ~300 +-50.  The larger buffer limits cases of more or less high amplitude chunking (sin(x) at [0 to Pi/4] vs [Pi/4 to Pi/2])
Magnitudes of ~150 +-100, in the bins near the signal wave freq.  The sweeping spreads out the bins a little, but the small buffer size (relative to sweep speed) helps.
Magnitudes of ~50 +-100, in the bins touched by the signal wave freq.  If the sweep is larger or faster expect the magnitudes to be even more spread out.  Note: the more bins touched, the lower the magnitudes overall.
Magnitudes of ~5 +-5, in all bins which happen to be strongly represented in this small buffer chunk
Magnitudes of ~1 +-1 in all bins.

So, it seems there is probably some form of log function which would help make the perceived levels of all these sounds not so distant when graphed by magnitude, but, overall, using normalization has worked in my project to pick out the most noteworthy note tones.  (btw notes and tones are anagrams)

P.S.  If any of this is wrong, please do make an answer covering and correcting these topics.  I will not be marking this as the accepted answer for at least a month because I'm not experienced enough to be an authority on this subject.  That being said, this is what I think may be true based off the results I've been getting, and it's better than nothing for a newcommer.


Help calculating / understanding the MFCCs: Mel-Frequency Cepstrum Coefficients I've been reading bits and pieces online but I just can't piece it all together. I have some background knowledge of signals / DSP stuff which should be enough prerequisites for this. I'm interested in eventually coding this algorithm in Java but I don't understand it fully yet which is why I'm here (it counts as math, right?).

Here's how I think it works along with the gaps in my knowledge.

Start with your audio speech sample, say a .wav file, that you can read into an array. Call this array $x[n]$, where $n$ ranges from $0, 1, \ldots ,N-1$ (so $N$ samples). The values correspond to audio intensity I guess - amplitudes.
Split the audio signal into distinct "frames" of 10ms or so where you assume the speech signal is "stationary". This is a form of quantization. So if your sample rate is 44.1KHz, 10ms is equal to 441 samples, or values of $x[n]$.
Do a Fourier transform (FFT for computation's sake). Now is this done on the entire signal or on each separate frame of $x[n]$? I think there's a difference because in general the Fourier transform looks at all elements of a signal, so $\mathcal F(x[n]) \neq \mathcal F(x_1[n])$ joined with $\mathcal F(x_2[n])$ joined with $\ldots \mathcal F(x_N[n])$ where $x_i[n]$ are the smaller frames. Anyways, say we do some FFT and end up with $X[k]$ for the rest of this.
Mapping to the Mel scale, and logging. I know how to convert regular frequency numbers to the Mel scale. For each $k$ of $X[k]$ (the "x-axis" if you'll allow me), you can do the formula here: http://en.wikipedia.org/wiki/Mel_scale. But how about the "y-values" or the amplitudes of $X[k]$? Do they just remain the same values but shifted to the appropriate spots on the new Mel (x-) axis? I saw in some paper there was something about logging the actual values of $X[k]$ because then if $X[k] = A[k]*B[k]$ where one of those signals is presumed to be noise you don't want, the log operation on this equation turns the multiplicative noise into additive noise, which hopefully can be filtered (?).
Now the final step is to take a DCT of your modified $X[k]$ from above (however it ended up getting modified). Then you take the amplitudes of this final result and those are your MFCCs. I read something about throwing away high frequency values.

So I'm trying to really iron out how to compute these guys step by step, and clearly some things are eluding me from above.

Also, I've heard about using "filter banks" (an array of band pass filters basically) and don't know if this refers to making the frames from the original signal, or maybe you make the frames after the FFT?

Lastly, there's something I saw about MFCCs having 13 coefficients?

Step by step...

1. & 2. This is correct. Note that the frames are usually overlapping, for example, frame 0 is samples 0 to 440; frame 1 is samples 220 to 660; frame 2 is samples 440 to 880 and so on... Note also that a window function is applied to the samples in the frame.

3. The Fourier transform is done for each frame. The motivation behind this is simple: a speech signal varies over time, but is stationary over short segments. You want to analyze each short segment individually - because on this segments the signal is simple enough to be described efficiently by few coefficients. Think of someone saying "hello". You don't want to see all phonemes collapsed into one single spectrum (FFT collapses temporal information) by analyzing all the sound at once. You want to see "hhhhheeeeeeeeeeelloooooooooo" to recognize the word stage by stage, so it has to be broken down into short segments.

4. "Mapping to the Mel scale" is misleading and that's probably why you're getting confused. A better description for this step would be: "Compute the signal energy through a bank of filters tuned to mel-scaled frequencies". Here is how this is done. We consider $N$ frequencies (a commonly used value is $N = 40$) equally spaced according to the mel scale, between 20 Hz (bottom of the hearing range) and the Nyquist frequency. Practical example: The signal is sampled at 8kHz and we want 40 bins. Since 4kHz (Nyquist) is 2250 mel, the filterbank center frequencies will be:
0 mel, 2250 / 39 mel, 2 x 2250 / 39 mel .. 2250 mel.

Once these frequencies have been defined, we compute a weighted sum of the FFT magnitudes (or energies) around each of these frequencies.

Look at the following picture, representing a filter bank with 12 bins:

The 8th bin has a center frequency of around 2kHz. The energy in the 8th bin is obtained by summing weighted FFT energies in the 1600 to 2800 Hz approximately - with the weight peaking at around 2kHz.

Implementation note: This bunch of weighted sums can be done in a single operation - a matrix multiplication of a "filterbank matrix" by the FFT energies vector.

So at this stage we have "summarized" the FFT spectrum into a set of 40 (12 in the illustration) energy values, each corresponding to a different range of frequencies. We take the log of these values.

5. The next step consists in taking the DCT of this sequence of 40 log-energies. This will yield 40 values. The first $K$ coefficients are the MFCC (Usually, $K = 13$). Actually, the very first DCT coefficient is the sum of all the log-energies computed at the previous step - so it is an overall measure of signal loudness and is not very informative as to the actual spectral content of the signal - it is often discarded for speech recognition or speaker id applications where the system has to be robust to loudness variations.


Explanation of fundamental filtering's consequences on signal Can anyone explain why exactly an "Overshooting" phenomena is observed when the fundamental harmonic is removed as seen on the figures? Is it technically right to call this "overshooting" at all ? If not, how is it referred to ?

Signals in time-domain:

Spectra:

As already mentioned by @Fat32 in the comments to your question what you're observing is not commonly referred to as overshoot.

As there are no units on both plots comparing the results in both figures is somewhat difficult. The following assumes that the horizontal axes in both plots are just the indices of either the time or frequency domain values (otherwise the graphs cannot be directly related without additional information and/or do not make sense).

The exact effects of what you're observing will generally depend on the frequency, amplitude and relative phase of the harmonics present in your signal.

A rough visual comparison of the green and blue lines in the first plot indicates that indeed the first harmonic has been removed. The difference between the green and blue lines is roughly a sine wave (phase zero) that fits into the window (i.e. no difference at the start, end and 1/2 of the window; largest difference at roughly 1/4 and 3/4 of the window).

Note that if the amplitude and/or phase of the first harmonic are different the time domain signal, resulting after subtraction of the first harmonic, will probably look completely different.


When calculating SNR, is noise included in signal? Given a signal $x(t) = s(t) + n(t)$ where $s(t)$ is the desired signal voltage and $n(t)$ is the noise, should the signal to noise ratio of this signal be 20log(xrms/nrms) or 20log(srms/nrms)? i.e., should I use the desired signal or the total signal for the signal to noise ratio?

The SNR is the ratio of $S$, the power in the signal $s(t)$, and $N$, the power in the noise $n(t)$. So, $$\text{SNR}_\text{dB} = 10\log_{10}\frac{S}{N}.$$

When using RMS values, we have $S = s_\text{RMS}^2/R$ and $N = n_\text{RMS}^2/R$. The resistances cancel out and then $$\text{SNR}_\text{dB} = 20\log_{10}\frac{s_\text{RMS}}{n_\text{RMS}}.$$


Can the Nyquist sample rate be extended to stochastic sampling? It appears there are lots of questions here about Nyquist, and a few questions about stochastic sampling here.  But I haven't found any that address quite what I'm after.  This is the closest I've found here, and searching around has led me to some work on Compressed Sensing, but I have no idea how to navigate that field.

So there's the Nyquist sampling rate: given a signal composed of frequencies lower than $F$, Nyquist gives a lower bound on sampling frequency necessary to reconstruct the signal.

But I think that's usually applied when trying to find an appropriate fixed sample rate.   So I'm curious if there is an extension to stochastic sampling: for the same signal as above, assume it is sampled such that the space between samples is distributed by Poisson $\mathcal{P}(\lambda)$ (or some other friendly distribution with parameters $\theta$.)

I'm hoping there exists a high probability bound on the ability to reconstruct the signal based on the sampling distribution.

Or, if signal reconstruction isn't well defined in this setting, perhaps there is an expression that can be derived that connects the source signal's frequency composition, the sampling distribution parameters, and the maximum mutual information between the source signal and the stochastically sampled signal.

But I have no idea if either actually exist.  Can anyone help me figure this out?

As long the mean sampling frequency fulfills the Nyquist criteria, you're good to go. So, for well behaved distributions with known mean, you won't need any extra math, aside from the required special handling of the non-uniform samples themselves.


Fresnel transform over image You can recommend the literature or tell you what the Fresnel Transform does to the image.
I realized that this describes the passage of a light wave through an optical element, but I need to understand more.

Within the context of optical image/signal processing, the Fresnel transform mathematically describes the manifestation of a physical phenomenon, known as diffraction, on the image formation process of a practical camera system using a finite aperture size.

Diffraction causes blurring on the image, and defines the physical limit of achievable image resolution, given an ideal lens system (which has no distortions of any kind on the image formation process). No better resolution is possible using that aperture size and optical imaging system.

Diffraction can be analysed through the Fourier transformation of the aperture function.


Implementing real-time signal processing filters I’m currently working on my capstone for school, and I’m looking for some guidance/tips when it comes to real-time filtering. The nature of the project is to monitor breathing patterns, send the aquired signals wirelessly, and then process/filter the data using a mobile app. My part of the project is the latter, and I would like some tips in terms of coding language/various articles that could guide me in the right direction! Thanks

There is a recent (2015) 150+ pages review called Smartphone-Based Real-Time Digital Signal Processing with lab codes. There are chapters like:

2: Android Software Development
3: iOS Software Development
6: Real-Time Filtering
7: Adaptive Filtering

  Real-time or applied digital signal processing courses are offered as
  follow-ups to conventional or theory-oriented digital signal
  processing courses in many engineering programs for the purpose of
  teaching students the technical know-how for putting signal processing
  algorithms or theory into practical use. These courses normally
  involve access to a teaching laboratory that is equipped with hardware
  boards, in particular DSP boards, together with their supporting
  software. A number of textbooks have been written discussing how to
  achieve real-time implementation on these hardware boards. This book
  discusses how smartphones can be used as hardware boards for real-time
  implementation of signal processing algorithms as an alternative to
  the hardware boards that are currently being used in signal processing
  teaching laboratories. The fact that mobile devices, in particular
  smartphones, have now become powerful processing platforms has led to
  the development of this book, thus enabling students to use their own
  smartphones to run signal processing algorithms in real-time
  considering that these days nearly all students possess smartphones.
  Changing the hardware platforms that are currently used in applied or
  real-time signal processing courses to smartphones creates a truly
  mobile laboratory experience or environment for students. In addition,
  it relieves the cost burden associated with using a dedicated signal
  processing board noting that the software development tools for
  smartphones are free of charge and are well-developed. This book is
  written in such a way that it can be used as a textbook for applied or
  real time digital signal processing courses offered at many
  universities. Ten lab experiments that are commonly encountered in
  such courses are covered in the book. This book is written primarily
  for those who are already familiar with signal processing concepts and
  are interested in their real-time and practical aspects. Similar to
  existing real-time courses, knowledge of C programming is assumed.
  This book can also be used as a self-study guide for those who wish to
  become familiar with signal processing app development on either
  Android or iPhone smartphones.


Processing 47-53 MHz signal sampled at 96000 Hz I am fairly new to the topic of DSP and FFT and I have hit a mental block. For a project, I am trying to pick out frequency's with amplitude peaks in a 47-53 MHz signal. This signal is being put through a 24 bit ADC with a sample rate of 96000 Hz. In order to successfully measure a sample, don't you have to have a sample rate >= to twice the desired frequency (Nyquist rate)? How can you get an accurate representation of a 47-53 MHz signal with such a low sample rate?

The second block I have hit involves the FFT of this signal. From my limited knowledge of FFT, the first "bin" (in the real part) represents 0 Hz and the last "bin" (in the real part) represents:
$$
F=i\times \frac{Fs}{N}
$$
Where $i$=Bin Index, $Fs$=Sample Rate, $N$=FFT Size

In my case with size 1024 and above sample rate, $511 \times 96000/1024 = 47906.25$ Hz, right? How can this range ($0$ Hz-$47906.25$ Hz) possibly represent my $47-53$ MHz signal? I get the feeling that I am missing something really crucial but I just can't figure it out.

Sorry if this is a really basic question and thanks for the help!

In short, what you're trying to do is impossible. The technique is known as band-pass sampling, or undersampling, but for it to work, you need that the signal's bandwidth be less than half your sampling frequency.

When you sample a signal with spectrum $S(f)$ and bandwidth $B$ at sampling frequency $f_s$, the sampled signal has spectrum $$S_s(f)=k\sum_i S(f-if_s)$$ where $k$ is a constant. If $B>f_s/2$ then the replicas $S(f-if_s)$ overlap and you end up with a distorted spectrum.

You are correct that the highest frequency bin in your FFT is 47906.25 Hz.


What does this notation $g(x) = S\{\delta(t) \}(x)$ mean? What does this notation mean $g(x) = S\{\delta(t) \}(x)$ mean, where $S$ is a "system", $\delta(t)$ is the impulse function and $g(x)$ the output function of the system? I am really not familiar with signals and signal processing. I know the definition of the impulse function

$$
\delta(t) =
\begin{cases}
1, \; t = 0,\\
0, \; t \neq 0
\end{cases}
$$

I understood that a system is a function that gets as input a function and returns another function (so it should be a functional). However, I do not get this notation $g(x) = S\{\delta(t) \}(x)$? For example, why the parentheses $\{ \}$? Why do we apply the $\delta(t)$ inside $S\{ \cdot \}$? What is the precedence of those function applications? How should I read that formula? Why do we even use an impulse here, what is the intuition behind this?

In my case, also have that

$$g(x) = S\{\delta(t) \}(x) = (A\delta(t))(x) = A_{x, 0}$$

Why are the notations $S\{\delta(t) \}(x)$ and $(A\delta(t))(x)$ equivalent? What is this $A_{x, 0}$? $A$ should be a Toeplitz matrix, but I don't understand why. How should this matrix $A$ look like?

$A_{x, 0}$ should represent coefficient associated to instant $0$ for computing the output at instant $x$. I also don't get this, even though I think we are trying to convert a continuous signal to a discrete one.

You can assume that $S$ is a linear shift-invariant system.

My understanding is that a system $S$ is a notation for   a generic object, which takes instances depending on its arguments, given between $\{\cdot\}$. Here, $S\{\delta(t)\}$ instantiates it as its effect on the zero-centered impulse response, with free variable $t$, taken on variable $x$.

Now, if this system $S$ is linear and shift-invariant, the operator given by $S\{\delta(t)\}$ defines $S$ it completely.

If the system applies on discretized inputs (that we can call $x$), it can be rewritten as a linear operator, represented by an infinite Toeplitz matrix. Each row implements the (linear) convolution, and the Toeplitz structure emulates the shift-invariance. I suggest to read Toeplitz and circulant matrices by R. Gray. The $0$ in $A_{x,0}$ may refer to the way an infinite matrix is matched with a signal centered at $0$, but I still wonder if this notation is clear.


nearest neighbor implementation matlab I am trying to implement this research paper:

I understand everything in the paper except the part about the nearest neighbors calculation (section 3).

The only MATLAB function I can think of that makes sense with this research paper squared euclidean distance is the pdist2 function. It would give me the D matrix.

Assuming I am correct, but what is the input to the pdist2 in terms of X and Y? Research paper mentions 100 nearest neighbors and sorted in ascending order, I am assuming this:

[D,I] = pdist2(X,Y,'squaredeuclidean','Smallest',100);

and X and Y is the magnitude spectrogram of the mixture signal. X = Y ?

Also, after I get D matrix, the paper mentions a P matrix of size n x p, I am not sure how to get that matrix from the previous calculation.

After I get that P matrix, I can do the rest of the research paper. So basically, I am just trying to figure out the nearest neighbors calculations D and P.

The author used MATLAB, and I would be happy if someone can provide MATLAB code for this nearest neighbor calculation, just like in the paper.

$D$ is $m \times m$, where $m$ is the total number of time frames in the spectrogram, it represents the squared euclidean distance of each frame to each other frame. In Matlab you could calculate it like this:

% Equation 10
for k = 1:m
    for l = 1:m
      D(k,l) = sum((X(:,k)-X(:,l)).^2);
    end
end

This is clearer from the text than the equation itself.

The $P$ matrix is calculated for each frame from the $p$ the nearest frames to it. The $Y$ matrix is then the median of $P$ at each frame. In matlab you could calculate it like this.

[~,I] = sort(D,2);
for k = 1:m
   pNearest = I(k,1:p);
   P = X(:,pNearest);
   Y(k,:) = median(P,2);
end


Hilbert Transform of a sine Function with Quadratic Argument $\sin(At^2 + Bt + \pi/4)$ I am looking for the Hilbert transform of the following function:

\begin{equation}
\mathcal{H}\bigg\{ \sin\Big(At^2 + Bt + \frac{\pi}{4}\Big) \bigg\}
\end{equation}

where $A$ and $B$ are constants with $A<0$ and $B>0$.

It is well known that $\mathcal{H}\{ \sin(Bt) \} = -\cos(Bt)$, which can easily be shown by rewriting the Hilbert transform as a convolution with $1/{\pi t}$ and using the spectral representation as shown hereinafter:

\begin{equation}
\mathcal{H}\{ \sin(Bt) \} = \sin(Bt) *\frac{1}{\pi t}
\end{equation}

where $*$ denotes the convolution operator. Consider the following Fourier pair:

\begin{equation}
\mathcal{F}\left\{\frac{1}{\pi t}\right\} = -\mathrm{j}\,\mathrm{sgn}(\omega)
\end{equation}

With this the problem can be solved in the spectrum as follows:

\begin{equation}
\mathcal{F}\big\{\sin(Bt) *\frac{1}{\pi t}\big\}
= \frac{\pi}{\mathrm{j}} \big(\delta(\omega-B) - \delta(\omega+B)\big) \; \big(-\mathrm{j}\,\mathrm{sgn}(\omega)\big)\\
%
= -\pi\, \big(\delta(\omega-B) + \delta(\omega+B)\big)
\end{equation}

Here, the two Dirac delta pulses are located at $+B$ and $-B$ angular frequency and hence the sign function is directly applied. Therefore, we get

\begin{equation}
\mathcal{H}\{ \sin(Bt) \} = \mathcal{F}^{-1}\Big\{ -\pi\, \big(\delta(\omega-B) + \delta(\omega+B)\big) \Big\} = -\cos(Bt)
\end{equation}

However, the same principle cannot be applied to $x(t) = \sin(At^2 + Bt + \pi/4)$, since its spectral function are two superimposed complex Gaussian functions also shifted to $+B$ and $-B$ angular frequency:

\begin{equation}
\mathcal{F}\{x(t)\} = \frac{\sqrt{-\frac{\pi}{A}}}{2 \mathrm{j}} \bigg(
\exp\Big(-\mathrm{j}\frac{1}{4A}(\omega-B)^2\Big) - \exp\Big(+\mathrm{j}\frac{1}{4A}(\omega+B)^2\Big)
\bigg) \quad \text{if} \quad A<0
\end{equation}

Each complex Gaussian function is defined for all frequencies and hence the application of the sign function does not simplify or solve the problem. I have also tried to directly solve the Hilbert transform integral with no success. I appreciate any help.

Directly calculating this seems difficult.

My argumentation is the following. For a signal $s(t)$, the so-called analytic signal $s_{\mathrm{an}}(t)$ can be obtained by
\begin{equation}
s_{\mathrm{an}}(t) = s(t)+\mathrm{j} \mathcal{H}\{s(t)\}\,\,, \text{where}\,S_{\mathrm{an}}(\omega) = 0\,\forall\, \omega < 0
\end{equation}
The analytic signal essentially corresponds to the spectral content of $s(t)$ in the positive frequencies only.

For your first example of a simple sinusoid, you can also come to the result for the Hilbert transform if you consider the analytic signal. The real sinusoid consists of the frequency components $\pm B$. The analytic signal should then be the component at $B$, which is then obviously a complex exponential, yielding the result for the Hilbert transform.

Now for your chirp signal $x(t)$, the situation is somewhat more complicated. If we think about a virtual "instantaneous frequency course" of the signal, it is
\begin{equation}
\omega_{x}(t) = \pm (2At+B)\,.
\end{equation}
Now this is somewhat strange, corresponding to two linearly changing components of opposite slope, crossing the zero-frequency point at $\omega_{x}(t=-\frac{B}{2A})=0$.

Now the analytic signal would have to represent the part of this frequency course above the zero-line of the $\omega-t$ plane (I might add some plots later). This means it would have to first have a negative slope, go down to frequency zero, and then abruptly change to a positive slope!

This means that the analytic signal would have to look something like
\begin{equation}
x_{\mathrm{an}}(t) = c_1 \exp{(-\mathrm{j} (At^2 + Bt + \frac{\pi}{4}))}\,\forall\,t < -\frac{B}{2A}
\end{equation}
and
\begin{equation}
x_{\mathrm{an}}(t) = c_2 \exp{(\mathrm{j} (At^2 + Bt + \frac{\pi}{4}))}\,\forall\,t \geq -\frac{B}{2A}\,\,,
\end{equation}
where the $c$ are some constant with $\vert c \vert = 1$.

Now we can determine the Hilbert transform of $x(t)$ by observing and checking in the equation for the analytic signal. This yields
\begin{equation}
\mathcal{H}\{x(t)\} = \cos{(At^2 + Bt + \frac{\pi}{4})}\,\forall\,t < -\frac{B}{2A}\,,\,\text{with}\,c_1 = -\mathrm{j}\,,
\end{equation}
and
\begin{equation}
\mathcal{H}\{x(t)\} = -\cos{(At^2 + Bt + \frac{\pi}{4})}\,\forall\,t \geq -\frac{B}{2A}\,,\,\text{with}\,c_2 = \mathrm{j}\,.
\end{equation}

One can probably also write these as one equation with the absolute value function. In any case, the point is the Hilbert transform seems to contain a discontinuity, which is what makes this particularly confusing to calculate I suspect.

I know it is somewhat "handwaivy", but I think the general idea/result is correct, so hope this helps!


Convolution: generating a continuous-time signal from a discrete-time sequence Given a discrete-time (DT) sequence $g[n]$, I want to represent it as a continuous-time (CT) signal. I can do this by representing this sequence as a weighted sum of Dirac delta impulses. Would it make a difference if I pass the DT signal through a DT filter first and then represent it as a weighted sum of Dirac impulses or pass the CT signal through a CT filter. The two cases are as follows:

Case 1: DT sequence converted first to CT signal and passed through CT filter.

The DT signal can be represented as a CT signal as:
$$g(t)=\sum_kg[k]\delta(t-kT)\tag{1}$$
If $g(t)$ in $(1)$ is the input to a continuous-time LTI system with impulse response $h(t)$, the output is given by
$$y(t)=\int h(\tau)g(t-\tau)d\tau = \int h(\tau)\left(\sum_kg[k]\delta(t-\tau-kT)\right)d\tau$$
or
$$y(t) =\sum_kg[k]\int h(\tau)\delta(t-\tau-kT) d\tau$$
or
$$y(t) = \sum_kg[k]h(t-kT)\tag{2}$$

Case 2: DT sequence passed through DT filter and then converted first to CT signal.

The DT signal is passed through a DT filter $h[n] = h(nT)$ to obtain a DT signal $z[n]$:
$$z[n]=\sum_kg[k]h[n-k]\tag{3}$$
This can be represented as a CT signal as
$$z(t)=\sum_kz[k]\delta(t-kT)\tag{4}$$

Is $z(t)$ same as $y(t)$? Can I express one in terms of the other? Any help wuld be greatly appreciated?
-ryan

I believe (2) and (4) are same. One can write $z(t)$ as

$$z(t) = \sum_k z[k]\delta(t-kT) = \sum_k \left(\sum_l g[l]h[k-l]\right)\delta(t-kT)$$

But $h[k-l] = h((k-l)T)$. Thus,
$$z(t) = \sum_l g[l] \left(\sum_k h((k-l)T) \delta(t-kT) \right) = \sum_l g[l] \left(\sum_k h(kT-lT) \delta(t-kT) \right) $$
or
$$z(t) = \sum_l g[l] h(t-lT)$$
which, by change of variables, is same as (2).


Intuition behind the Continuous Wavelet Transform? I was thinking sometime back about how to explain the Continuous Wavelet Transform ELI5. So this is what I came across.

The correlation of two exact signals is 1. So if I have an input signal $f(x)$ made up of an array of frequencies, how can I find out what frequencies exist at what points? Well, slide a signal $m(y)$ where $-\infty > y > +\infty$ over $f(x)$, and at those points where the co-relation of these signals is 1, well, those frequencies are present, at those times. This if of course a Continuous Wavelet Transform. Am I correct?

I think that the best way to explain the CWT is to start by explaining the Fourier Transform, then move on to explaining the Short-Time Fourier Transform, and then finally explain the CWT as a variation of the STFT.

The Fourier Transform exploits the fact that any decently behaved function can be represented as a sum of sinusoids (i.e. a Fourier series) and that the sinusoid basis posses the property of orthogonality:

$$ \int_{-\infty}^{\infty}  \sin(nx)\,\sin(mx)\,dx = \int_{-\infty}^{\infty}  \cos(nx)\,\cos(mx)\,dx = \begin{cases}
1,  & \text{if $n=m$} \\
0, & \text{if $n\neq m$}
\end{cases}$$

So, since:

$$ e^{iax} = \cos(ax)+i\,\sin(ax)$$

the Fourier Transform is simply doing this integration for all frequencies and keeping track of which outputs are zero (i.e. that frequency is not in the signal) and which are non-zeros (i.e. that frequency is in the data and its output is scaled by how much of it is in there):

$$ S(f) = \int_{-\infty}^{\infty} s(t)\,e^{i2\pi ft}dt $$

In this case you are doing this integration over the entire signal so you can't really tell if the frequency content is changing from the beginning of the signal to the end. One way around this is to compute the Short-Time Fourier Transform: i.e. window the signal, calculate the Fourier transform of the windowed signal, store it, then shift the window down a bit and repeat for all shifts:

$$ S(\tau,f) = \int_{-\infty}^{\infty} w(t- \tau)\,s(t)\,e^{i2\pi ft}dt  $$

where

$$ w(\tau,t)=\begin{cases}
1,  & \text{if $\tau\approx t$} \\
0, & \text{if $\tau \not\approx t$}
\end{cases} $$

The key thing here is that you are calculating the typical Fourier transform but of a new signal that only exists in a localized part of the t-axis. To emphasize this, you can see the new signal whose Fourier transform we are calculating by associating:

$$ S(\tau,f) = \int_{-\infty}^{\infty} [w(t-\tau)\,s(t)]\,e^{i2\pi ft}dt  $$

And here's a graphical example of this showing the new signal for different $\tau$ values and on the right are a few sinusoids to represent what we are using to decompose the signals (i.e. our basis, or kernel).

But, we can also change the association as such without changing the outcome:

$$ S(\tau,f) = \int_{-\infty}^{\infty} s(t)\,[w(t-\tau)\,e^{i2\pi ft}]dt  $$

So this means that instead of windowing our signal, we are windowing our basis functions. But here's the kicker, if we are windowing our basis functions, we don't have to use a constant-size window since we know that a basis function of high frequency will need a shorter window than a basis function of low frequency. This is the whole point of the CWT. It is a decomposition of a signal by "wavelets" (i.e. windowed sinusoids in this case) where the windowing is adaptive to the sinusoid frequency. If the we choose a Gaussian window (as I have chosen in these examples), then our wavelets are called Morlet wavelets (or Gabor wavelets, in some literature).

Finally, you can generalize this for any choice of wavelet that you want. In that generalization, you can describe your wavelet basis functions as being "stretched" and "squeezed" versions of some arbitrary "mother" wavelet, $ \psi $. And so the previous equation can now be written as the final form of the CWT:

$$ S(a,b) = \frac{1}{\sqrt b}\int_{-\infty}^{\infty} s(t)\,\psi (\frac{t-a}{b})dt $$

where $a$ is now what used to be our $\tau$ (i.e. a time shift), and $b$ is called the "scale" which is just a parameter to stretch and squeeze the wavelet (similar to what our parameter $f$ was except now the interpretation is more difficult). And the only reason you have the $\frac{1}{\sqrt b}$ up front is to normalize the wavelets so that they have the same "energy" and you end up comparing apples to apples in your time-frequency representation.

I hope this helps!

-Antonio


Minimum phase All-pass Why is it that reflecting any poles or zeros of a rational function gives a minimum phase system? And why is doing that make a unique minimum phase system?
I understand the all-pass function absorbs the poles or zeros when you reflect it, not sure if that has something to do with it. Also, wouldn't reflecting any poles or zeros potentially cause some poles and zeros that were before in the stability region to be out of it?

An LTI system is said to be minimum-phase if the system and its inverse system are causal and stable. That's implying that all poles and zeros must be strictly inside the unit circle.
An all-pass filter's zero-pole pairs are reflections of each other across the unit circle so that an all-pass filter absorbs the poles or zeros. For a stable all-pass filter, its zeros are outside the unit circle and its poles are inside the unit circles.
A causal and stable mixed-phase system has all poles inside the unit circle and zeros which could be inside or outside the unit circle. Therefore, reflecting zeros outside the unit circle by cascading an all-pass system makes all zeros inside the unit circle, which means you get a minimum-phase system.
Once you get a minimum-phase system, the only thing you can do to prevent changing its magnitude response is to cascade an all-pass system. That brings additional zeros outside the unit circle and the system is no longer minimum-phase. Therefore minimum-phase system with a specific magnitude response is unique.


Minimizing Time Sidelobes with Pulse Compression I am trying to compute a compression filter function to minimize the error from a a desired pulse compression result.

The usual approach to doing this is to minimize the RMS error of the the compressed signal while it's represented in linear terms (i.e., in volts). This method handles complex signals and pulse compression for my nonlinear FM signal, and appears to work as advertised. (The method is described in the end notes of this question.)

While the method minimizes the RMS error, it leaves some individual samples larger than I'd like them to be, while producing some that are much lower than they need to be. What I'd like is a time sidelobe response that is more uniform in the time domain. Something like the sidelobes of a Chebyshev Type II filter function is what I'm looking for, but two-sided and in the time domain, not the frequency domain:

This would be much more useful for my purpose (which is for a radar application).

My first thought is to perform the RMS calculations in decibels, rather than as linear values. If I could do that, the same optimization function can be used to solve the liner equations for minimizing the error. It may work or it might not, but that's all I can think of at the moment. One problem is that I don't know how I would handle the complex signals and the complex compression function if I do this in decibels. So I'm stuck at the moment.

Can someone give my some guidelines or a solution to achieve what I'm looking for?

End Notes:

Here is the general method for minimizing the RMS error:

Given a complex pulse
$$x(n), \ \ 0 \leq n \lt N$$
and a pulse compression function $$\omega(k), \ \ 0 \leq k \lt K$$
we have the compressed pulse
$$z(n)=\sum_{k=0}^{K-1}x(n-k)\omega(k), \ \ 0\leq n \lt (N+K-2)$$

If the desired compression results is
$$d(n), \ \ 0\leq n \lt (N+K-2)$$
and the error is defined as
$$\epsilon(n)=z(n)-d(n)$$
the residues is deifned as
$$\epsilon^2=\sum_{n=0}^{N+K-2} |d(n)-z(n))|^2 = \sum_{n=0}^{N+K-2} |d(n)-\sum_{k=0}^{K-1}x(n-k)\omega(k)|^2$$

If $Z$ and $D$ are column vectors for $z()$ and $d()$, we can write
$$\epsilon^2=(D-Z)^{*}\cdot(D-Z)= D^*D-D^*Z-Z^*D+Z^*Z$$

$Z$ can be written as
$$Z=\begin{bmatrix}
X & 0 & 0      & 0 & 0 \\
0 & X & \vdots & 0 & 0 \\
X & 0 & \vdots & 0 & 0 \\
0 & 0 & \vdots & X & 0 \\
0 & 0 & 0      & 0 & X
\end{bmatrix}\cdot W = M \cdot W$$
$X$ and $W$ are column vectors for $x()$ and $\omega()$ and $M$ is a $(N+K-1)$ by $K$ matrix.

I'll omit the intermediary steps, and conclude that $\epsilon^2$ is minimized with respect to $w$ when
$$W = (M^{*} \cdot M)^{-1} \cdot M^{*} \cdot D$$

FOLLOW UP
Not sure if this is the best place to put this, but the recommendation from Matt L. to use iterative reweighted least squares did produce some improvement.

Using a weight vector $p(n)$ (I already am using $w$ and $W$), this converts to the following solution in matrix form:
$$W = (M^{*} \cdot P \cdot M)^{-1} \cdot M^{*} \cdot P \cdot D$$

And here is an example of the result for a specific case, using the matched filter, the first iteration, and the 50th iteration. What I see is great improvement in the lobes I wanted to reduce. But I had hoped it would level out the sidelobes to something more even. A few more iterations led to a singular matrix problem, so I didn't go any further.

You could define a weight sequence $w(n)$ that weighs the error differently for different time indices $n$:

$$\epsilon^2=\sum_{n=0}^{N+K-2} w(n)\big|d(n)-z(n)\big|^2\tag{1}$$

Samples that are too large get a higher weight, whereas samples that are smaller than necessary can get a lower weight.
Of course, you cannot know in advance which samples should get which weight, so you need to iteratively solve the RMS problem minimizing $(1)$, and adapt the weight sequence $w(n)$ in each iteration depending on the magnitude of the error $e(n)=d(n)-z(n)$.

A possible update rule for the weights $w(n)$ is

$$w^{(k)}(n)=\frac{\left|e^{(k)}(n)\right|}{\displaystyle\frac{1}{N+K-1}\sum_n\left|e^{(k)}(n)\right|}\cdot w^{(k-1)}(n)\tag{2}$$

where $k$ is the iteration index. You probably need to experiment a bit with the update equation to adapt it to your desired result, but the general idea is to increase the weight where the resulting error is too large and decrease it where it is smaller than needed.

This type of algorithm is called iteratively reweighted least squares.


Passband vs Baseband Bandwidth Bandwidth is the difference between the upper and lower frequencies in
  a continuous band of frequencies. A key characteristic of bandwidth is
  that any band of a given width can carry the same amount of
  information, regardless of where that band is located in the frequency
  spectrum.

So I understood bandwidth to be a measure of how much space it takes up in the frequency domain. Am I right in saying that because of symmetry about the Y axis, the same information is repeated for negative frequency values, this is why the "bandwidth" of the bandlimited signal given below is not $2B = B - (-B)$ but $B$.

1) Since the left is a mirror image of the right part only for real signals, if the signal was not real, can we still say bandwidth is $B$?

2) Why do we have two different definitions for passband and baseband bandwidth? Is there any other reason?

  Passband bandwidth is the difference between the upper and lower cutoff
  frequencies of, for example, a band-pass filter, a communication
  channel, or a signal spectrum. Baseband bandwidth applies to a
  low-pass filter or baseband signal; the bandwidth is equal to its
  upper cutoff frequency.

3) Am I right in understanding that in case of either type of signals, alias-free signal sampling is possible at $2B$ samples per second?

You can unify the definitions for baseband and passband bandwidths by saying that bandwidth is only measured at positive frequencies. So for low pass signals, the lower frequency is zero and the upper frequency is the cut-off frequency. For bandpass signals, you have the difference between upper and lower cut-off frequencies (at positive frequencies).

For complex signals, a reasonable definition of bandwidth should take into account that the negative frequencies are not redundant anymore. So a complex low pass signal extending from $f=-B$ to $f=B$ should have bandwidth $2B$.

Concerning the choice of sampling frequency for bandpass signals, take a look at this answer. It is not true that a sampling frequency equal to twice the bandwidth is always sufficient. Yet you can usually do better than just sampling with twice the upper cut-off frequency.


Support of the convolution of two rectangular signals I'm trying to convolve two rectangular signals in the frequency domain $$H_1(\omega) = u[\omega +.2\pi] - u[\omega -.2\pi]$$ and $$H_2(\omega) = u[\omega +.1\pi] - u[\omega -.1\pi]$$

My result is a trapezoid that ranges from $[-.2\pi, .4\pi]$, but I am expecting a trapezoid centered at 0 with range $[-.3\pi, .3\pi]$.

I am pretty sure I'm getting something wrong with the bounds of integration, but I'm not sure what I'm doing wrong. My strategy is to flip $H_{2}$ and slide it from negative infinity through $H_{1}$. What am I missing in determining the bounds correctly?

We have rectangular pulses

$$H_1 (\omega) = u \left(\omega + \frac{\pi}{5}\right) - u \left(\omega - \frac{\pi}{5}\right)$$

and

$$H_2 (\omega) = u \left(\omega + \frac{\pi}{10}\right) - u \left(\omega - \frac{\pi}{10}\right)$$

and would like to convolve them. When convolving piecewise constant functions, a useful "trick" is to differentiate, convolve and then integrate. Differentiating $H_1$,

$$H_1 ' (\omega) = \delta \left(\omega + \frac{\pi}{5}\right) - \delta \left(\omega - \frac{\pi}{5}\right)$$

and, convolving with $H_2$,

$$(H_1' * H_2) (\omega) = H_2 \left(\omega + \frac{\pi}{5}\right) - H_2 \left(\omega - \frac{\pi}{5}\right) = \begin{cases} \,~~ 1, & \text{if } \omega \in \left[ -\frac{3 \pi}{10}, -\frac{\pi}{10} \right]\\ -1, & \text{if } \omega \in \left[ \frac{\pi}{10}, \frac{3\pi}{10} \right]\\ \,~~0, & \text{otherwise}\end{cases}$$

Finally, integrating,

$$(H_1 * H_2) (\omega) = \begin{cases} \, \frac{3 \pi}{10} + \omega, & \text{if } \omega \in \left[ -\frac{3 \pi}{10}, -\frac{\pi}{10} \right]\\ \qquad\frac{\pi}{5}, & \text{if } \omega \in \left[ -\frac{\pi}{10}, \frac{\pi}{10} \right]\\ \, \frac{3 \pi}{10} - \omega, & \text{if } \omega \in \left[ \frac{\pi}{10}, \frac{3\pi}{10} \right]\\ \,\, \qquad 0, & \text{otherwise}\end{cases}$$


Unexpected result after calculating the response of a system in frequency domain I'm implementing a Python script that calculates the response of a system given an input. The system is $y(t) = x(t-2)$, i.e. it delays the signal by $2$, and the signal is $x(t) = \sin(3t)u(t)$.

As I'm learning the discrete Fourier transform (and it's inverse), I'm doing all the calculation in frequency domain. The procedure is very simple:

Sample $x(t)$ between $[-10, 10]$ where the interval between each sample is $0.01$
Calculate the DFT of $x(t)$ sampled (the result is $X(j\omega)$)
Multiply the $X(j\omega)$ by $e^{-j2\omega}$ to shift the signal (the result is $Y(j\omega)$)
Calculate the inverse DFT to obtain $y(t)$ sampled

In the end, the result is plotted and it's straightforward to check its correctness for $t > 0$. However, for $t < 0$, the output has some kind of strange behavior, as you can see below:

Am I missing something? Can someone please explain me what causes this?

That's not strange – that's expected!

The DFT is a circular operation – by shifting to the right, your signal "wrapped" around.

What you need to do is called fast convolution based on the DFT/FFT. The Wikipedia has a nice article on the Overlap-Add and the Overlap-Save methods, so I'd recommend starting with these, and then going for the literature that suits your understanding best.


SIFT: why it is invariant to affine tranlsation? I am studying SIFT. I am a little confused in some of Lowe's algorithm.

In order to have a descriptor invariant to translations like scaling, we must get rid of the scaling factor in LoG (Laplacian of Gaussian). I know that the LoG is represented like this:

Then the scale invariant LoG would look like this:

But how the DoG can be like this?:

In other words, how by approximation,using DoG, we get rid of that scale factor in SIFT. I know how DoG is computed, but the thing is I cannot understand how we get rid of that factor to have a invariant-to-scale LoG approximation.
Any precise explanation is really appreciated.

EDIT: some more questions, Please answer precisely. Thanks :-)
1- At the final step of SIFT, we calculate feature vector using Gradient Orientation. to do do, the keypoint's rotation is subtracted from each orientation. I cannot understand why we subtract them?

2- What about Illumination dependence?

Anyone can explain how we remove Illumination and Orientation dependency?

I consider that you already know the proof of how we get DoG from LoG, so I won't go further with it.

According to your question "how by approximation,using DoG, we get rid of that scale factor in SIFT". DoG works as a feature detector to detect blobs in the given image (blob like shapes) which you can say it is a bit kind of edge detectors here is a photo of some blobs detected using DoG:

These blobs detected at different scales using the scale space pyramid. Using larger values of sigma you get bigger blobs. Experimentally DoG give you the best notion of scale as you can see here :

So, experimentally DoG gives you the most robust features that are scale invariant by constructing the scale space and computing the DoG and finally extracting local maxima points.

Useful links:

1- Lecture PPT

2- Best known video lecture for SIFT

NOTE: This not a professional answer and I hope to hear more from other fellows


Echo detection for localization I was reading about active sonar and couldn't quite understand how the echos are recognized.

My doubt is, when the sound signal is reflected back, it is greatly attenuated and distorted due to many factors.

Do they use MFCC(which I think is only good for audible range) for identifying the echo or are there any other better techniques? I feel that the problem is modeled and the echo data is fit to the model but I am not sure.

My doubt is, when the sound signal is reflected back, it is greatly attenuated and distorted due to many factors.

That's true!

  Do they use MFCC(which I think is only good for audible range)

I don't see how MFCC come into play here, at all. Radar/sonar applications typically depend on some kind of comparison between TX and RX, and a power spectrum only observation can't help you at all here.

  for identifying the echo or are there any other better techniques?

yes, there are – for example, simple cross-correlation would do.


How to reconstruct signal of its phase and magnitude functions? I have two continuous periodic (a period of $2\pi$) functions which belong to the phase and magnitude of Fourier transform of a signal, how can I reconstruct the original signal? What kind of transform should I use? (Basic Fourier transform or DTFT)

You can look at the domains covered by the Fourier transform (FT) or the DTFT (Discrete Time Fourier Transform). The DTFT applies to discrete signals, and yields a (possibly $2\pi$) periodic transformation. Here, $\omega$ denotes the continuous normalized radian frequency variable. So:

$$ (s[k])_{k\in \mathbb{Z}} \stackrel{\mathrm{DTFT}}{\longrightarrow} S(\omega) = \sum_{k=-\infty}^{\infty} s[k] z^{-k}, \quad z=e^{\imath \omega}$$
while for the FT:

$$s(t) \stackrel{\mathrm{FT}}{\longrightarrow} S(\omega) = \int_{-\infty}^{\infty} s(t) e^{- \imath \omega  t}dt$$
with no specific periodicity in general.

So most likely, you are asked to recover a discrete signal. The complex quantity $S(\omega)$ can be recovered with  the amplitude and the phase you have. You can then recover the discrete signal using the inverse DTFT over one period:

$$s[k] =  \int_{-\pi}^{\pi} S(w)e^{\imath \omega k} d\omega \,.$$

You definitely need to check the  Discrete Time Fourier Transform (DTFT) and the associated pages.


Keypoint position and scale in SIFT If I understood it right, a keypoint is a tuple $$(x,y,\sigma, r),$$ where $x,y$ define the position of the keypoint and $r$ an orientation (given by the most domiant gradiant vector around the keypoint). I'm still not sure about the intention of $\sigma$.

What is done in the first part of SIFT is to consider $$L(x,y,\sigma):=G(x,y,\sigma) \ast I(x,y),$$
where $I(x,y)$ denote the value of the image at position $(x,y)$, $\ast$ the convolution and $G(x,y,\sigma)$ the gaussian filter with standard-deviation $\sigma$, evaluated at position $(x,y)$.

Now DoGs are considered, for different sizes of the image and for different blurring-values $\sigma$ and interesting points are collected.

What confuses me most about this procedure is that Lowe denotes $\sigma$ as the scale of the image. However, $\sigma$ is just a blurring-parameter.

Hence, I define a blured image of a scaled image by
$$L(x,y,\sigma,s):=G(x,y,\sigma) \ast I_s(x,y),$$ where $I_s$ is the Image $I$ scaled by factor $s$ (So multiplying both dimensions of the image by s, and rescaling it). Thus, an interesting point in a DoG is a point $(x,y)$ in the image $L(x,y,k\sigma,s)-L(x,y,\sigma,s)$.

In Lowe's paper, he never considers the parameter 's' but just $\sigma$. Why is that possible ? How do I relate $(x,y)$ of $L(x,y,\sigma,s)$ to a position in the image $I$ and why do I save the blurring-parameter $\sigma$?

In what way is the keypoint now scale (resize)-invariant ?

Should there be in the resized image at some DoG a keypoint, which corresponds to the former keypoint by having more or less the same orientation and magnitude of its surrounding gradiants ? Is that the invariance Lowe is talking about ?

To find a keypoints, a local extrema of $L(x,y,\sigma)$ around a point $(x,y,\sigma)$ is considered. What is idea behind using interpolated extrema via Tayler expansion ?

Thanks Michals

The $\sigma$ parameter is both. The Gaussian function can generate a scale-space where $\sigma$ is the scale parameter. It doesn't mean the image is scaled, instead it is the scale at which the features are being evaluated.

For example, with higher $\sigma$ the image is more blurred and therefore only larger image features contribute to the gradient histogram. SIFT is made scale-invariant by choosing the greatest feature response from a set of responses from different scales. The Gaussian function is used because with increasing scale one can guarantee no new local extrema (features) are introduced.

I'm not sure about the interpolated extrema, but I guess it is to obtain sub-pixel accuracy.


Stable gain nyquist plot From Modern Control Engineering 5th Edition page 469

The transfer function of a plant controller by a proportional controller is given by

$$G(s) = \frac{K(s+0.5)}{s^3 + s^2 + 1}$$

In the book $G(j\omega) / K$ is plotted and from that they infer that for $ 2 < K$ the system is stable.

I plotted the nyquist plot of the plant itself

and tried getting the answer consistent with the procedure from our lecture.

I know that I need 2 counter clockwise encirclements of $(-1,0j)$ so I want to see for which gain the plot is shifted past -1.

$$\frac{-1}{K} > -0.5$$

$$K > 2$$

So far so good. They did not compute the upper limit of $K$.

I think I can see that the encirclements occur between -0.5 and 0 from the plot.

So next I need to calculate

$$\frac{-1}{K} < 0$$

$$\frac{1}{K} > 0$$

$$K > \infty$$

and this can't be. I assumed to get to $$K < \infty$$

I think my main problem is that I don't exactly know when to use < and > and which points (-0.5,0,0.5) I should choose for my calculations. Also how would I proceed if there was only 1 RHP and I'd need only 1 ccw encirclement.

Could someone clarify the general procedure for this or guide me to a ressource that explains it for beginners?

There is a mistake in the last step. When you invert both sides of the inequality, the $>$ inverts as well. So, for $K>0$ (as you found that $K>2$, so it must be positive):

$$\frac 1 K > y \iff K<\frac 1 y$$

If you do this when $y\to 0$, then you get the desired result.

Also, note that $\frac 1 K >0$ would hold for any $K>0$, as the only way to get a negative number from $1/K$ is letting $K$ be negative.

Regarding your other question, it would completely depend on the transfer function and its respective Nyquist. But don't worry, the procedure is exactly the same. For example, try doing the Nyquist plot for $-\frac{K}{1-s}$ and see what you get.