Correct method for drawing waveforms I need to draw waveforms for biometric data like ECG and EEG signals. When I have more samples than pixels at the X-axis, I need to draw a vertical line between the MIN and MAX sample-value for that period.

Now I was wondering, suppose the signal is highly volatile, and the MAX for the next point is lower than the MIN from the previous point, there will be a 'gap' in the waveform.

What is the correct way to connect these points? When drawn using standard polyline routines the pen will draw from the last position (MAX) to the start of the new line (MIN), which is a much larger distance, causing the signal look thicker than it would be if I just connected the previous MIN to the new MAX (if the trend is downwards), or the previous MAX to the new MIN (if the trend is upwards).

Is there a standardized way to do this, or is it just a matter of taste?

Remember that the signal is not the samples; the signal passes through the samples and exists between them, too.  So the maximum value of one chunk cannot be less than the minimum value of the next; they can only be equal, as the signal passes from a lower sample in one chunk to a higher sample in the next chunk.

Plotting a solid line from the min of a chunk to the max sample value of a chunk is not a very accurate representation of the signal anyway (because of intersample peaks and not showing the density of the wave).  If that's good enough of a representation for you, I'd just approximate it by drawing to the midpoint in this case:

find max of chunk 1
find min of chunk 2
if max(chunk 1) < min(chunk 2):
    in column 1: draw line from min of chunk 1 to average(max(chunk 1), min(chunk 2))
    in column 2: draw line from average(max(chunk 1), min(chunk 2)) to max(chunk 2)

If you want to get fancier, you could upsample the signal to get a better idea of the intersample peaks, and then do the same thing based on those.

And the best is to use brightness to show the density of the wave at different amplitudes, but that's significantly more complicated: https://flic.kr/p/7S8oHA


The effect of upsampling on DFT coefficients I am learning DSP by myself and I encountered a problem that bewilders me.
If I have a sequence of length N, and I upsample it by a factor of 3. How would the DFT change or related?

For example:

x = [1 2 3 4 5 6 7];

y = [1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7];

What is the difference between $X[k]$ and $Y[k]$?

Note that if you upsample a sequence $x[n]$ by a factor $M$ you have

$$y[Mn]=x[n]\tag{1}$$

The DFT of $x[n]$ is

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}nk}\tag{2}$$

The DFT length of $y[n]$ is $MN$, so the DFT of $y[n]$ is

$$Y[k]=\sum_{n=0}^{MN-1}y[n]e^{-j\frac{2\pi}{MN}nk}\tag{3}$$

Since only every $M^{th}$ sample of $y[n]$ is non-zero, $(3)$ can be written as

$$Y[k]=\sum_{n=0}^{N-1}y[Mn]e^{-j\frac{2\pi}{MN}Mnk}=\sum_{n=0}^{N-1}y[Mn]e^{-j\frac{2\pi}{N}nk}\tag{4}$$

With $(1)$ we get

$$Y[k]=X[k]\tag{5}$$

Note that from $(2)$ $X[k]$ is $N$-periodic, so the length $MN$ DFT of $y[n]$ is just an $M$-fold repetition of the DFT coefficients of $x[n]$.

A simple example:

x = [1,1];
X = fft(x)
ans =

   2   0

y = [1,0,1,0];
fft(y)
ans =

   2   0   2   0

y=[1,0,0,1,0,0];
fft(y)
ans =

   2   0   2   0   2   0


How can I prove stability of a biquad filter with non-zero initial conditions Ok, so the situation is that I have a DFII biquad with some filter coefficients:

  \begin{align}
w[n] &= x[n] - a_1*w[n-1] - a_2*w[n-2]\\
y[n] &= b_0*w[n] + b_1*w[n-1] + b_2*w[n-2]
\end{align}

While the filter is running, I change the coefficients and I believe my filter is going unstable.  I want to build in a check that it is going to be stable.  Currently, I simply check that the poles lie inside the unit circle by checking my $a_n$ coefficients.  But it doesn't work.  I believe that the problem is that I have non-zero values in my states (i.e. $w[-1]$ and $w[-2]$) and my stability check doesn't even consider this.  I cannot reset my states while the filter is running without producing an audible click.

The questions: How do I check for stability of a filter with non-zero initial conditions?

I have a couple thoughts for solutions:

Simply take the $\mathcal Z$-tranform of the impulse response (of the feedback portion), which looks something like this:
$$
h[n] =
\begin{cases}
    0 & n < 0\\
    1 - a_1*K_1 - a_2*K_2    & n = 0\\
    0 - a_1*h[0] - a_2*K_1   & n = 1\\
    - a_1*h[n-1] - a_2*h[n-2] & n \ge 2
\end{cases}
$$
Solve the difference equation à la differential equations:
$$
h[n] + a_1*h[n-1] + a_2*h[n-2] = 0, n \ge 2
$$
Or prove that the sequence $h[n]$ (from above) converges
$$ \sum_{n=0}^{\infty}\lvert h[n]\rvert  < \infty $$

Am I on the right track? I can't seem to figure out the solution for any of these methods.

Tricky problem. Not sure I can answer this but here are a few pointers:

Direct Form II is the worst biquad for audio processing. The transfer function between the input and the state is given just by the poles. The gain can be really large, I have seen gains of in excess of a 100 dB for reasonable audio filters.
This makes real time switching of Direct Form II very tricky. The magnitude of your state can change 10s of dB with seemingly small overall filter changes. This creates massive clicks and pops and occasionally unstable filters.
Much better choices are Direct Form I or transposed Form II. For these filter types the state variables are always in the same order of magnitude than input and output. Much less potential for loud pops and clicks or instability
There is not "sure fire" way of updating an IIR filter in real time without the potential for audible artifacts. Some potential methods are running two filters in parallel and slowly cross fading the output or updating pole and zero locations really slowly with very small steps at the time
If you must update in real time, make sure that you don't update one coefficient at a time but ideally both poles and zeros in one go. At least do the poles all at ones (i.e. always update a1 and a2 at the same time)
Update lower gain first. If the pole section of the new filter has lower gain, update the poles first and then the zeros. If the pole section of the filter has higher gain, update zeros first and then poles
If you need to stick with Direct Form II, it's probably best to zero out the state during update. This will still create an artifact but it's probably less than if you don't and at least it guarantees stability.


Difference between autocorrelations M. H. Hayes calls in his book "Statistical digital signal processing" autocorrelation sequences $r[k]$.

For optimal filters the desired autocorrelation has the index d -> $r_d[k]$.
However, often there is also an $r_{dx}[k]$. Are they the same? Or is there an assumption where you can interchange them?

Because in some cases it seems like he does just that. No clue why the index changed in some formulas.

This is how d and x are connected:

All these notation should be considered within the context of the development.

So $r_d[k]$ may refer to the auto-correlations sequence of the desired signal whereas the $r_{dx}[k]$ is the cross correlation between the input $x[n]$ and desired signal $d[n]$. So they are two different things.

The latter is used in Wiener optimal filtering to solve the Wiener-Hopf equations (aka The Normal Equations) but I cannot recall where the former is used, can you cite the formula, chapter etc. that makes use of it?


What does "convergence" mean for Thinning? I've been using Thinning recently on line-like structures in an attempt to generate skeleton-like chains of pixels connected through their Moore neighborhoods.

The description of Thinning says that the algorithm "finds the skeletons of foreground regions in image by applying morphological thinning until convergence".  What does "convergence" mean?  If I start with a cluster of pixels connected through their Moore neighborhoods, what am I guaranteed after "convergence"?  Will the product always be a chain of pixels also connected through their Moore neighborhoods?

The thinning operation become idempotent after a number of iterations, which means the result does not change anymore - in other word, the algorithm have converged to a solution.


Is jointly wss (wide sense stationary) a transitive relation? I've been try to either prove or find a counter-example to the idea of jointly-wss being transitive.
In other words: does ($x$ and $y$ are jointly wss) $\wedge$ ($y$ and $z$ are jointly wss) imply that $x$ and $z$ are jointly wss?
Obviously it boils down to only asking about the cross-correlation condition, so one can ask:
Does $\forall t, \Delta t : C_x,_y(t, t+\Delta t)=C_x,_y(0, \Delta t) \ and \   C_y,_z(t, t+\Delta t)=C_y,_z(0, \Delta t)$ imply that $\forall t, \Delta t : C_x,_z(t, t+\Delta t)=C_x,_z(0, \Delta t)\ $?
My basic intuition is that this is not the case, but I cannot find an appropriately pathological set of signals...Nor have I managed to prove transitivity.
Any help would be very much appreciated!
Thank you.

The simplest counter example I can think of (and it is certainly pathological) is $y(t)=0$ for any $x,z$ which are not jointly WSS.


Fast & accurate convolution algorithm (like FFT) for high dynamic range? It seems that FFT-based convolution suffers from limited floating-point resolution due to evaluating everything around the roots of unity, as you can see in the $10^{14}$-factor error in this Python code:

from scipy.signal import convolve, fftconvolve
a = [1.0, 1E-15]
b = [1.0, 1E-15]
convolve(a, b)     # [  1.00000000e+00,   2.00000000e-15,   1.00000000e-30]
fftconvolve(a, b)  # [  1.00000000e+00,   2.11022302e-15,   1.10223025e-16]

Are there any fast convolution algorithms that do not suffer from this problem?
Or is direct (quadratic-time) convolution the only way to get an accurate solution?

(Whether such small numbers are significant enough not to chop off is beside my point.)

Disclaimer:
I know this topic is older, but if one is looking for "fast accurate convolution high dynamic range" or similar this is one of the first of only a few decent results. I wanna share my insights I got on this topic so it might help somebody in the future.
I apologize if I might use the wrong terms in my answer, but everything I found on this topic is rather vague and lead to confusing even in this thread. I hope the reader will understand anyway.

Direct convolution is mostly accurate to machine precision for each point, i.e. the relative error is usually roughly or close to 1.e-16 for double precision for each point of the result. Each point has 16 correct digits. Rounding errors can significant for untypically large convolutions though, and strictly speaking one should be careful with cancellation and use something like Kahan summation and high enough precision data types, but in practice the error as almost always optimal.

The error of a FFT convolution apart from rounding errors is "global relative" error, meaning the error in each point depends on the machine precision and the peak value of the result. For example if the peak value of the result is 2.e9, then the absolute error in each point is $2\cdot10^9\cdot10^{-16} = 2\cdot10^{-7}$. So if a value in the result is supposed to be very small, let's say $10^{-9}$, the relative error in that point can be huge. FFT convolution is basically useless if you need small relative errors in the tail of your result, e.g. you have a somewhat exponential decay of your data and need accurate values in the tail. Interestingly if FFT convolution is not limited by that error, it has much smaller rounding errors compared to direct convolution, since you obviously do less additions/multiplications. This is actually why people often claim FFT convolution is more accurate, and they are almost right in some sense, so they can be quite adamant.

Unforunately there is no easy universal fix to get fast and accurate convolutions, but depending on your problem there might be one... I have found two:

If you have smooth kernels which can be approximated well by a polynomial in the tail, then the black-box Fast Multipole Method with Chebyshev interpolation might be interesting for you. If your kernel is "nice" this works actually perfectly: you get both linear (!) computational complexity and machine precision accuracy. If this fits your problem you should use it. It's not easy to implement however.

For some specific kernels (convex functions I think, usually from probability densities) you can use an "exponential shift" to get optimal error in some part of the tail of the result. There is a PHD thesis and a github with a python implementation using that systematically, and the author cals it accurate FFT convolution. In most cases this is not super useful however, since either it regresses back to direct convolution or you can use FFT convolution anyway. Although the code does it automatically, which is nice of course.

--------------------EDIT:--------------------

I looked a little bit at the Karatsuba algorithm (I actually made a small implementation), and to me it looks like it has usually similar error behavior like the FFT convolution, i.e. you get an error relative to the peak value of the result. Due to the divide and conquer nature of the algorithm some values in the tail of the result actually have better error, but I don't see an easy systematic way to tell which ones or in any case how to use this observation. Too bad, at first I thought Karatsuba might be something useful in-between direct and FFT convolution. But I don't see common use cases where Karatsuba should be prefered over the common two convolution algorithms.

And to add to the exponential shift I mentioned above: There are many cases where you can use it to improve the result of a convolution, but again it's not a universal fix. I actually use this together with FFT convolution to get pretty good results (in the general case for all inputs: at the worst same error as normal FFT convolution, at best relative error in each point to machine precision). But again, this only really works nicely for specific kernels and data, but for me both kernel and data or somewhat exponential in decay.


Impulse response of a second order LTI I have a set of measurement which I want to model with 2nd-order difference equations (first order eqs don't model well enough).

The equation is

$$y[n] = \alpha_1 y[n-1] + \alpha_2 y[n-2] + \beta_0 x[n] + \beta_1 x[n-1]$$

I want a generic formula for the impulse response.

I did a derivation, but being a total DSP beginner, I'm not sure whether it's correct or I'm missing some bits:

The transfer function is:

$$\begin{align}
H(z) &= \displaystyle \frac{Y(z)}{X(z)} \\
 &= \frac{\beta_0 + \beta_1 z^{-1}}{1 - \alpha_1 z^{-1} - \alpha_2 z^{-2}} \\
 &= \frac{\beta_0 z^2 + \beta_1 z}{z^2-\alpha_1 z -\alpha_2 } \\
 &= \frac{A_1 \, z}{z-p_1} + \frac{A_2 \, z}{z-p_2} \\
 &= A_1\frac{1}{1-p_1 z^{-1}} + A_2\frac{1}{1-p_2 z^{-1}} \\
\end{align}$$

where $p_{1,2}$ are the poles and $A_1$ and $A_2$ can be calculated using partial fraction expansion.

Therefore the impulse response is

$h[n] = \big(A_1 p_1^n + A_2 p_2^n \big) u[n], $

where $u[n]$ is the discrete step function.

Is this ok?

Your procedure is correct. But, you need to consider two things:

(i) Pole/zero locations.

(ii) Whether the system is causal or non-causal.

Your final answer is OK provided you consider that all the zeros and poles are distinct and the underlying LTI system is causal.

You will get two more solutions if the LTI system is considered non-causal, even if you assume to have distinct zeros and poles.

Hope this helps.


Initial state recognition in HMM I am building a speech recognition system using Hidden Markov Model in python. I referred to this and this question and its answers, which were very helpful.

In my approach, I split the continuous speech into separate words. I am thinking of using HMM to detect each word. So my states of HMM will be phones.

What I understood so far is that HMM estimates next state based on current state(phone). But I don't get how to estimate first state of HMM(i.e. the first phone of the word).

Can you suggest the best approach to use HMM to achieve this?

Also states of HMM will be phones, but I am not getting what can be observation in problem? There are multiple frame in a single phone and there is a feature vector corresponding to each frame. What should I use as observation?

The Baum-Welch algorithm uses the EM (Expectation Maximization) algorithm to estimate the model parameters $(T, E, \pi)$, where:

$T$: the transition probabilities
$E$: the emition probabilities
$\pi$: probability distribution on the states

Some years ago, I made the following quick-and-dirty implementation (may be fairly broken now), for the discrete case.
Hope this helps.


Get the transfer function when using the "designfilt" command in MATLAB To put things into perspective, I initially had these four graphs in MATLAB, whose data were sampled with a rate of 10 samples/sec :-

I wanted to remove some of the spikes on the blue graphs, so I thought about making a lowpass-filter that would take care of that. I used the command d = designfilt('lowpassfir','Filterorder', 5, 'CutoffFrequency',0.1,'SampleRate',10) . I then put the data from the blue graph through the filter, and plotted again. Now they look like this:-

As can be seen, some of the spikes have been removed and the graph have been smoothed out, which is what I desire. However, now I want to get the transfer function for the filter, but the only output I get in MATLAB is
d =

 digitalFilter with properties:

         Coefficients: [0.0287 0.1430 0.3283 0.3283 0.1430 0.0287]

   Specifications:
    FrequencyResponse: 'lowpass'
      ImpulseResponse: 'fir'
           SampleRate: 10
      CutoffFrequency: 0.1000
          FilterOrder: 5
         DesignMethod: 'window'

So my question is, can I just translate the coefficients directly into a transfer function of 5th order and with no zeroes? Like this:-
$$H(s) = \frac{1}{0.0287s^5 + 0.1430s^4 + 0.3283s^3 + 0.3283s^2 + 0.1430s + 0.0287} $$
If not, how can I get the transfer function from my designed filter?
Also, I didn't rigorously derive the arguments for the designfilt  command (filterorder, cutoff frequency) I just guessed at some appropriate values.

The filter is a discrete FIR filter, so the transfer function would be given as:
$$H(z) = .0287 + 0.1430z^{-1} + 0.3283z^{-2} + 0.3283z^{-3}+0.1430z^{-4}+0.0287z^{-5}$$
With the frequency response given by using the unit circle for the complex variable $z$, as in $z = e^{j\omega}$ for $\omega \in [0, 2\pi)$, with the sampling rate normalized to $\omega = 2\pi$.
You can use the freqz command in Matlab to plot or return the frequency response:
freqz([0.0287 0.1430 0.3283 0.3283 0.1430 0.0287])

Which creates a plot similar to what I have done using Octave below for these coefficients:

The frequency axis will be normalized to the radian angular sampling frequency, but you can refer to Matlab's documentation for other parameters that can be passed in to freqz to scale the frequency otherwise.


What is a low-pass filter for a sequence of tuples? I am currently reading SEGMENTATION AND RECOGNITION OF SYMBOLS WITHIN HANDWRITTEN MATHEMATICAL EXPRESSIONS. There is one sentence I don't understand:

  The preprocessing stage is subdivided into three steps:

  smoothing the input data by lowpass-filtering.

In this context your data usually looks like this:

[[(x1, y1), (x2, y2), ...],
 [(x12, y12), ...],
 ...
]

So you have list of coordinates that are ordered by time. Each sublist is exactly one line or a single point.

I've only heard of low-pass filters in the context of acoustics / electronics so far (I don't know much about these topics). In those contexts it's obvious what it means to let low frequencies pass / cut high frequencies.

So what is a low-pass filter in the context of handwritten symbol recogntion?

'Frequency' for a sequence of values is the measure of how sharp the values in that sequence change.

For example, the sequence $x_1(n)=[\begin{array}[cccccccc]&\cdots &1, &1, &1, &1, &\cdots\end{array}]$ has a frequency 0 but the  $x_2(n)=[\begin{array}[cccccccc] &\cdots &1 &1 &7 &1 &\cdots\end{array}]$ is a  higher frequency signal since there is a sharp change in sequence value (1 to 7 and 7 to 1).

A low pass filter is a system which attenuates the high frequency components. OR it reduces the sharpness (or introduces blur) in case of an image. OR a low pass filter does not allow sudden changes in values and smoothens these changes. Figure below shows low passing various sequences ($\style{color:blue}{blue:input}$ $\style{color:red}{red:output}$).

Consider a system whose output is a sequence $y(n)$ where n$^{th}$ value in $y(n)$ is the average of n$^{th}$, (n-1)$^{th}$ and (n+1)$^{th}$ value of input sequence is a simple low pass filter. If $x_2(n)$ is fed as input to this filter, then output, will  be $[\begin{array}[cccccccc] &\cdots &1 &3 &3 &3 &1 &\cdots\end{array}]$ and it is clear that the magnitude of change between consecutive values are reduced. or the high frequency components are attenuated.

Since you have a tuple, you have to separate it into two sequences ([x1,x2,....] and [y1,y2,....]). The low pass filtering should be done on these sequences separately.

In character recognition algorithms, low pass filtering is used in pre-processing to remove noises or unwanted small dots, having a width of one or two pixels. see an example of low pass filtering used for noise removal in image processing.


How to plot the periodic digital spectrum? In a previous answer, Dan Boschen wrote (emphasis mine):

  First let me explain the "unfolded" digital spectrum: If you allow the
  frequency axis of the sampled signal to extend to $\pm \infty$,
  instead of limiting to the unique digital frequency range of $\pm
> F_s/2$ (where $F_s$ is the sampling rate), you will see replicas of
  the original spectrum that is centered about 0 (DC) to also be
  similarly centered around every multiple of $F_s$. This is why we only
  need to show the spectrum from $\pm F_s/2$ (or even $0$ to $F_s/2$ for
  real signals) since this replicates everywhere else.  However, I find
  this visualization helps immensely in understanding many concepts in
  multi-rate signal processing as well as bridging analog and digital
  systems.

Indeed this was an(other) eye-opener to me.  Unfortunately, all the GUI tools I have at hand at this time only plot digital spectrum in the $0$, $F_s/2$ range--and I always thought it was the full spectrum. But as mentioned by Matt in a comment below "the spectrum is periodic".

Do you know how I can plot several periods of the unfolded digital spectrum using Matlab/Octave or Python Scipy?

As a simple example, consider the discrete-time signal

$$x[n]=a^nu[n],\qquad |a|<1\tag{1}$$

where $u[n]$ is the unit step function. The discrete-time Fourier transform (DTFT) of $x[n]$ is given by

$$X(e^{j\omega})=\frac{1}{1-ae^{-j\omega t}}\tag{2}$$

Note that $\omega$ is a normalized angular frequency:

$$\omega=2\pi f/f_s\tag{3}$$

where $f_s$ is the sampling frequency.

Also note from $(2)$ that $X(e^{j\omega})$ is a $2\pi$-periodic function in $\omega$, i.e., it is periodic with period $f_s$ in the frequency variable $f$. This periodicity of the spectrum is a general property of discrete-time signals.

Now you can plot the spectrum of $x[n]$ (i.e., $X(e^{j\omega})$) as a function of the frequency $f$, and depending on the chosen range of $f$ and the chosen sampling frequency, you will plot a certain number of periods. This is shown in the following Octave/Matlab script:

a = 0.9;
fs = 1;
fmin = 0;
fmax = 2.5;
f = linspace(fmin,fmax,1000);
X = 1 ./ (1 - a * exp(-1i*2*pi*f/fs) );

plot(f,abs(X))
xlabel('f')
title('|X(f)|')
grid on


Bounds of higher order derivative of bandlimited signals Here is a bandlimited function f(t) with bandwidth Ω:

The function f(t) is bounded in [-A,A].
Then the bound of the derivative of f(t) is bounded as:
|f'(t)|≤2πΩA.
So, what is the bound of its n-th order derivatives?
Is it as follows?

where n is the order of derivative.
If yes, then are the bounds of higher order derivatives much larger than the bound of f(t)?
I tried MATLAB simulations but the results were contrary.

What simulations did you try? I think what you mention is right.
If $\mathcal{F}[f(t)] = F(\omega)$, according to the inverse Fourier transform, we have
$$
f(t) = \frac{1}{2\pi}\int_{-\infty}^\infty F(\omega)e^{j\omega t}d\omega
$$
$$
f'(t) = \frac{d}{dt}\Big(\frac{1}{2\pi}\int_{-\infty}^\infty F(\omega)e^{j\omega t}d\omega\Big)= \frac{1}{2\pi}\int_{-\infty}^\infty j\omega F(\omega)e^{j\omega t}d\omega
$$
So the Fourier transform of derivative should be
$$\mathcal{F}[f'(t)] = j\omega F(\omega)$$
and the Fourier transform of $n$-th order derivative is
$$
\mathcal{F}[f^{(n)}(t)] = (j\omega)^nF(\omega)
$$


How to calculate time domain SNR using known sequence I am using the following formula to calculate SNR of a real world complex baseband signal sampled at 1x Nyquist.
SNR      = Rxy(tm)^2 / [ Px*Py - Rxy(tm)^2 ]
SNR (dB) = 10*log10(SNR)

where
Rxy(tm) = peak of the cross correlation at time delay, tm
Px = power in reference signal
Py = power in received signal

I verified proper implementation of the formula using simulated real-valued and complex-valued signals with and without noise.
On real data the SNR estimates using the above formula are too low (by 10.0+ dB). I manually verified the actual SNR a few different ways. I used spectral analysis to visually measure the signal power to the noise floor. I also measured the signal power to noise power (when signal is off), and both of those techniques give me an answer closer to what I expect.
I am flummoxed as to why this equation is not working on real-world signals. Do I need to take the DC bias (mean of data) into account and add that back to the SNR estimate? If I do that then I get values closer to what I expect.
Reference: Formula came from Principles of Communications (Tranter, Ziemer) textbook

The peak of the cross correlation should be the transmit signal power times the channel's attenuation.
Realizing that, it's really just
\begin{align}
\text{SNR} &= \frac{P_\text{signal}}{P_\text{noise}}\\
&= \frac{P_\text{signal}}{P_\text{received} - P_\text{signal}}\\
&= \frac{P_\text{tx}\cdot a_\text{channel}}{P_\text{received} - P_\text{tx}\cdot a_\text{channel}}\\
&= \frac{P_{crosscorr,max}}{P_\text{received}-P_{crosscorr,max}},
\end{align}
which is
=Rxy(tm)^2 / (Py-Rxy(tm)^2)

in your notation.
So, either that book is wrong or you're missing something there.


Modeling spring oscilliations with electric noise (50 hz) I have Analog to Digital converter from which I get data of measurements of 1kg load.
The aim is to model the simple spring-damper system, so that after entering correct initial data (x[0] = a, x'[0] = b) the model will give almost the same data as the original data from converter.

It is my first modeling attempt of a such real system. The problem is that the data from converter contains the 50Hz noise (electricity frequency) and if I use FIR filter I am afraid that the initial values needed for modeling the same wave data will be corrupted. How this process is done usually? As for modeling equation I use mx'' = -kx - c x'. And I have made numerical solution for this equation.

Edit 1: How about not filtering but modeling the system with taking into account the electrical noise?

The problem is that just taking the x[0] and x'[0] of your observed data and plugging them into the differential equation will be extremely sensitive to any kind of noise ; not only the 50Hz hum. You are trusting those two samples (rather than the entirety of your observed data) for your estimation, and you will get incorrect results from that, whether the 50Hz noise is present or not.

How to do it right? First, get rid of the 50Hz noise with a notch filter tuned to 50Hz and maybe to a few harmonics.

Then, what you want to do is fit the parameters $A$ and $\phi$ of a $x(n) = A e^{-\delta n} sin(\omega n + \phi)$ model to your data $y(n) = x(n) + noise$. You already know $\delta$ and $\omega$ since they depend on your known system parameter (damping and oscillation frequency). What you really want to know is the best estimate of $A$ and $\phi$ to reconstruct the "theoretical" part (in proper speak: signal subspace projection) of your observation.

Here is how you can estimate $A$ and $\phi$ without lame methods involving looking at peaks and zerocrossing:

You can write your signal model as $x(n) = Re (\alpha z^n) = \frac{1}{2}\alpha z^n - \frac{1}{2} \bar{\alpha} \bar{z}^n$ where $z$ is the known complex pole $z = e^{-\delta + j\omega}$, $\alpha$ is the complex amplitude $\alpha = Ae^{j \phi}$ and $\bar{z}$ denotes the conjugate.

Lay out a $N\times2$ matrix where N is your number of observation, the first column being $\frac{1}{2} [1, z, \ldots, z^N]$ and the second $-\frac{1}{2} [1, \bar{z}, \ldots, \bar{z}^N]$. Do a least square regression of that matrix (which is a basis of your signal subspace) onto the data matrix $[y(1), \ldots, y(N)]$. This will give you LS estimates of $\alpha$ and $\bar{\alpha}$ as the two solutions. Take the modulus and angle of the first to get $A$ and $\phi$. Plug this back into $A e^{-\delta n} sin(\omega n + \phi)$ to get the best fit theoretical model from your data.

# Estimate amplitude and phase of an exponentially damped model
# given observations y and system parameters delta and omega
n = numpy.arange(y.shape[0])
z = numpy.exp(-n * delta + n * omega * 1j)
z_conj = numpy.vstack((0.5 * z, -0.5 * numpy.conj(z))).T
s, _, _, _ = numpy.linalg.lstsq(z_conj, y)
print "Estimated amplitude: ", numpy.abs(s[0])
print "Estimated phase: ", numpy.angle(s[0])

If you want to also estimate $\delta$ and $\omega$ the procedure is a tad more complex and involves an EVD or SVD of the correlation matrix of your data to get a signal subspace basis, one more EVD to get the estimate of $z$, then the procedure sketched above to retrieve $A$ and $\phi$. You will have a robust estimate of all the parameters of the theoretical model from your data. Search for ESPRIT or MUSIC.

Note that this approach can be used to estimate the phase / frequency of a mixture of any number of damped sine waves. So for example, you could use this to retrieve the amplitude and phase of the 50Hz hum itself. Just put two additional columns in the matrix used for the least square procedure, with powers of $e^{j\omega_{noise}}$ and their conjugates. This will give the best fit amplitudes and phases of an exponentially damped sine wave (described by omega and delta) + 50 Hz hum model.


Why is carrier wave affected by percentage modulation? In amplitude modulation, Power of a carrier wave is represented by : $$P_c = \frac{A_c ^2}{2R}$$

Which is independent of modulation index!

Question: At what factor is the power of carrier wave dependent on!?

A conventional AM signal can be written as

$$x_{AM}(t)=A_c\big(1+\alpha\; m(t)\big)\cos(\omega_ct)\tag{1}$$

where $\alpha$ is the modulation index and $m(t)$ is the message signal. In Eq. $(1)$ it is assumed that the message signal $m(t)$ is normalized, i.e., $\max_t|m(t)|= 1$. The average power of $x_{AM}(t)$ is

$$\overline{x^2_{AM}(t)}=\frac{A_c^2}{2}\left(1+\alpha^2\overline{m^2(t)}\right)\tag{2}$$

If you assume now that the power $(2)$ is fixed to some value $P$, then the total available power must divided between the carrier term and the message-bearing term in $(2)$. Equating $(2)$ with the available power $P$, and expressing the carrier term $A_c^2/2$ in terms of the modulation index and the average power of the message signal gives

$$\frac{A_c^2}{2}=\frac{P}{1+\alpha^2\overline{m^2(t)}}\tag{3}$$

In the figure in your question it looks like they assumed $\overline{m^2(t)}=\frac12$, which for a maximum modulation index $\alpha=1$ results in

$$\frac{A_c^2}{2}=\frac23 P\tag{4}$$

This value corresponds to the right-most point of the curve labeled "carrier".


Invertibility of an ideal differentiator Is the system $y(t)= dx(t)/dt$  invertible or not?
If yes, please determine the inverse system for it.

It is not. It suffices to find a counter example. Let me tell you that "I am flat". Could you derive my actual value?

So, any constant function $f(t) =c$ is differentiable, and have the same derivative, $f'(t)=0$. Only from knowing that  $f'(t)=0$, you cannot recover the original (constant) function. Hence the system is not invertible.


Understanding Spatial Frequency of Images (using matlab) I wish to isolate a crosswalk from its surrounding, and I thought that as crosswalks have a repetitive pattern of thick white stripes, I might be able to isolate them by measuring their frequency range, and using a band pass filter with the right cut-off frequencies.

The problem is that I have no idea how to analyze the frequency of patterns in images. I know I should run fft2 on an isolated crosswalk, but I don't know what to do with the result.

If it is the repetitive pattern, you may observe some strong magnitude at some position in frequency domain. Use imagesc(abs(fftshift(fft2(img)))), you probably would see:

Those small white spots in the image correspond to the repetitive texture that is superimposed on the image. You may then use ordfilt2(F,roi^2,ones(roi)) to find out the roi * roi local maximum of the image. Here F = abs(fftshift(fft2(img))), roi is the region size of interest. Then change those values to zero:

With the real(ifft2(fftshift(F0))) you will recover the image in spatial domain. F0 is the frequency distribution after you implement the white spots removal on  fftshift(fft2(img)))


Why is the lowpass-filter coming right after the local oscillator when mixing an RF? In the book, instead of having a band- or lowpass-filter before mixing, it's coming right after mixing, which makes sense, considering there's gonna be a high-frequency component, which has to be filtered out.

But the lowpass-filter is called RF Image Rejection Filter, which irritates me, since that would be a band-pass instead and not something done after, but before mixing.

Curiously, afterwards there's an IF Band Selection Filter, which is meant to filter images?

Today I received an answer from the authors of the book. It turns out that figure 1.7a was wrong, and they sent me the updated one.

The corrected diagram will be added in this Summer's update!


Forward/reverse median filter. Reference needed I am using the double median filtering of my data. For this, I perform median filtering of the signal in forward direction, then in reverse direction (from end to start). Results are averaged. Example of the MATLAB code:

Signal=rnd(1,100);
Mfw=medfilt1(Signal,10);
Mrev=medfilt1(Signal(end:-1:1),10);
M=mean(Mfw,Mrev(end:-1:1));

I am happy with this approach, but now it is time to write a report...
I understand that I am not the first who use this kind of filtering.

Could you help me to find some references to this filter?

The problem is I even don't know the official name for this approach. "bidirectional median filter" shows me something else.

I've held off on answering this because I don't have MATLAB to test any of this.

I've employed the same trick with exponential smoothing.  The purpose of this is to have the lags in the opposite direction cancel each other out, so the resultant smoothed signal is centered instead of shifted from the source signal.  For sinusoids, the lag cancellation is exact.

So, I thought this might be the same motivation in this.  It does not seem to be the case.  (Side note:  In my searching I found the MATLAB function "movmedian".  I don't immediately see how it is that different from "medfilt1".)

I agree with Laurent Duval that if an odd number is used that the forward and backward filtered signals will be the same.  However, OP is using an even number.  It seems to me that the backward filtered signal will be just one index off of the forward one since the window size is the same and the median for each window should give the same result.  Thus, excepting the end cases, you should be able to just do the forward one and then define:

$$ M[x] = ( Mfw[x] + Mfw[x+1] ) / 2 $$

Which is a very simple FIR.

Two cents minus inflation worth.

Ced


Background subtraction for object tracking i am working on object tracking. i want to use background subtraction. i have tried to find out papers on background subtraction but i am really confused to use which one. I have already done tracking using just gray scale subtraction and thresholding. but now i want to model the background and then use it for tracking. i do not know much about background modelling.So, if anyone can suggest me some techniques and good papers on background subtraction and modelling, i will be greatly thankful.

The seminal paper on background modeling is Adaptive background mixture models for real-time tracking by Stauffer and Grimson.

If you have Matlab and the Computer Vision System Toolbox, then you can try using vision.ForegroundDetector, which implements a version of their algorithm.


What does it mean that the DFT equals the Complex Fourier Coefficients for even frequencies? I have a periodic signal, with period $1$

$$x(t) = \begin{cases}
1 \qquad & 0 \le t - \lfloor t \rfloor < \tfrac12 \\
0 \qquad & \tfrac12 \le t - \lfloor t \rfloor < 1 \\
\end{cases}$$

$\lfloor t \rfloor = \operatorname{floor}(t)$ is the floor() function, returning the largest integer no greater than the argument $t$.

$$ x(t+1) = x(t) \qquad \forall t \in \mathbb{R} $$

The complex Fourier series for $x(t)$ is

$$ x(t) = \sum\limits_{k=-\infty}^{\infty} c_k \ e^{i 2 \pi k t} $$

The complex Fourier coefficients are

$$\begin{align}
c_k &= \int_{-1/2}^{1/2} x(t) \ e^{-i2\pi k t} \ \mathrm{d}t \qquad \qquad k \in \mathbb{Z} \\
    &= \int_{0}^{1/2} 1 \ e^{-i2\pi k t} \ \mathrm{d}t \\
    &= \tfrac{1}{-i2\pi k } \big( e^{-i\pi k} - 1 \big) \\
    &= \tfrac{i}{2\pi k } \big( (-1)^k - 1 \big) \\
\end{align}$$

and are $0$ for even $k$.

Now, we sample $x(t)$ at $N\in 2\mathbb{N}$ time values,

$$\begin{align}
x[n] &= x(t_n) \\
     &= x\left(\tfrac{1}{N}n\right)
\end{align}$$

where $\tfrac{1}{N}$ is the sampling period, $N$ is the sampling frequency, and $t_n = \frac{n}{N}$, with $n=0, N-1$.

Let $X[k]$ denote the DFT of this finite sequence $x[n]$.

$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \ e^{-i2\pi nk/N} $$

One can show that $X[k] = 0$ for even $k$ as well, so $X[2k] = c_{2k}$, but $X[2k+1] \neq c_{2k+1}$.

I'm trying to come up with an intuitive explanation as to why $X[2k] = c_{2k}$, but $X[2k+1]\neq c_{2k+1}$. Obviously, we shouldn't really expect them to be equal in general since $X[k]$ is really just a Riemann sum approximation to $c_k$ with $N$ intervals, but in this case it seems there may be an explanation, since $X[2k] = c_{2k}$.

So far, all I can come up with is that since $x(t)$ isn't bandlimited, the DFT of the discrete sample of $x(t)$ is essentially trying to "fit" a bandlimited signal to the samples $x[n]$, and so for some reason this doesn't have any content at frequencies $2\pi (2k) = 4\pi k$ for any $k\in\mathbb{N}$.

Is there any specific reason for this? From the above argument, I have a feeling it has to do with aliasing, but I can't exactly make the connection.

okay, i am using these definitions for the continuous Fourier Transform

$$ \mathscr{F}\Big\{x(t)\Big\} \triangleq X(f) \triangleq \int_{-\infty}^{+\infty} x(t) \, e^{-i 2 \pi f t} \, \mathrm{d}t  $$

and inverse

$$ \mathscr{F}^{-1}\Big\{X(f)\Big\} \triangleq x(t) = \int_{-\infty}^{+\infty} X(f) \, e^{i 2 \pi f t} \, \mathrm{d}f $$

$f$ is "ordinary" frequency, not angular frequency (which for the continuous-time signals, the notation we use is $\Omega \triangleq 2 \pi f$).  this is a preferred representation of the Fourier Transform for electrical engineers because it is the easiest manner to apply the Duality property (lotsa similarity between the FT and the inverse).

in this answer is shown how sampling in one domain (e.g. the "time" domain) causes periodic extension (repetition and overlapping and adding) in the reciprocal domain (e.g. the "frequency" domain).  this is pointed out without regard to bandlimiting, so if the input is not sufficiently bandlimited, there will be aliasing.  given the assumption of a real input to start with (which means the Fourier transform or two-sided spectrum has symmetry about $f=0$) frequency components above Nyquist (which is $\tfrac12$ the sample rate) are folded back to locations below Nyquist and added to whatever frequency components existed there.

now, your continuous-time input

$$ x(t) = \sum\limits_{k=-\infty}^{\infty} c_k \ e^{i 2 \pi k t} $$

has a spectrum (which is not necessarily bandlimited, and is actually not bandlimited since the $c_k$ coefficients are not zero above some finite $k$) of

$$ X(f) = \sum\limits_{k=-\infty}^{\infty} c_k \ \delta(f - k) $$

first defining the Dirac comb:

$$ \operatorname{III}_T(t) \triangleq \sum_{n=-\infty}^{\infty} \delta(t-nT) $$

now, the ideally sampled signal (still represented in the continuous-time domain) is

$$\begin{align}
x_\text{s}(t) &\triangleq x(t) \cdot T \cdot \operatorname{III}_T(t) \\
 &= x(t) \cdot T \sum_{n=-\infty}^{\infty} \delta(t-nT) \\
 &= T \sum_{n=-\infty}^{\infty} x(t) \ \delta(t-nT) \\
 &= T \sum_{n=-\infty}^{\infty} x(nT) \ \delta(t-nT) \\
 &= T \sum_{n=-\infty}^{\infty} x[n] \ \delta(t-nT) \\
\end{align}$$

and also this

$$\begin{align}
x_\text{s}(t) &\triangleq x(t) \cdot T \cdot \operatorname{III}_T(t) \\
 &= x(t) \cdot T \cdot \sum_{m=-\infty}^{\infty} \frac1T e^{i 2 \pi m \frac{t}{T}} \\
 &= \sum_{m=-\infty}^{\infty} x(t) \ e^{i 2 \pi m \frac{t}{T}} \\
\end{align}$$

now from the latter, we know that the spectrum of the ideally sampling signal is this:

$$ \mathscr{F}\Big\{x_\text{s}(t)\Big\} = X_\text{s}(f) = \sum_{m=-\infty}^{\infty} X\left(f - \tfrac{m}{T} \right)
 $$

which is periodic in the frequency domain with period $\frac{1}{T}$ where your sampling period $T=\tfrac{1}{N}$.  that means your sample rate is $\tfrac{1}{T}=N$ and your Nyquist frequency is $\tfrac{N}{2}$.  that is the "folding frequency" or "foldover frequency".

Okay, winding up, you plug the expression of $X(f)$, the non-bandlimited Fourier series, into the equation for $X_\text{s}(f)$, which is the spectrum of the ideally-sampled signal.  And you plug in $N$ for $\tfrac1T$.

$$ X(f) = \sum\limits_{k=-\infty}^{\infty} c_k \ \delta(f - k) $$

$$\begin{align}
 X_\text{s}(f) &= \sum_{m=-\infty}^{\infty} X\left(f - \tfrac{m}{T} \right) \\
  &= \sum_{m=-\infty}^{\infty} X(f - N m) \\
  &= \sum_{m=-\infty}^{\infty} \sum\limits_{k=-\infty}^{\infty} c_k \ \delta\big((f - N m) - k\big) \\
  &= \sum\limits_{k=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} c_k \ \delta\big(f - (k+Nm) \big) \\
  &= \sum\limits_{k=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} c_{k-Nm} \ \delta(f - k) \\
  &= \sum\limits_{k=-\infty}^{\infty} \tilde{c}_k \ \delta(f - k) \\
\end{align}$$

where

$$ \tilde{c}_k \triangleq \sum_{m=-\infty}^{\infty} c_{k-Nm} $$

so the coefficients of the Fourier series of the sampled periodic signal $x_\text{s}(t)$ are different from those of $x(t)$.  even in the baseband where the DTFT (and following, the DFT) sees the spectrum.


Comparing spectra from hyperspectral images I extracted the mean spectra of regions of interests of different tissues from  hyperspectral images. How to compare those spectra by the difference of their reflectance and their spectral behavior  (I mean the spectra look similar but are shifted to the left or right).

Example of the mean spectra of two regions:

I would like to know the similarity between the two spectra. I tried cross correlation and Root mean square error but the results were not satisfying.

The goal is to predict a genetic distance of species from the spectra. For instance if the genetic distance of species 1 and species 2 is very height I would expect a different spectra. If the distance is 0 or very low,m the spectra are the same /very similar.

The goal is to predict a genetic distance of species from the spectra. For instance if the genetic distance of species 1 and species 2 is very height I would expect a different spectra. If the distance is 0 or very low,m the spectra are the same /very similar.

In that case, defining the distance between two spectra based on the spectra makes no sense – you want to define it based on the genetic distance.

So, you want to infer a model that describes genetic distance based on spectra. That's fine – but such a model will be very complex, and it's (as far as I can tell) not possible to find an easy function that will just give you what you want.

Unless you have a large dataset of genetic distances, I'll even go as far as saying that this is scientifically a very dangerous thing to do: your spectra are themselves large sets of data (a lot of points, say 1000). If you have, say, only 100 genetic distances, then you'll always find something that correlates well with the genetic distance within your large number of points, without having anything to do with it: the 101. genetic distance would not be represented by that model at all, i.e. it would only fit to the data set you had when you formed that model.

Basically, what you're describing is a problem that these days would usually be solved using neural networks. Again, this all depends on you having enough data to train your network, or else your network will not generalize to new spectral measurements and be useless outside your training data set.


What image-processing techniques are ideal for this particular shift/scale invariant template matching? A problem I had originally discussed here has evolved, and might have gotten a little simpler as I studied it in some more and attained new information.

Bottom line, I would like to be able to detect this pattern shown here, using computer-vision/image-processing techniques. As you can see, the ideal pattern is made up of four 'pings'. The object-recognition should be:

Shift Invariant

Horizontally, the image will be cyclical. (ie, Push to the right, comes out on the left, and vice-versa).
(Fortunately) Vertically, it is not cyclical. (ie,Push to either top or bottom and it will stop).

Scale Invariant (The pings can vary in 'thickness' as you can see.)

I could go on and on about it but I have attached images covering what I mean, please see below:

Of course, they can also be at a different 'scale', as can be seen from this family:

And finally, here are some 'realistic' scenarios of what I might actually receive, where there would be noise, the rows might 'fade' as you go towards the bottom, and of course, the  image will have plenty of spurious lines, artifacts, etc.

And of course, as a grand finale, there is the distinct possibility of this 'extreme' scenario:

So once again, I would like to ask for some guidance on what computer-machine-vision techniques I should utilize here in order to best detect the occurance of my pattern, where I need to be shift and scale invariant as you can see, and also be able to get decent results for the realistic scenarios. (Good news is that I DONT need it to be rotationaly invariant). Only idea I have been able to come up with thus far is a 2-D correlation.

I should add, that in reality, I do NOT have colored images - I will just get a martix of numbers, so I suppose we are talking about 'greyscale'.

Thanks in advance!

P.S. For what its worth I will probably be using open C V.

EDIT # 1:

Based on the comments, I have added the details you requested here:

For characteristics defining the data, we can assume the following:

Horizontal length of each ping can vary, but I know the upper and lower bounds of it. YES for anything within this bound, NO for anything outside. (Example, I know the length of the pings can be anywhere between 1 and 3 seconds, for example).
All pings need to be 'visible' for YES, however, last row can be missing, and still want to say 'YES'. Otherwise NO.
Vertical length, (the 'thickness') of each ping can vary, but again, also know upper and lower bound. (Similar to what you see in those images). YES for anything within that bound. NO for anything outside.
Height between each ping should always be the same for YES. If they are not, then NO. (Example, you can see how all the pings are of the same height relative to each other, ~110 on the vertical axis). So so 110 +/- 5 can be a YES, anything else has to be NO.

I think thats about it - but let me know what else I can add... (Also, everything shown here should register as a YES, btw).

Broadly speaking there are two primary approaches to solving this problem:

a. template matching or
b. matching with a parameterized model.

Personally, given the simplicity of the model, i would prefer the second approach for this problem.

Step 1: Identify the ping blobs

As a first step, extract the "Boxes" - essentially the yellow/blue squares. A simplest method here could be to just threshold the image. There doesn't seem to be much variations in the background except that as the noise increases it tends to be quite similar to the pings. Hence, the threshold can be global across the image - and i think you will be able to get a decent threshold that would work across the images. The threshold, however, should be smart so you can use something like Otsu's method. See this and this for more reference.

Improvement of the blocks

One of the good thing you can do to remove spurious points which looks similar to ping foreground where as some pixels inside the pings are also blue. You need here a morphological operation called "Opening". Here is one of the reference from HIPR. These type of operations requires smart shape that helps preserving similar shapes and removing others. In your case you can choose horizontal lines.  By the end of this, you have background and foreground pixels neatly separating the pings without touching each other.

Step 2: Mark/Register blocks

Now that you have got the cleanest possible blobs, time to mark them as 1,2,3,4 or A,B,C,D etc. This is usually accomplished by what is called some simple algorithm. Run through each pixel and if it is touching with any marked region, and is also foreground then assign it to that segment else create a new one. If you happen to find more than 4 important segments you need to choose the most relevant ones. Apply some heuristic here rather than any theoretical algorithm.

Step 3: Modeling of parameters

Now, once we get most likely top 4 candidates of the 4 corresponding pings, you can identify the bounding boxes - essentially the top left and bottom right corners (or left and right most pixels, as well as top and bottom rows of the blob). Find the bounding box, you can fined centroid the of each ping box.

Hence you can have four centroids $C_i \text{ where } i \in \{1,2,3,4\}$
Based on the description, the absolute location of $C_1$ can be arbitrary. However, following to $C_1$, the other centroids should be within certain range,

hence inter centroid distances can now be calculated:

$$

C_{1i}[x,y] = \{ C_1[x] - C_i[x],  C_1[y] - C_i[y] \}

$$

The other critical parameter is the length of the pings - which is $L_i \text{ where } i \in \{1,2,3,4\}$

Essentially you have 8 parameter vecoter:

$$
1.    \{ C_1[x],C_1[y] \}  \\
2..4  C_{1i}[x,y]  \\
5..8  L_i  \\
$$

Step 4: Classification
In the last step, now we need to have classification whether there is a ping-pattern yes or no.

For something like this, you can have a nice Bayesian classifier framework.

There are two Hypothesis you have

  H0 : The ping pattern exist.
  H1 : The ping pattern doesn't exist

Now, i am putting only 1 parameter $C_{1,i}[x]$ case here but you can extend this.

First of you study a lot of images where you know they belong to the image. Now, you can plot various histograms or apply some method of density estimation method. Read Pattern Classification by Duda Hart.

For simplicity you found that in cases where pings are present, the density function of parameter $C_{1,i}[x]$ is Gaussian with $ \mu = 20, \sigma = .1$ The sigma being small here could imply that there cannot be so much variations from the original.

This is actually conditional probability of $C_{1,i}$ given $H0$ is true. So you have values $P[ C_{1,i}[x] | H0 ]$

Now, from the test data

$P[ H0 | C_{1,i}[x] =x_1]$ is = given the reading of $C_1[x]$ being $x_1$, what is the probability of $H0$ ?

Hence, you can compute

$$

P[ H0 \space | C_{1,i}[x] ] = { P[ C_{1,i}[x] | H0 ] * P[H0] \over {P[ C_{1,i}[x]]  } }

$$

You need to combine this expression for the vector expression and need to fill in a lot of effort for exactly putting the parameter estimation to be done. But i am leaving you with basic approach here.

Please do your own math, this will be lengthy, but should be still quite intuitive.

Once, you combine expression, if $P[H0 | \text{all parameters}] > 1/2 $ then the ping pattern exists.

EDIT
Since you are already defining the specific criteria of 110 +/- 5 pixels, the last step can be simpler. You probably don't quite need to calculate all these probabilities if your classification criteria is fixed.


What framesize is usually used for audio effects? I want to implement an audio effect that reacts to an external parameter that varies over time. The book I have describes an effect in a function that processes a certain framesize of the audio input. But using 1024 samples or more might be too slow to react on the external parameter, since it would only change all 1024 samples. In the MATLAB simulation I am passing a pre-computed array to the function. Now I want to implement this on a DSP.
One solution I can think of is using a global variable and a multihreaded application. One thread writes the current paramenter to the global variable and the other thread processes the effect filter only with that current variable. That way I would not have to pass an array every time and always have the up-to-date parameter.
Another way would be calling the function every few samples. Since it is based on a second-order filter, the minimum I need is three samples. This would also be more than enough for the parameter.
Is there a standard way to go about in digital audio processing?

Typically this requires a trade off between CPU efficiency and latency. Frame based signal processing incurs at least one frame of latency and it you are planning to play live through the effect, you need to keep this quite low otherwise it get's distracting.
On the other hand, frame-based processing is way more efficient than sample by sample processing, since you can use DMA, vector instructions, highly optimized loops etc and you can amortize setup (coefficient calculation and loading, loading and storing of states, etc.) over the entire frame.
For live audio effects which moderate algorithmic complexity, something like 32 or 64 samples at 44.1 kHz or 48 kHz typically is a good starting point.


What is a "square window function"? When I read previous studies it mentioned a method of "square windowing function". I searched but I did not find anything talking about it. Can someone help to explain what is "square windowing function"? How to do it in matlab?

The description in the article is like this:

  A Cooley-Tukey FFT algorithm was implemented with a square
  windowing function and a 50% overlap of data segments to reduce
  spectral variance (Press et al., 1989).

This is from:

  Hamill et al., 1995. Shock attenuation and stride frequency during
  running. Human Movement Science. 14, 45-60.

A square window seems to be another name for a rectangular window.  A rectangular window is the default window used to apply an FFT to a subset of a longer stream of data and without any additional scaling or weighting of any of the subset's samples.


Same PSD from different FFTs Is it possible to obtain same PSD from different FFTs?
Can anyone explain it mathematically?

Sure.

Roughly speaking the PSD is simply the magnitude squared of the FFT. So any set FFT's that have the same magnitude but different phases will result in the same PSD. For example white noise and a delta impulse have the same PSD.

The following figure shows the impulse responses of a filter that has the same magnitude frequency response.


How to get the background noise of microphone? I am working on the automatic volume control for television which should take care of background noise. For this, I need to first find out the background noise. I have read in research papers that we need to subtract original audio signal which can be obtained from audio amplifier present in tv from the recorded signal by the microphone. Here microphone will be recording both original audio signal plus noise. But they have also mentioned about phase and amplitude correlation. I do not understand how to do that.

For a given TV, if you send a single short pulse (with an arbitrary amplitude of '1') out the speaker. The microphone will receive a slightly distorted pulse a short time later. (Typically you would send a pulse with enough volume that the microphone gets a clear reception of it, then reduce the microphone signal proportionally to simulate an original speaker signal with amplitude of '1'.) That distorted and delayed microphone signal is called the "impulse response" of that particular TV system. That impulse response (including the delay from the instant the signal was sent to the speaker) will need to be stored in your software for that TV.

Then to remove the distorted speaker signal from your microphone signal, you need to take the original electronic speaker signal and multiply it by that impulse response, and then subtract the result from your microphone signal. You do all this in the frequency domain because its far simpler and more efficient than doing it in the time domain. You don't even have to turn the final microphone signal back into time domain since you only wanted the ambient volume, anyway.

Look up "convolve impulse response".

EDIT:
To clarify (per your comment): The impulse response of the speaker-to-microphone path is a measure of the delay and distortion created by that path--and that's what you initially experimentally measure in the lab. But you don't know the impulse response of the noise-source-to-microphone path, so you can't correct for that distortion. All you can do is note the power spectrum of the noise (after subtracting the speaker contribution), and then adjust the speaker volume and maybe spectrum in order to help compensate for that noise. By the way, be sure to do the arithmetic (the multiply and the subtract) in the complex plane, because the FFT data is complex.


Why do most digital pll's only use the positive zero crossing? I've been trying to get my head around digital plls from knowing nothing and seen that most examples only use the positive zero crossing, not the negative. I've tried to find the answer but the best i could find is a textbook that said "complications arise when using the negative zero crossing so this type of PLL isn't as common" without further explanation.

The only reason that i can think of is that it would require more processing, but surely this would be worthwhile as you would be able to achieve faster lock times?

Getting accurate lock when using both edges or zero-crossings requires that the duty cycles of both the PLL local oscillator and the reference signal match and/or be as stable as the frequency.  If the duty cycles don't match when using both edges, then you would have to filter out some twice per cycle jitter which could well slow down lock speed.


Is there such a thing as a moving bandpass filter algorithm? I am looking for an algorithm, preferably available in a Python library, which passes only frequencies of interest, with the center frequency continuously changing as it moves through an array of audio response data.

I need this because I am analyzing the result of an audio sine wave sweep from 20Hz to 20000Hz and the recorded data contains noise.  Since I am driving the system with a known frequency and am only interested in the system's response near that frequency, I would like to have a tracking bandpass filter.

I have Googled terms such as "moving bandpass filter" "tracking bandpass filter" to no avail. I would appreciate any suggestions.

Depends on what exactly you want to do but easiest would be Heterodyne detector: multiply with a local oscillator (cos() and sin()) with variable frequency followed by a lowpass with a fixed bandwidth. The cutoff frequency of the lowpass is equivalent to the bandwidth of the equivalent bandpass. The frequency of the local oscillator determines the center frequency.

If you need non-uniform bandwidth (like 3rd octave for example), things get a little more complicated. You can do this primarily through frequency "warping".

Finally, you can always table and/or interpolate filter coefficients.


Intervals of convolution product Given are two Signals $x_1[n]$ and $x_2[n]$. $x_1[n]$ is in Intervall $[0,2]$ different than null and $x_2[n]$ is in Intervall $[0,3]$. The convolution product is Null outside intervals: $[0,6]$ and $[-2,10]$. Can someone explain me why? I understand for 0,6 it's simply multiplication but why $[-2,10]$?

An easy way to deduce the nonzero range of discrete-time convolutions is the following:
Consider convolution of $x[n]$ with $y[n]$: $$z[n] = x[n] \star y[n] .$$
Let the shorther sequence $x[n]$ has a nonzero range of $M_1 \leq n \leq M_2$ , while the nonzero range of longer $y[n]$ be $N_1 \leq n \leq N_2$ .
1- Expand $x[n]$ as weighted impulses: $x[n] = \sum_{k=M_1}^{M_2} x[k] \delta[n-k]$ .
2- Distribute the convolution over expanded $x[n]$:
$$ (\delta[n-M_1]+...+\delta[n-M_2]) \star y[n] = \delta[n-M_1] \star y[n] + ... + \delta[n-M_2] \star y[n] $$
3- Consider the first impulse $\delta[n-M_1]$, and the last impulse $\delta[n-M_2]$.
4- Apply shifting property for those two impulses :
$$ \delta[n-M_1] \star y[n] = y[n-M_1]  ~~~,~~~  \delta[n-M_1] \star y[n] = y[n-M_2] .$$
5- Determine the nonzero ranges of $y[n-M_1]$ and $y[n-M_2]$:
$y[n]$ is nonzero over $N_1 \leq n \leq N_2$, then $y[n-M_1]$ will be nonzero for $$ N_1 \leq n-M_1 \leq N_2  ~~~~~~\Rightarrow ~~~~~~  N_1 + M_1 \leq n \leq N_2 + M_1 .$$
Similarly for $y[n-M_2]$ the nonzero range will be: $$ N_1 \leq n-M_2 \leq N_2  ~~~~~~\Rightarrow ~~~~~~ N_1 + M_2 \leq n \leq N_2 + M_2 .$$
One can see that the output $z[n]$ has a nonzero range that begins at the leftmost sample of $y[n-M_1]$, and ends at the rightmost sample of $y[n-M_2]$. Yielding:
$$N_1+M_1 \leq n \leq N_2 + M_2 .$$
Length of $z[n]$ can be found as: $$L_z = L_x + L_y -1 = N_2+M_2 - N_1-M_1 + 1 .$$
Applying this to yourcase, $N_1 = 0$, $M_1 = 0$, $N_2 = 3$ and $M_2=2$ which yields the nonzero range to be
$$ 0 + 0 = 0 \leq n \leq 5 = 2 + 3 .$$


What is causing my ANC LMS update to diverge? I'm trying to implement a simulation of an ANC system with python, using this model here.
My simulation keeps diverging, and I honestly don't know why. I'm using a source for LMS adaptive filter from Mathworks here. When I comment out the LMS update function, the system is stable, and the output is exactly the input. But when I plug in the adaptive filter update, the system starts to diverge. So I thought the problem is with my LMS update. And I implemented the LMS adaptive filter with padasip package. The system is still diverging. Now I honestly don't know how I copied the model from Matlab wrong. Can someone help? Though It's not a syntax error, I pasted my code below for reference. I used all difference equations for every system.
Some values I used in this,
The input I'm using is a white noise file I generated from Matlab, 1dB power and sampling frequency of 16000 Hz.
b for S(z) is [0.5,0.5,-.3,-.3,-.2,-.2,]
b for S_est(z) is [0.466,0.533,-0.257,-0.274,-0.231,-0.175]
b for Main path S(z) is [0.0500,0,0.0200,0,-0.0000,0,-0.1250,0,-0.0500,0,0.0750,0,0.0300]
All I copied from the Simulink simulation of the same system
#ANC Simulation Main

import pyaudio, wave, struct, math
import numpy as np
from matplotlib import pyplot as plt
import padasip as pa

## Variables Setup
mu = 0.1
MAXVALUE = 2**15-1  # Maximum allowed output signal value (because WIDTH = 2)

# Initialization of adaptive weight w
w = np.zeros(13)

# Main Path Filters
order_path = 12
b_path = np.array([0.0500,0,0.0200,0,-0.0000,0,-0.1250,0,-0.0500,0,0.0750,0,0.0300])
x = np.zeros(13)
y_main_path =0

#Adapt filter
order_adapt = 12
y_adapt =np.zeros(13)

filt = pa.filters.FilterLMS(13, mu=mu)

#second path
order_sec_path = 12
b_sec_path = [0.5,0.5,-.3,-.3,-.2,-.2,0,0,0,0,0,0,0]

y_secpath = 0

#est sec path
order_est_sec_path = 12
b_est_sec_path = np.array([0.466,0.533,-0.257,-0.274,-0.231,-0.175,0,0,0,0,0,0,0])

est_sec_out = 0
y_estsecpath = np.zeros(13)

# File names
wavfile = 'matlab_1db_wn.wav'
output_wavfile = 'ANC_Result.wav'

## Read WAV file
wf = wave.open(wavfile,'rb')

CHANNELS        = wf.getnchannels()     # Number of channels
RATE            = wf.getframerate()     # Sampling rate (frames/second)
signal_length   = wf.getnframes()       # Signal length
WIDTH           = wf.getsampwidth()     # Number of bytes per sample

print('The file has %d channel(s).'            % CHANNELS)
print('The frame rate is %d frames/second.'    % RATE)
print('The file has %d frames.'                % signal_length)
print('There are %d bytes per sample.'         % WIDTH)

# Read first BLOCKLEN
binary_data = wf.readframes(1)

## Output WAV file
output_wf = wave.open(output_wavfile, 'w')
output_wf.setframerate(RATE)
output_wf.setsampwidth(WIDTH)
output_wf.setnchannels(CHANNELS)

## Open audio stream
p = pyaudio.PyAudio()
stream = p.open(
   format      = p.get_format_from_width(WIDTH),
   channels    = CHANNELS,
   rate        = RATE,
   input       = False,
   output      = True )

## Main Loop
while len(binary_data) > 0:

   # convert binary data to numbers
   input_block = struct.unpack('h', binary_data)
   input_value = input_block[0]

   x = np.delete(x,-1)
   x = np.insert(x,0,input_value)

   y_main_path = np.dot(b_path,np.transpose(x))

   y_adapt = filt.predict(input_value)

   #y_adapt = np.delete(y_adapt,-1)
   #y_adapt = np.insert(y_adapt,0,adapt_out)

   est_sec_out = np.dot(b_est_sec_path,np.transpose(x))
   y_estsecpath = np.delete(y_estsecpath,-1)
   y_estsecpath = np.insert(y_estsecpath,0,est_sec_out)

   #sec path
   y_secpath = np.dot(b_sec_path,np.transpose(y_adapt))

   output = y_main_path - input_value

   # LMS update
   filt.adapt(y_main_path,input_value)

   output = np.clip(output, -MAXVALUE, MAXVALUE)
   output = output.astype(int)
   binary_date = struct.pack('h', output)
   stream.write(binary_data)
   output_wf.writeframes(binary_data)
   binary_data = wf.readframes(1)

stream.stop_stream()
stream.close()
p.terminate()
# Close wavefiles
wf.close()
output_wf.close()

Turns out it's just a scaling issue. Because the struck.unpack() gives out a signed 16-bit value, the step size is too large for it. Two ways to go around this, you scale down your input and scale back up when outputting, or you just adjust the step size mu.
I used a mu value of 10e-8 and seems like this value works differently in different environments, reasons unknown. With jupyter-notebook is around 10e-7, python idle is 10e-8, and Matlab is 10e-9. Just tweak around that range with the code I posted above and should give you a converged result.
Hope this helps whoever is reading.
edit; for MatLab, I scaled the input by 32768 so the mu value is 10e-9, if max amplitude is 1, use 0.1.


Kaiser window approximation Are there any "good" low-computational-cost approximations for parameterized Kaiser window generation suitable for small systems/languages that do not include Bessel functions in their standard library?

(maybe only to 8 to 12-bit accuracy, and something that can be run in Basic on an Apple II, et.al.)

just implement the Modified Bessel function.  it's easy.

i always like my window definitions centered about zero, since pretty much all of them are even symmetry.

i'll do this in discrete-time, but it's essentially the same thing in continuous-time.

Kaiser window:

$$ w[n] \triangleq \begin{cases} \frac{1}{I_0(\beta)} I_0\left(\beta \sqrt{1 - \left(\frac{n}{M/2}\right)^2} \right) \quad \quad & |n| \le M/2 \\
0 & |n|>M/2 \\
\end{cases} $$

$I_0(x)$ is the 0th-order modified Bessel function of the 1st kind. $M+1$ is the number of non-zero samples or FIR taps (the FIR filter order is $M$ and, in my centered and symmetrical case, must be even). $\beta$ is a "shape parameter" and O&S recommend this heuristic:

$$ \beta = \begin{cases} 0.1102 \cdot (A-8.7) & A>50 \\
0.5842 \cdot (A-21)^{2/5} + 0.07886 \cdot (A-21) \quad & 21 \le A \le 50 \\
0.0 & A<21 \\
\end{cases}$$

$$ M = 2 \left\lceil \frac{A-8}{4.57 \cdot \Delta\omega} \right\rceil $$

$A$ is the desired stopband attention in dB and $\Delta\omega$ is the desired width of the transition band in normalized angular frequency.

finally, the Bessel is evaluated as:

$$ I_0(x) = 1 \ + \ \lim_{K \to \infty} \ \sum\limits_{k=1}^{K} \left(\frac{x^2}{4}\right)^{k} (k!)^{-2} $$

when you evaluate this with a computer, pick a $K$ decently large (my guess is that $K=32$ is good enough) and evaluate the summation starting with $k=K$ and work it backwards to $k=1$ to keep numerical accuracy.  you might want to use Horner's method.

$$ I_0(x) \approx  1 + x^2\left( \tfrac{1}{(1!)^2 \, 4^1} + x^2\left(\tfrac{1}{(2!)^2 \, 4^2} + x^2\left(... + \, x^2\left(\tfrac{1}{((K-1)!)^2 \, 4^{K-1}} +  x^2 \tfrac{1}{(K!)^2 \, 4^K} \right) \right) \right) \right) $$

you can evaluate all of the $(k!)^{-2}$ in advance with a short table.

Alright, someone made me do some work.  This took about 45 minutes to code up and debug.  So here is my MATLAB code for implementing the 0th-order Bessel function of the first find (which is $I_0(x)$ above):

function y = mybessel(x)
%
%   Computes the 0th-order Modified Bessel function of the first kind
%

    K = 32;

    bessel_coef =  zeros(1,K);

    kfac = 1;
    two_to_the_k = 1;

    for k = 1:K
        kfac = kfac * k;
        two_to_the_k = two_to_the_k * 2;
        bessel_coef(k) = 1/(kfac*two_to_the_k)^2;   % compute power series coefficients in advance
    end

    x = x.^2;

        y = x .* bessel_coef(K);
    for k = K-1:-1:1
        y = x .* (bessel_coef(k) + y);              % Horner's method
    end

    y = 1 + y;

end

and here is the test code:

x = linspace(-16, 16, 32*4096+1);

I_0 = besseli(0, x);            % MATLAB's modified bessel

y = mybessel(x);                % my bessel

figure(1)
plot(x, I_0)
hold on
plot(x, y)
hold off

figure(2)
plot(x, y - I_0)                % plot error

with results:

and error:

for $|x| \le 16$, the relative error is less than $10^{-15}$.  the error increases with increasing $|x|$.


Intuitive explanation of cross-correlation in frequency domain According to the cross-correlation theorem : the cross-correlation between two signals is equal to the product of fourier transform of one signal multiplied by complex conjugate of fourier transform of another signal. After doing this, when we take the ifft of the product signal, we get a peak which indicates the shift between two signals.

I am not able to understand how this works? Why would i get a peak which indicates the shift between two signals. I got the math from : http://mathworld.wolfram.com/Cross-CorrelationTheorem.html
but i am not able to understand as to what this means intuitatively. Can somebody please provide some explanation or point me to the right documents?

Thanks!

The concept is based on the convolution theorem, which states that for two signals $x(t)$ and $y(t)$, the product of their Fourier transforms $X(f)$ and $Y(f)$ is equal to the Fourier transform of the convolution of the two signals. That is:

$$
\mathcal{F}\{x(t) * y(t)\} = \mathcal{F}\{x(t)\}\mathcal{F}\{y(t)\}
$$

You can read more on the derivation of this theorem at the above Wikipedia link. Now, convolution is a very important operation for linear systems in itself, so the theory on its properties is well-developed.

However, what you're looking for is the cross-correlation between $x(t)$ and $y(t)$. Here's the key: you the cross-correlation integral is equivalent to the convolution integral if one of the input signals is conjugated and time-reversed. This allows you to utilize theory developed for evaluating convolutions (like frequency-domain techniques for calculating them quickly) and apply them to correlations.

In your example, you're calculating the following:

$$
\mathcal{F}\{x(t)\}\left(\mathcal{F}\{y(t)\}\right)^*
$$

Recall that in the Fourier domain, complex conjugation is equivalent to time reversal in the time domain (this follows directly from the definition of the Fourier transform). Therefore, using the first equation given above, we can state that:

$$
\mathcal{F}\{x(t) * y^*(-t)\} = \mathcal{F}\{x(t)\}\left(\mathcal{F}\{y(t)\}\right)^*
$$

If you then take the inverse Fourier transform of this equation, the signal you're left with is the cross-correlation between $x(t)$ and $y(t)$.


Applications of Correlation in Signal Processing There have been discussions on the differences between convolution and correlation, autocorrelation, cross-correlation, et al.

In spite of looking at the visualizations, and understanding the formula for calculation of convolution and correlation, I don't understand the need for correlation - is it required for identifying a signal uniquely?
What is the distinction between correlation and cross-correlation (if any exists)?

I get that convolution can be important in determining system responses to any arbitrary input signals when the impulse response is known, but what are the practical scenarios where correlation is helpful?

I know the theoretical differences between the concepts and this link isn't helpful for my question.

A relatively simple application where correlation is used is in determining timing of a communications signal.  A known synchronization signal will be sent periodically so that the receiver can use it as a point of reference.  How does the receiver use it as a point of reference?  It correlates the incoming signal stream against the known sequence, and when the correlation peak is above a certain threshold, it will use the peak to establish symbol timing and make decisions on the symbols to extract a binary sequence.

Correlation usually means cross-correlation, but I think the latter term only exists to be distinguished from auto-correlation.


Conversion into the Q31 Format I have a Setup, where I receive 24-bit ADC samples in the two's complement in 8-bit chunks. Usually I just "converted" them into a signed 32-bit integer like this:

sample=((byte1<<24)|(byte2<<16)|(byte3<<8))>>8;

Now, for educational purposes I moved the project to an Arduino Due and want to use the CMSIS-DSP Library of the Cortex-M3.

The DSP Methods Support the Q1.31 Format. I've read a lot of different sources but don't quite understand what this means for me.

As I understood, the Q1.31 format can have values between -1 and 1. Am I supposed to scale my samples somehow? Or is it kind of implicit scaling by just using it with my samples?

Furthermore the documentation provides the following information:

RFFT Size: 512
Input Format: Q1.31
Output Format: Q9.23
Number of bits to upscale: 8

I don't understand what this means for the output I would get.

Can anyone point out to me what this means?
I'm quite new at this :-)

This is actually pretty complicated (if you haven't done it before) and there isn't enough space here to do go through it thoroughly. Start here https://en.wikipedia.org/wiki/Q_(number_format)

Quick (simplified) overview: We use 8 bit unsigned numbers as an example. Signed makes it more complicated.

When you multiply two integers you basically get twice the number of bits. The result of an 8 bit number multiplied with another 8 bit number requires 16 bit for accurate representation. So after a few multiples you end up with hundreds of bits which simply doesn't work. So you need to decide which bits to keep and which ones to throw away.

The Q notation is a bookkeeping process to help with this. You assume that somewhere there is a decimal point. "Q" denotes where the decimal point is. Q=0 means all the way to the right, Q=8 means all the way to left. To convert from integer to floating point you simply multiply by $2^{-Q}$.

The largest number you can represent this way is 255 for Q0, $\frac{255}{256} = 0.966$ for Q8, and $\frac{255}{16} = 15.9375$ for Q4. The smallest non zero number you can represent is 1 for Q0,  $\frac{1}{256} = 0.004$ for Q8 and   $\frac{1}{16} = 0.0625$ for Q4. So you basically trading the absolute range against resolution.

Let's say you want to represent the number 2.5. You can't do this in Q0 (you could chose 2 or 3 but there is nothing in between) and you can't do it in Q8 either since it's larger than 1. Q1, however would work. 0b00000101 in Q1 = 2.5 Q1. Q3 would work as well: 0b00010100 in Q3 = 2.5. So any floating point number can be represented as different bit pattern, depending on what Q you chose.

If you want the same floating point number in a different Q you need to shift. shifting. Shifting to the left increases Q, shifting to the right decreases Q.

In Q8 the largest number you can represent is smaller than 1. That means it's impossible that any multiplication would give a result that's larger than any of its inputs and hence it can never overflow. To produce an 8 bit result you keep the MSBs of the 16 bit result and throw away the LSBs.

For a Q0 multiply you keep the LSB and throw away the MSBs. This can and will overflow. 2*2 is 4 but 16*200 will overflow since the result cannot be represented in 8 bits.

If you multiply two numbers, both the number of bits and the Qs will add. So if you multiply an 8-bit Q6 with an 8-bit Q5, you get a 16-bit Q11. To convert 16 bit back to 8 bit you have to throw some bits away. Throwing away LSBs (least significant bits) reduces the number of bits and the Q by the same amount. In our example, removing the lower 8 bits and keeping the upper 8 bits would result in an 8-bit Q3.

So what does this all mean in practice? It basically determines how you implement multiplies and which bits you keep. In general a multiply looks like this

Multiply 2 N bit numbers into a 2N wide accumulator or register
Shift the result depending on the desired output $Q_C$. $S = Q_C - (Q_A + Q_B - N)$.
Keep the upper N bits and transfer to result variable

The tricky part here is to figure which Q to chose for each variable. If the Q is too low, you can overflow and if it's too high you can get very bad signal to noise ratio. Often this requires careful analysis of the specific algorithm and the properties of the input signals.


What is the difference between pulse shaping and equalizers when used to cancel Inter-symbol Interference (ISI)? According to Nyquist, to transmit R symbols/sec without ISI, the minimum bandwidth required is R/2 Hz, and the sinc pulse shape satisfies this condition.

If we use equalizers instead of pulse shaping, do we need to keep this condition ( minimum R/2 Hz to transmit R symbols/sec)?
If we use equalizers, should the total system transfer function be of the sinc function shape, thus equivalently achieve the Nyquist condition?

1.
This is a question investigated by many researchers decades ago. They discovered that the bandwidth limitation of R/2 Hz is not a fundamental limit set by nature. It is just a criterion if we don't want to do any further processing at the Rx (that also doesn't mean that we can infinitely break this limit).

As long as we are willing to pay the cost in terms of Rx complexity, we can go higher than the rate R within the same bandwidth R/2. Mazo found that 25% faster than the rate R can be accommodated without any penalty other than the higher computational complexity. Pulse shapes are still the same but they are packed more closely.

It was mostly an academic thing before the advent of turbo codes. After the turbo codes, iterative receivers became a reality and then people started implementing turbo equalization where the decoder and the equalizer exchange extrinsic information to improve their belief of each bit.

Then, incorporating a higher symbol rate became possible and this is known as Faster than Nyquist (FTN) signaling. The idea was further enhanced to incorporate frequency domain packing as well (think OFDM) and now people break the limits in both dimensions and still recover the information.

Coming back to your question, no we don't strictly need to obey this limit as long as we are aware of the price we have to pay.

2.
I assume that you are asking about the system cumulative impulse response (and not the transfer function). When employing trellises, Viterbi algorithms and iterative processing, the concept of a pulse shape does not go along that far. The samples affected with inter-symbol interference are taken and this interference is canceld out in some fashion, e.g., sequence estimation or iterative receivers.


From transfer function to differential equation I have the below detailed solution (boxed in blue) that I don't understand completely:

I can reconstitute the differential equation from:
$$ (1+Ts) X(s) = K_v U(s) $$
$$ x(t) + T\dot x(t) = K_v u(t) $$
but I don't get how do they get $x(t)$ from the initial conditions?
I realize that my question is more about differential equations but it is something I obviously missed before...

Given a transfer function
$$G_v(s) = \frac{k_v}{1 + sT} \tag{1}$$
the corresponding LCCDE, with $y(t)$ being the solution, and $x(t)$ being the input, will be
$$ T ~\dot{y}(t) + y(t) = k_v ~x(t) \tag{2} $$
Your formulation replaces $x(t)$ with a unit-step $u(t)$, and $y(t)$ with $x(t)$, yielding
$$ T ~\dot{x}(t) + x(t) = k_v ~u(t) \tag{3} $$
or equivalently
$$ \dot{x}(t) + \frac{1}{T} x(t) = \frac{k_v}{T} ~u(t) \tag{4} $$
Eq.4 represents a typical first order, constant coefficient, linear, ordinary differential equation (abbr LCCDE) whose solution procedure is as follows:
First, find the homogeneous solution to the Eq.4 with RHS being zero, as
$$ x_h(t) = A e^{-t/T} \tag{5} $$
where $A$ is an unknown coefficient to be determined by the initial condition $x(0)$.
Second, find the particular solution of Eq.4 with the RHS being $u(t)$. For this you can use the method of undetermined coefficients, by assuming a solution of the type
$$ x_p(t) = B u(t) \tag{6}$$
an letting $t \to \infty$, will yield in Eq.4 as
$$ 0 + \frac{B}{T} = \frac{k_v}{T} \tag{7} $$
where you find $B$ as
$$ B = k_v \tag{8}$$
Then you find the complete solution by adding the homogeneous part to the particular part, and solving for the unknown constant $A$ from the initial condition $x(0)$
$$ x(t) = x_h(t) + x_p(t) = A e^{-t/T} + k_v ~~~~,~~~~ t >  0\tag{9} $$
The initial condition $x(0)$ is given as $0$ in the problem, hence
$$ x(0) = 0 = A + k_v \implies   A = -k_v \tag{10} $$
And the overall solution of the initial value problem in Eq.4 with $x(0)=0$ becomes:
$$ x(t) = -k_v e^{-t/T} + k_v = k_v (1 - e^{-t/T}) u(t) \tag{11} $$
Note that we look for a causal solution, unless otherwise stated, hence assume $ t>0$.


How to convert NIST to RIFF wav file I have a dataset of audios(namely TIDIGITS, collected at Texas Instruments, Inc.) and they are in .wav format. However, when I try to read them from python, I get the following error:
Error: file does not start with RIFF id

Error is thrown from the following line:
import wave
...
f = wave.open(/path/to/audio/example.wav, "rb")

I searched on the Internet, and it looks like the problem is that "wave" library does not support NIST formats and the audio are in the NIST SPHERE format, according to the dabaset's website:
https://catalog.ldc.upenn.edu/ldc93s10
So I transformed them into .wav using "sox" program:
find . -name '*.wav' | xargs -I@ bash -c "sox @ @"

The command reads though all of my ".wav" files and then applies
sox name.wav name.wav

to transform them to RIFF format.
However, now when I play the files in my media player, I get no sound.
The waveform of the audio becomes as follows:

Where it should be something like this:

I think the problem lies at the sox command.
Can anyone help me to convert all of my data to a RIFF format, so that I could read them from the "wave" library?
EDIT: Also, now I noticed, that when I apply the sox command to only one file, the file is converted to "RIFF" format and the waveform is similar to the second image, however when I apply them to all of them with the above described command, the audios loose their sound.
Thank you in advance!

OK, I solve the problem. The problem was that I was overriding the files when I used the "sox" command, as I want them to convert into RIFF format. However, "sox" didn't convert them correctly when using the same file name, i.e. when overriding the files. When I changed the file name of the overridden files, it converted them into RIFF format.


Pulse shaping and Baseband filtering I'm slightly confused about baseband pulse shaping.
Let's assume I have a complex data vector in an arbitrary complex constellation (QAM for example). I would like to pass this complex vector through a pulse shaping filter (RRC or RC).

As far as I know RC/RRC filters have a real impulse response, thus I can generate a sampled version of the filter.
Am I correct to assume the filter is real ? Do I have to represent the filter as an analytic function (ie: $h_{RC}(t)=h_{RC}+j\hat{h}_{RC}$ ) before convolving my vector with it ?

Your complex baseband signal before modulation is given by

$$s(t)=\sum_{m=-\infty}^{\infty}a_mh(t-mT)\tag{1}$$

where $a_m$ are the complex symbols, $h(t)$ is the impulse response of the transmit filter, and $T$ is the symbol period. As you correctly assumed, $h(t)$ is usually real-valued, so you need to filter the real and imaginary parts of your symbols with the same filter $h(t)$ to get the complex baseband signal $(1)$:

$$s(t)=\sum_{m=-\infty}^{\infty}\text{Re}\{a_m\}h(t-mT)+j\sum_{m=-\infty}^{\infty}\text{Im}\{a_m\}h(t-mT)\tag{2}$$

There is no need to modify the real-valued impulse response $h(t)$.


Interpretation of system response to unit impulse I want to understand the meaning of unit impulse response when time flows. An example is given below:

In the example the left signal depicts a $\delta[n]$ function value of which is $1$ only at $n=0$. The LHS signal is response of the system to the unit impulse.

I want to understand how to interpret this when time flows, i.e. at $n=-4,\ n=-3$, what is response of the system?
Or at $n=0,\ n=5$, etc. what is response of the system?

In the figure that you have posted, the right plot shows the impulse response $h[n]$ according to the convention that only its nonzero samples are displayed. And therefore, it's understood that all the remaning samples of $h[n]$ are zero unless otherwise stated explicitly.

This means that $h[0]=0$, $h[1]=0$,...,$h[12]=0$. Which can be stated in the following:

$$ h[n] = \begin{cases}
0 &,& \text{ for } n \le 4 \\
\alpha_n &,& \text{ for } 5 \le n \le 11 \\
0 &,& \text{ for } n \ge 12 \\
\end{cases}
$$

In this setting the domain of support for the finite length sequence $h[n]$ is stated as $n \in [5,11]$ , which means that $h[n]$ is certainly zero for $n$ out of this interval, while it can have nonzero as well as zero values inside this interval.


Convolution of $f(2x)$ and $g(3x)$ As I know, convolution is defined as $f(x)*g(x) = \int_{-\infty}^{+\infty}f(\tau)g(x-\tau)d_{\tau}$, but what if we want to convolve $f(2x)$ and $g(3x)$? It should be like $f(2x)*g(3x) = \int_{-\infty}^{+\infty}f(2\tau)g(3x-\tau)d_{\tau}$ or $f(2x)*g(3x) = \int_{-\infty}^{+\infty}f(2\tau)g(3x-3\tau)d_{\tau}$ or anything else?

Replace $x$ by $x - \tau$ . So option b i.e $3x-3\tau$.


Confusion about adjoint of convolution I have the following problem with convolution operators. I know that if I have a filter $g$, then the adjoint filter would be given by $h(x,y)=g(-x,-y)$, so that $\langle g * u,v \rangle = \langle  u,f *v\rangle $ holds for the $\mathbb{R}^2$ scalar product.

I also know that $Ff(-\omega) = \overline{Ff(\omega)}$ for real $f$. So it seems I should hold that

$ h(x,y) = F^{-1}\overline{Fg} = F^{-1}Fg(-x,-y) = g(-x,-y)$

However when I compare

g = ifft2(fft2(g))

( for some gaussian kernel, e.g matlab: g = padarray(fspecial('gaussian',11,2),[5,5]) )

and

h = ifft2(conj(fft2(g))

I get the expected result for g

But when I compute h, which should be equal to g because g is gaussian

the result is shifted by (1,1)

which is somewhat unexpected.
I wonder if there is some subtlety in the discrete <-> continuous conversion that I have missed ?

Complex conjugation in the discrete Fourier domain should result in spatial domain horizontal and vertical mirroring around the origin. The origin happens to be at row 1 column 1, as is evident from MATLAB's definition of FFT, with Y = fft(X):

$$Y(k)= \sum^n_{j=1} X(j) e^{−2πi(j−1)(k−1)/n}$$

For $j = 1$ the phase is zero for all $k$ so anything on column 1 or on row 1 will stay there through complex conjugation. Column 11 and row 11 are each 10 steps from the origin horizontally or vertically, as are column 12 and row 12 if you measure the distance across the border. So the reversal swaps the data between columns 11 and 12, and between rows 11 and 12. Figure 1 below illustrates the mirroring.

Figure 1. Row or column indices and the corresponding phases of frequency $k = 2$ plotted as angles on the circle.

If the dimensions of your data were even rather than odd as they are now, there would be a middle row and a middle column that would be at equal distance from the origin both directly and across the border.


Derive Frequency Representation of Impulse Train Function I want to walk through the derivation of the frequency representation of an impulse train.

The definition of the impulse train function with period $T$ and the frequency representation with sampling frequency $\Omega_s = 2\pi/T$ that I would like to derive is:

\begin{align*}
  s(t) &= \sum\limits_{n=-\infty}^{\infty} \delta(t - nT) \\
  S(j\Omega) &= \frac{2\pi}{T} \sum\limits_{k=-\infty}^{\infty} \delta(\Omega - k\Omega_s) \\
\end{align*}

Using the exponential Fourier series representation of the impulse function and applying the Fourier transform from there results in:

\begin{align*}
  s(t) &= \frac{1}{T} \sum\limits_{n=-\infty}^{\infty} e^{-jn\Omega_s t} \\
  S(j\Omega) &= \int_{-\infty}^\infty s(t) e^{-j\Omega t} dt \\
  S(j\Omega) &= \int_{-\infty}^\infty \frac{1}{T} \sum\limits_{n=-\infty}^{\infty} e^{-jn\Omega_s t} e^{-j\Omega t} dt \\
  S(j\Omega) &= \frac{1}{T} \int_{-\infty}^\infty \sum\limits_{k=-\infty}^{\infty} e^{-j(k\Omega_s + \Omega) t} dt \\
\end{align*}

To get from there to the end result, it would seem that the integration would need to be over a period of $2\pi$. Where $\Omega = -k\Omega_s$, the exponent would be $e^0$ and integrate to $2\pi$ and for other values of $\Omega$, there would be a full sine wave that would integrate to zero. However, the limits of integration are negative infinity to positive infinity. Can someone explain this? Thanks!

@MattL suggested a nice, simple way to see the above result.

But if you want to see the result in the normal analysis equations you mentioned, you can do like below.

Say S(t) is a periodic train of impulses.So S(t) can be written as

$$\ S(t)= \sum_{n=-\infty}^{\infty} \delta(t-nT)$$

Now if you take the fourier series of S(t),you can write S(t) as

$$S(t) =\sum_{n=-\infty}^{\infty} C_ke^{jnw_ot} $$

Where $C_n$ are exponential fourier series coefficients and $w_o$ is the fundamental frequency .

So from exponential fourier series we know that

$$C_n= (1/T)\int_{-T/2}^{T/2} S(t)e^{-jnw_ot} dt$$

Now in the above expression substitute the value of S(t) from the first expression.

So $$C_n = (1/T)\sum_{n=-\infty}^{\infty}\int_{-T/2}^{T/2} \delta(t-nT)e^{-jnw_ot} dt$$

Now, you have to make an observation, if you observe the integral, it's from -T/2 to +T/2. During this integral period, observe that only a single impulse $\delta(t)$ exists.All the other impulse functions in the summation occur after T/2 or before -T/2. So in total the above equation for $C_n$ can be written as

$$C_n = (1/T)\sum_{n=-\infty}^{\infty} \delta(t)e^{-jnw_ot}$$

From sifting property we can write the above as

$$C_n = (1/T)e^{-jw_on(0)} = (1/T)$$

Now put this value of $C_n$ in the first S(t) equation

$$S(t) = (1/T)\sum_{n=-\infty}^{\infty} e^{jnw_ot}$$

Now find the fourier transform of above equation

$$1\Longleftrightarrow 2\pi\delta(w)$$

$$e^{jw_ot}\Longleftrightarrow 2\pi\delta(w-w_o)$$

So the fourier transform is
$$S(jw) = (2\pi/T)\sum_{n=-\infty}^{\infty} \delta(w-nw_o)$$

This should help.


How would I go about designing a difference equation for a given type of filter? I have experience using design tools such as Xilinx's CoreGen and MatLab's DSP tool kit, but I have not really had much experience designing filters from the ground up. Reinventing the filter as it were.  I prefer to use FIR filters because they are easier to implement and take less fabric.

I know that the filter command from MatLab implements a difference equation and give the filtered results based from the inputs given, but how does it calculate this difference equation?

How would I go about designing say, a high pass filter when I know the frequency I want to cutoff at, and the frequency I am sampling at?

I get a good majority of it, but I am still learning, so my question is :

How do I design a high pass filter with a difference equation so that I could implement it in something like C, and could I do this in the time domain, or would it have to be frequency domain only?

Matlab's filter() function implements a difference equation by calculating

$$\begin{align}a_0y[n]=&b_0x[n]+b_1x[n-1]+\ldots +b_Mx[n-M]\\&-a_1y[n-1]-a_2y[n-2]-\ldots-a_Ny[n-N]\end{align}\tag{1}$$

for every sample index $n$, where $y[n]$ is the output sequence, and $x[n]$ is the input sequence.
The coefficients $a_i$ and $b_i$ are stored in the vectors a and b. Note that normally the first element of the vector a is chosen as $a_0=1$. From (1) you see that in order to implement the difference equation you need to store past input and output values. This is usually done with a circular buffer.

Normally you don't talk about "designing a difference equation", but you design a filter, i.e. you compute the coefficients $a_i$ and $b_i$ in such a way that the filter performs some chosen task. There is a large variety of filter design algorithms available. Take a look at this answer for some basic background information on filter design (this is a huge topic in itself!).

In most cases you would use some filter design package (e.g. available in Matlab or Octave), design the filter (i.e. compute the filter coefficients), and then you can implement the filter according to Eq. (1). Writing your own filter design routine can be a nice exercise, but since filter design is a very mature topic, you'll basically be re-inventing the wheel. Implementing the filter is a different thing, and it can be useful to write your own implementation if you want to perform some real-time filtering.

Note that a direct implementation of Eq. (1) is only one of many ways to implement a filter. There are many possible implementations, all with their advantages and disadvantages. They all implement Eq. (1), but they perform computations in different orders, and they use different intermediate variables, which results in very different numerical behaviors. Of course, they would all be equivalent if it were possible to perform exact computations without any numerical errors. Especially for fixed-point implementations it is very important to choose an appropriate structure, because a direct implementation of Eq. (1) will be disastrous for larger filter orders (i.e. larger values of $M$ and $N$).

Note that Eq. (1) is a time-domain implementation. If you're interested in frequency domain implementations, search for "overlap-add" and "overlap-save" methods. However, I would start with implementing a time-domain filtering routine if you're working on your first implementation. Also note that for frequency domain filtering you need an FFT routine.

If you can use floating point arithmetic, I would as a first step try to implement Eq. (1) directly, and compare your results with Matlab's filter() function. For computing the coefficients you should use Matlab's filter design functions. If you then have any more specific questions you could come back and formulate a new question.