Classifying a time series with a Lorentzian power spectral density I was wondering if there is a general classification for time-varying signal with a Lorentzian shaped power spectral density, specifically with some nonzero center frequency. In terms of equations, what I mean is a power spectral density of the form $S(f) = \frac{a}{1+((f-f_0)/b)^2} + \frac{a}{1+((f+f_0)/b)^2}$ where $a$ is its maximal amplitude, $f_0$ is the center frequency and $b$ is the full width at full maximum. Note that I am talking about a real signal for all that follows.

The reason I am interested in this is because I know that if one has $f_0 = 0$, the timeseries could be seen as an instance of the Ornstein-Uhlenbeck (OU) process. However, I do not know of any process that generates a displaced Lorentzian power spectral density.

I tried thinking about this in terms of the autocorrelation function, but this has not led me very far. For OU this is an exponential decay, and by the Fourier transform this leads to the Lorentzian PSD. But conversely, a displaced Lorentzian again leads to exponentially decaying autocorrelations, multiplied by some phase factor that does not show up in the actual autocorrelations.

However, I am pretty sure there is a physical difference between the two power spectral densities (think about processes where the energy of the Lorentzian is relevant), and I would be interested in looking into the $f_0 \neq 0 $ case a bit further. So I am therefore wondering if there is some OU-type process that generates this, or some classification such as white noise, pink noise, etcetera, that fits this type of spectrum.

While it doesn't exactly give the power spectral density formulated by you, a bimodal (when considering also negative frequencies) power spectral density is given by the time evolution of the position or displacement of a harmonic oscillator driven by Brownian motion. The only modification compared to the Ornstein–Uhlenbeck process is inclusion of inertia, which enables oscillations when the harmonic oscillator is underdamped. An Ornstein–Uhlenbeck process is governed by the differential equation:

$$\dot x(t) = -\theta x(t) + F(t) \tag{1}$$

The dot above the $x(t)$ indicates the derivative of $x(t)$. In an infinitesimal time step $dt$, displacement $x(t)$ relaxes a step $\theta\,dt\,x(t)$ towards zero and is also affected by white noise external forcing $F(t)$. There is no inertia: neither the relaxation term or the external forcing term describe "memory of past velocity". I have poor differential equation handling skills, but my hunch is that if a bit of the past velocity $\dot x(t-dt)$ was mixed into the current velocity, it would give this modified form with inertia:

$$\ddot x(t) = -\gamma_m \dot x(t) - f_m^2 x(t)+F(t) \tag{2}$$

Here $\ddot x(t)$ denotes the second derivative of $x(t)$. Now there are terms both for damping of velocity and for relaxation of the displacement. For the modified process, the power spectral density is of form (Gröblacher 2012, Coffey & Kalmykov 2012):

$$S(f) = \frac{a}{\left(f_m^2-f^2\right)^2 + f^2\gamma_m^2}, \tag{3}$$

where $\gamma_m$ is the width of the resonant peak at half maximum, and $a$ is a normalization constant which has also absorbed the power spectral density of the driving force which is white (constant in frequency) for a Brownian heat bath. The true peak frequency $f_\text{res}$ is (Gröblacher 2012):

$$f_\text{res}=f_m\sqrt{1-\frac{\gamma_m^2}{2f_m^2}}.$$

Figure 1. Power spectral density of the displacement of a massive (as opposed to massless) harmonic oscillator in a heat bath, with $a = 1$, $f_m = 1$ and $\gamma_m=\frac{1}{2}$.

$S(f)$ of Eq. 3 is identical in form to the square of the magnitude frequency response of a two-pole resistor-inductor-capacitor (RLC) filter:

$$|H(if)|^2 = \frac{a}{(f^2 - f_0^2)^2 + f^2\,(f_0/Q)^2}$$

If the filter is fed white noise, the power spectrum of the output will be identical to $S(f)$. It is characterized by a -12 dB/octave (exactly -40 dB/decade) asymptotic slope:

Figure 2. Square of the magnitude frequency response of the two-pole filter, with logarithmic axes. The constants are equivalent to those in Fig. 1.

I'll return to the "just add inertia" argument once more. Consider a discretization of the Ornstein–Uhlenbeck process defined in Eq. 1:

$$\dot x[k] = c_0 x[k] + F[k]$$

We can define velocity and acceleration as:

$$\dot x[k] = x[k] - x[k-1]$$
$$\ddot x[k] = \dot x[k] - \dot x[k-1] = x[k] - 2 x[k-1] + x[k-2]$$

If we add damped inertia to the process as a term $c_1 \dot x[k-1]$:

$$\begin{array}{rrcl}
&\dot x[k] &=& c_0 x[k] + c_1 \dot x[k-1] + F[k]\\
\Rightarrow& x[k] - x[k-1] &=& c_0 x[k] + c_1 x[k-1] - c_1 x[k-2] + F[k]\\
\Rightarrow& 0 &=& (c_0 - 1) x[k] + (1 + c_1) x[k-1] - c_1 x[k-2] + F[k]\\
\Rightarrow& c_1 \ddot x[k] &=& (c_1 + c_0 - 1) x[k] + (1 - c_1) x[k-1] + F[k]\\
\Rightarrow& c_1 \ddot x[k] &=& (c_1 + c_0 - 1) x[k] - (c_1 - 1) x[k-1] + F[k]\\
\Rightarrow& c_1 \ddot x[k] &=& (c_1 - 1) \dot x[k] + c_0 x[k] + F[k]\\
\Rightarrow& \ddot x[k] &=& \frac{c_1 - 1}{c_1} \dot x[k] + \frac{c_0}{c_1} x[k] + \frac{1}{c_1}F[k]\\
\end{array}$$

We have arrived at the discrete equivalent of Eq. 2.

References

Gröblacer, Simon (2012) "Quantum Opto-Mechanics with Micromirrors: Combining Nano-Mechanics with Quantum Optics" DOI: 10.1007/978-3-642-34955-3, chapter 2, page 7.

Coffey, William T.; Kalmykov, Yuri P. (2012) "The Langevin Equation: With Applications to Stochastic Problems in Physics, Chemistry and Electrical Engineering" 3rd ed. ISBN: 978-981-4355-66-8, page 241.


CZT: Do Bluestein algorithm perform backward transformation? I am in discussion with one of my colleague that he says changing sign of the exponent in Bluestein's algorithm for Chirp-Z transformation, transforms the spectrum from frequency to time domain?
I am sort of confused, since they have invented ICZT very recently and it works in a different way.

It is easy to calculate IDFT using DFT, by applying the property:
$$\operatorname{IDFT}(x) = \frac{1}{N}\big(\operatorname{DFT}(x^*)\big)^*,\tag{1}$$
where $^*$ denotes elementwise complex conjugation. You can find this equation in a more verbose form in Eq. 13-55 of Richard G. Lyons "Understanding Digital Signal Processing" 3rd edition.
Probably this was used along with conjugation properties of the exponential.


MATLAB: Fourier transform vs $\tt lsim$, results don't match My (simplified) problem:
Suppose I have a signal $u$, which passes through system/filter $G$, which gives output signal $y$. Then I have a theorem that states:

If a signal $u(t)$, with Fourier transform $U(e^{j \omega})$ is applied to system $G$, then the Fourier transform of the output $y(t)$ is given by:

$$Y(e^{j \omega})=G(e^{j \omega})U(e^{j \omega})$$

I try to verify this in MATLAB, but something goes wrong.

Implementation:

$U(e^{j \omega})$ is obtained easily enough using fft(u).
For $G(e^{j \omega})$ I use freqresp(sys,$\omega$) (and build
$\omega$ in such a way that it is in the same format as the fft
results: [dc positive-freqs negative-freqs]). I verify this by comparing the result to the bode plot (see commented lines in script)

To verify I perform the inverse Fourier transform y = ifft(Y) and compare the obtained $y(t)$ to the lsim result, which is obtained with $y_{prime}(t)$ = lsim(G,u). Unfortunately these dont match. What am I missing?

clear
close all

G = tf([4700 4393 3.245e08],[1 7.574 1.202e5 0 0]);

f_sh = 1000;                                    %sampling frequency

Gd = c2d(G,1/f_sh);                             %plant discretization

t_end = 10;                                     %[s] simulation time end
t = linspace(0,t_end,t_end*f_sh+1);               %time vector
f_ref = 0.5;                                    %reference freq
u = sin(2*pi*f_ref*t);                        %reference signal

L = length(u);
dF = f_sh/L;                                    %frequency bin step

w = dF:dF:f_sh/2;                               %one-sided frequency vector in Hz

w_eval = 2*pi*[0 w -fliplr(w)];                 %frequency vector for freqresp, formatted to match fft results

H = freqresp(Gd,w_eval);                        %frequency response of system wrt freq vec

H = transpose(squeeze(H));                      %remove singelton dimensions

% bode(Gd) %should be equal to the semilog plot
% figure
% semilogx(w,20*log10(abs(H(2:(L+1)/2))))             %plot in Hz

U = fft(u);
Y = H.*U;
y = ifft(Y);                                    %frequency domain response

yprime = lsim(Gd,u);                              %simulation/time domain resonse

figure
plot(t,yprime,t,y)                           %y should equal yprime
legend('lsim response','fft(sys)*fft(ref) response')

The lsim results appear to have some integrator effect in them???

To perform the FFT approach correctly, you will need the FFT length to be at least $N_U + N_H - 1$ where $N_U$ is the length of $u$ and $N_H$ is the length of $h$.

It looks like $h$ is an IIR filter, so the FFT approach probably won't work well.

The Mathworks has a nice example comparing circular with linear convolution.

IIR filters are best implemented in the time domain. For their specs, they are usually much lower order than FIR filters and so we can afford to implement them directly.


Record Signal after Threshold Reach in GNU Radio I would like to use GNU radio to begin recording a predetermined number of complex samples to a file sink automatically as soon as a signal magnitude threshold is reached. For example, as the gnu radio flow graph is running, it should be monitoring a signal and waiting for the amplitude to reach a threshold. As soon as the threshold is reached, the next 10,000 complex samples should be written to a file sink. No more samples should ever be written to the file.

What is the best way of doing this?

Is this possible to accomplish using a custom embedded python block that has one input and one output that only allows a predetermined number of samples through to the output after the threshold has been reached and never allows any more samples through?

Sounds like you have a very workable approach:

Write a GNU Radio, Embedded Python, out-of-tree Python or C++, doesn't matter, which:

is a general block (not a sync_block)
has a member property triggered or similar, which is initialized to False in the constructor
has a member property threshold or similar, which is initialized to the value passed as constructor argument in the constructor
If triggered, then

copy input to output (specifically, min(len(input_items[0], noutput_items)), because in a general block, there's no guarantee that in- and output buffers have the same free space)
consume_each the number of items copied (tell GNU Radio you've processed this much input), and
return the number of items copied (you produced this much output).

else,

look for a super-threshold value in the input (numpy.any(numpy.abs(input_items[0]) > threshold)
if no values are found:

consume_each(len(input_items[0]))
return 0 (you didn't produce any output)

if a value is found:

triggered = True
find the position of the first value larger than the threshold
consume_each as much items as are before that position
return 0 (you didn't produce any output, but your work method will be called again, with triggered set)

... and connect its output to a "Head" block, set to pass 10000 items; connect that to your file sink.


Implementing a finite impulse respone (FIR) filter for computing the derivative of a discrete signal I am trying to reimplement an algorithm on my own. In the description of the implementation, it's written that they compute the derivate of a series of value using a [-1/2, +1/2] finite impulse response filter for obtaining an array of the original length.
Basically, a FIR filter should be a convolution but I don't understand this description applied to this case.
I've tried to develop a simple algorithm using an interpolation:
import numpy as np
from scipy.interpolate import interp1d
x = np.arange(len(signal))
signal  = [2,2,2,3,4,5,5,5,3,4,2,2,2]
# classical differentiation
diff = np.diff(signal)

f = interp1d(x, signal,fill_value='extrapolate')
newx = np.arange(-0.5, len(signal)+0.5)
ynew = f(newx)
intdiff = np.diff(ynew)

And to compare it with the classical differentiation:
import matplotlib.pyplot as plt
plt.plot(x,signal,marker='x',color='g',label='original signal')
plt.plot(x[:-1]+0.5,diff,marker='x',color='r',label='differenciated')
plt.plot(x,intdiff,marker='x',label='interpolation differenciated')
plt.legend()

I am wondering how to implement a derivate of a series of value using a [-1/2, +1/2] FIR filter? And what are the advantages compared to a derivative compute using I[n+1] - I[n].

FIR filter with impulse response $h[n] = {-\frac{1}{2}, \frac{1}{2}}$ means that :
$$y[n] = \frac{1}{2}x[n-1] - \frac{1}{2}x[n]$$
This is in some sense a mirror operation of Moving Average of two samples. This is Moving Difference (samples reversed). A High Pass Filter.
So, each output sample is the difference of current input sample with previous input sample. And, this is the result of convolution operation between $x[n]$ and $h[n]$.
But I am not sure that this is will correctly give the Derivative, since for a positively increasing discrete ramp, you would expect the discrete derivative to give positive result. Because for a increasing discrete ramp $x[n]>x[n-1]$. And that FIR filter will give negative result and in addition to negative slope it will also scale the result by $\frac{1}{2}$, which is unwanted.
Correct definition of a Discrete Derivative would be in my understanding, the following:
$$\dot{x} = x[n] - x[n-1]$$
And this can be implemented using FIR filter with impulse response $h[n] = 1,-1$


Filter common signal from image In a set of images taken with the same infrared camera, I have a noise signal presumably originating from the camera chip. How would I go about extracting said signal to be able to substract it from all images I have taken?

Basically you have to do a shading correction. This process removes known artifacts from an image (or a series of images) based on a reference image (or "compensation image" as it is sometimes called) containing only those artifacts. Such a reference images can be obtained by photographing a uniform object (i.e. a perfectly white wall or perfectly black darkness) under good conditions so that only the signal from the defect camera chip is recorded.

I have done this example with Mathematica, the basic principle can be implemented anywhere.

Given a set of images (I used grayscale images since I have no infrared images) without the defective signal

img1=

img2=

img3=

and an artificial defect camera chip signal aka the reference image (I pulled this out of my hat).

noise = GaussianFilter[Binarize[Image[RandomReal[1, {299, 499}]], 0.95], 5]

Now we generate the noisy images.

img1n = ImageAdd[img1, noise]
img2n = ImageAdd[img2, noise]
img3n = ImageAdd[img3, noise]

Then we simply subtract the reference image from the noisy images and get the actual images without the noise.

ImageSubtract[img1n, noise]
ImageSubtract[img2n, noise]
ImageSubtract[img3n, noise]

If you happen to have a noisy image and a perfectly fine image with exactly the same scene (and of course image sizes)

img1n =

img1=

then you can simply generate the reference image yourself by subtracting the latter from the former.

ImageSubtract[img1n, img1]


Calibration - acoustics I have one recorded signal (.wav) and the calibration signal (.wav - 94dB @ 1kHz).

I want to know how to calibrate the recorded signal using the calibration signal, for example, in Matlab.
My tasks are

estimating background noise level,
equivalent levels

Aeq,
Ceq,
Zeq,

max-values of

S and
F time

integrations

ASmax,
AFmax and

peak values

Apeak,
Cpeak,
Zpeak.

First you need to know something about your calibration signal. Lets say you know that it is a 1kHz sine with an RMS amplitude of 94dB SPL. If you are confident that this is all your calibration recording includes (which is likely), then your calibration factor, is the difference between the actual dB SPL level and the 'digital' dBFS level. In matlab code, it looks like this:

caldb = 94 - 20*log10(rms(calibration_signal);

It could be that the calibration level is specified as peak SPL, then you would just substitute the peak value for the rms value above.

Now you can calculate all your values (LAeq etc) in dBFS (i.e. -inf to 0dB) and simply add the calibration factor, to get a calibrated SPL.


contrast curves problems I am doing some tests about local/global contrast algorithms and need a simple equation to define an arbitrary S-shaped contrast curve. It's simple if the turning point is assumed to be (0.5, 0.5), but couldn't find a straight equation if the turning point has to be arbitrarily set.

These are the typical contrast curves I achieved for turning point in the middle (0.5, 0.5):

And the simple equations follow this rule:

y = x^2*2 for x<=0.5

y = 1-(1-x)^2*2 for x>0.5

I would like to be able to allocate the turning point to any (x,y) position and still have a smooth S-shaped curve where the slope is continuous on both sides of the turning point. I would like to keep the equation as simple and straightforward as possible, avoiding the use of splines or any other more complex math.

Any ideas? I tried to google I found some equation like this.

y = [ (1-s)x + sa*(x/a)^E ] ^ ( log(b)/log(a) )

Parameters:
(a,b): turning point

s: slope at ends

E: contrast strength

But it does not work like this

Would you please let me know how do I make the equation like above.

One good way is to utilize Bezier curves (or Catmull Rom splines). This way, you could define the control points in the way you like. Positioning of the control points would then allow you to adjust the histogram to the desired output.

I would suggest Catmull Rom splines if you want them to pass through the control points. It gives you better control over the curve behavior. Use 5 such points: 2 for the end points, 1 for the middle, and 1 on each side of the center to give you certain curvature.

If you don't want to use splines, the black curve looks like a Sigmoid function, while the yellow is a third order polynomial.

You can checkout this code to see different implementations.

If you are keen on using Logistic family, then I would suggest the generalized logistic growth curve with horizontal asymptotes at $−1$ and $1$:

$y(t)=-1+\frac{2}{[1+ke^{-\beta t}]^{1/v}}$

where $k > 0$, $β > 0$ and $ν > 0$. The parameter $k$ can be adjusted by a time shift, $β$
corresponds to a scaling of time; $ν$ is the key parameter that determines the shape of the
growth, the asymmetry depending on whether it is less than or greater than $1$. This formula reduces to the standard logistic growth for $ν = 1$, $k = 1$, $β = 2$. The turning point you refer to is commonly known as the inflection point, where the derivative changes sign.

This is explained in detail Ayse Humeyra's article and here. Here is a MATLAB code:

function [y]= general_logistic(t, k, beta, v)
y = -1 + 2./((1+k*exp(-beta.*t)).^(1./v));
end

Finally, in Dino Dini's blog, he explains the normalized tunable sigmoid functions. There, the solution to your problem is addressed. Moreover, a similar question is already posted and answered here. Last but not least, polynomial shaping functions, which would serve the same necessities are discussed in Golan Levin's page.

I hope all these would help.


Amplitude Spectrum, Nyquist Frequency, mixed/min/max wavelets The problem is here. Now I know the definition of mixed/max/min phase wavelets, whether the roots lie within the unit circle or not. Starting from n = 1, let
$$ x_t = ( 5, 6) $$
$$ X(z) = 5 + 6z $$
This would have root -5/6, making it a max phase wavelet (since it's within unit circle). If we switch 5 and 6 around in xt, we would get root -6/5, which would make it a min phase wavelet.

I don't know how I would show they have same amplitude spectrum. I know |X(f)| amplitude spectrum would lie between -fc and fc (Nyquist frequency) and fc = 1/2*delta t. Not sure how I would go about proving it.

Note that moving between minimum phase, maximum phase, and mixed phase involves basically two operation in the time domain:

shifting
time-reversal

If $X(e^{j\omega})$ is the Fourier transform of the sequence $x(n)$ then the transforms of the shifted and the time-reversed sequences are given by

$$x(n-k)\Longleftrightarrow e^{-jk\omega}X(e^{j\omega})\\
x(-n)\Longleftrightarrow X(e^{-j\omega})$$

Since $|e^{-jk\omega}|=1$, the magnitudes of the spectra of $x(n)$ and $x(n-k)$ are identical. Furthermore, if $x(n)$ is real-valued then $X(e^{-j\omega})=X^*(e^{j\omega})$ and consequently the spectra of $x(n)$ and $x(-n)$ also have the same magnitudes.


Direct form I and Direct form II I was studying about realization structures of digital filters. Is it mandatory to have the order of numerator must be less than that of  denominator of transfer function for realization of filter using direct form I and II?

no.  but, for a causal LTI digital filter with a rational transfer function, it is necessary that the number of zeros, $M$, not exceed the number of poles, $N$.

$$\begin{align}
\\
 H(z) \triangleq \frac{Y(z)}{X(z)} & = H_0 \frac{(z-q_1)(z-q_2)(z-q_3)...(z-q_M)}{(z-p_1)(z-p_2)(z-p_3)...(z-p_N)} \\
\\
 & = H_0 \frac{z^M (1-q_1 z^{-1})(1-q_2 z^{-1})(1-q_3 z^{-1})...(1-q_M z^{-1})}{z^N (1-p_1 z^{-1})(1-p_2 z^{-1})(1-p_3 z^{-1})...(1-p_N z^{-1})}  \\
\\
 & = z^{M-N} \ \frac{H_0 \ (1-q_1 z^{-1})(1-q_2 z^{-1})(1-q_3 z^{-1})...(1-q_M z^{-1})}{(1-p_1 z^{-1})(1-p_2 z^{-1})(1-p_3 z^{-1})...(1-p_N z^{-1})}  \\
\\
 & = z^{M-N} \ \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + ... + b_M z^{-M} }{ 1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + ... +  + a_N z^{-N} }
\end{align}$$

this gets manipulated:

$$\begin{align} Y(z) & \left(  1 + a_1 z^{-1} + a_2 z^{-2} + ... + a_N z^{-N}  \right) \\
& \quad = X(z) z^{M-N} \left( b_0 + b_1 z^{-1} + b_2 z^{-2} + ... + b_M z^{-M} \right) \end{align}$$

or

$$\begin{align} Y(z) = \ & b_0 z^{M-N} X(z) + b_1 z^{M-N-1} X(z) + b_2 z^{M-N-2} X(z) + ...  \\
 & \quad ... + b_M z^{-N} X(z) - a_1 z^{-1}Y(z) - a_2 z^{-2}Y(z) - ... - a_N z^{-N}Y(z) \end{align}$$

inverting the Z Transform, this becomes

$$\begin{align} y[n] = \ & b_0 \ x[n+M-N] + b_1 \ x[n+M-N-1] + b_2 \ x[n+M-N-2] + ...  \\
 & \quad ... + b_M \ x[n-N] - a_1 \ y[n-1] - a_2 \ y[n-2] - ... - a_N \ y[n-N] \end{align}$$

for a causal system or filter, we cannot look into the future for any input sample $x[n]$.  that means $M \le N$ so that there is no negative delay on the $x[n+M-N]$ term.  a negative delay means the LTI filter is peeking into the future, which is not causal.

but if you remove that leading delay factor and have this:

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + ... + b_M z^{-M} }{ 1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + ... +  + a_N z^{-N} } $$

there is no reason that $M$ cannot be larger than $N$. but the number of poles will still be as large as $M$. if $M>N$, then some of the poles will have to be at the origin $z=0$. in fact, if $N=0$, then you have an FIR filter and all of the poles are at the origin.


How to evaluate fixed-point implementation of LMS filter is correct? I am having an LMS block with 6 filter coefficients.
The value of filter coefficients are
0.0001
0.00045
0.2535
0.546536
0.0000243
0.3423

I have tried to implement the LMS algorithm in floating-point with $\mu=0.01$.
LMS algorithm implemented in floating point is based on the below equations:(from figure) :
v(n)=wT(n)u(n),
e(n)=d(n)-v(n),
w(n+1)=w(n)+mu e(n)*u(n) ,which is the weight updation equation.
Here d(n) is the desired signal,u(n) is the input to lms block,v(n) is the lms filter output ,w(n) is the adaptive filter weights and e(n) is the error.
I have given a set of desired samples and input samples to the LMS filter and the LMS filter output data samples were analyzed. Now, to implement the same LMS filter in fixed-point representation, I will have to do scaling and rounding. Is scaling of filter coefficients by some factor $2^{s_1}$ enough and $\mu$ value be represented like $2^{-6}$ or should we consider any other factors while doing scaling operations. Finally, how to check the fixed-point implementation of LMS filter worked correctly?

I'm not an expert on the LMS algorithm. Perhaps you should add a link to the algorithm description so we can help you. However, I have adapted a lot of algorithms to fixed-point implementation, so I think I can help you.
I'm gonna give you a few pointers :
1 - If your algorithm is not recursive, quantization errors can affect the accuracy of your results but the fixed-point implementation will always be stable.
2 - If your algorithm is recursive, it is possible that quantization will radically change the behaviour of your algorithm, it could become unstable, you could be stuck with a limit cycle, etc.
3 - Create a "perfect" model. In your case, the floating-point implementation could be your perfect model. Then, convert your model to fixed-point and compare the results of your fixed-point model. You can toy around with the number of bits for your coefficients, rounding versus truncation, etc.
Regarding your algorithm, you have 6 coefficients, 3 of them are significant (coefficients 3, 4 and 6) while 3 of them seem insignificant (1,2,5). Is this ok?


Autocorrelation sequence in terms of Fourier transform of the underlying signal Let $x(n)$ be a sequence of length $N$, which is zero outside the interval $(0,N-1)$. Let $X(k), k=0,1,\cdots,N-1$ be the FFT coefficients of $x(n)$, that is, $X(k)=\sum_{n=0}^{N-1}x(n) \exp\left( -\frac{j2\pi k n}{N}\right)$. How to relate the autocorrelation sequence  $y(l)=\sum_{n=0}^{N-1-l}x(n)x^{*}(n+l)$ in terms of $X(k)$?

My try:

\begin{equation}
\begin{split}
y(l)& =\sum_{n=0}^{N-1-l}x(n)x^{*}(n+l)\\
& =\sum_{n=0}^{N-1-l}x(n) \left[  \sum_{k=0}^{N-1}X^*(k) \exp\left( -\frac{j2\pi (k+l) n}{N}\right) \right]\\
& = \sum_{k=0}^{N-1}X^*(k) \left[\sum_{n=0}^{N-1-l}x(n)    \exp\left( -\frac{j2\pi (k+l) n}{N}\right) \right]\\
\end{split}
\end{equation}

How to proceed from here?

There is, unfortunately, no simple way of getting to where you want to from what you have.  As Marcus Muller's comment says, the convolutions and correlations that are associated with the DFT (FFT if you like) are cyclic or periodic whereas the autocorrelation that you seek is the linear or aperiodic autocorrelation, and it is difficult to derive the latter from the former.  Note for example that cyclic shifts e.g. $(x_k, x_{k+1},\ldots, x_{N-1},x_0, x_1, \ldots x_{k-1})$ instead of $(x_0, x_1, \ldots, x_{N-1})$ have the same periodic autocorrelation but quite different aperiodic autocorrelations. The classic example of this is a PN sequence (maximal-length sequence) of length 7 which has the classic "inverted thumbtack" periodic autocorrelation but quite different aperiodic autocorrelation function depending on which cyclic shift we are using (a.k.a. what is the initial loading of the linear feedback shift register (LFSR) that is generating the PN sequence. In particular, one of the cyclic shifts of the PN sequence is a Barker sequence whose aperiodic autocorrelation function has maximum magnitude $1$ only $\big($that is, all the autocorrelation values are in $\{-1, 0, 1\}\big)$ while other cyclic shifts have larger autocorrelation magnitudes.

So, what to do?  The standard trick is to zero-pad the sequence.  If you want to find the aperiodic autocorrelation function of $(x_0, x_1, \ldots, x_{N-1})$, begin with the sequence
$$\hat{\mathbf x} = (x_0, x_1, \ldots, x_{N-1}, \underbrace{0, 0, \ldots, 0}_{\text{at least} ~N-1~0\text{s}})$$ of length $M \geq 2N-1$, compute its length-$M$ DFT etc. There are plenty of answers here explaining the zero padding technique and so I won't bother explaining further.


What happens when you wrap a 3x3 kernel at the corners of an image? I am puzzled by what actually happens when you process border cases when applying an image filter. For instance if I apply a 3x3 Gaussian kernel to the top left pixel[0][0], the remaining pixels that are unable to convolute will go to the right and bottom of the image. But do you multiply them again or u replace the original pixels?

[1, 2, 1,
 2, 4, 2,
 1, 2, 1] // 4,2 & 2, 1 will convolute with the pixels in image but how about row 1,2,1 and col 1,2,1 ?

There are three main attitudes when dealing with borders of images and filter kernels.

The simplest is to only process the pixels that have enough neighboring pixels. As a result, some pixels on the sides are left unprocessed, and the filtered image can be cropped.

The shape of the  filter can be adapted (notably in size) to fit the neighborhood of of the running pixel.
The image is extended to the appropriate size. There are many extension methods, from filling with zeroes or constant values to clever extrapolation.

Traditionally, symmetric or antisymmetric extensions are used in many cases, since this is quite consistent with the symmetries or antisymmetries of filter kernels themselves, and this can be performed efficiently. Wraparound is also performed, and sometimes implicit when using FFT filtering or discrete wavelet restoration. Due to discontinuities (the far ends of the image can be quite different), a preprocessing is advised (windowing, smoothing, background removal).


DSP based FMCW radar I have an SDR receiver I'm interested in making into a FMCW radar. The bandwidth is somewhat narrow as it only samples at 64 MSPS.

There is a problem in that without a bunch of embedded programming the RX and TX chains won't have their timing synchronized. Basically the up-conversion oscillator is synchronized but the DSP portion on the TX and RX lines won't be.

What I'm wondering is, is there a way to use the initial signal sent from the transmitter and received by the receiver through direct LOS to establish timing offset for the reflected signal? Given that the RX and TX antenna will be close together you could use this to synchronize the system.

Also I'm not sure what the processing would consist of though(or if it's even viable), most FMCW radar mix the transmitting signal with the reflected signal and then downconvert to produce to a local beat. I would be attempting to do this through DSP and have a much narrower bandwidth useable.

We did the same a couple of years back with what we called "echotimer", but not form FMCW, but for OFDM radar: basically, if you have a very large target that's very close and at 0 doppler, that's your crosstalk, and you can just subtract the distance of that.

Same applies here: with TX and RX hopefully sharing the LO, you're looking for a zero-velocity target that's close. You determine your range based on that, it represents your RX-TX time offset, and you just subtract.


Alternatives to Hough transform for detecting a grid-like structure I have an image which is composed of multiple 'angles' which form a grid shape:

After some searching, Hough seemed like a good fit, because it isn't troubled by breaks in lines. However, the problem I have is that my lines are 'fat', and any edge detection I run (Canny in this case), picks out the edges of the line, and not the middle:

This means that the Hough transform ends up picking out one (or both) of the 'sides' of my grid lines, and not the middle.

Given that I know what I'm looking for (this grid-like shape, always in roughly the same orientation), is there a better way to perform the 'edge detection' part to give me the centre lines, or should I be looking at a totally different method?

An alternative to the Hough Transform would be the Radon Transform (1, 2). A rough description of an algorithm to detect a grid-like structure could look like this:

1. Perform Radon Transform from 0 to 180 degrees.
2. Find the two highest maxima in the angle bins.
3. For the two angles with maximal amplitude find the local maxima within the bin.
4. You can use the constraint that the maxima should have even spacing to deal with outliers.

EDIT:

Here is a small matlab snippet to illustrate step 1-3:

im = imread('grid.png');
[R, xp] = radon(im, 0:180);
imagesc(0:180, xp, R)

plot(max(R)) % the two maxima are at 65 & 117 degrees

plot(R(:, 65))

plot(R(:,117))

To answer your question from the comments: It appears to me from the one example image that you provided, that this method is more robust against small defects in the detected grid. Skeletons will rarely yield straight lines which might be a disadvantage for the subsequent Hough Transform.


How to set a Carrier Wave for DBPSK Modulators in gnuradio? for a Project I need to build an Audio Modem in GNU Radio.

I decided, that it would be the best to use DBPSK Modulation, because that is fast and relatively noise resistant.

But as it seems, the DBSK Modulators are just modulating the Signal on the Baseband.

That's not what  I want. Would it be enough, to just multiply the signal with a signal source? And when I do that, do I have to meet any preconditions so that the demodulator doesnt expect a Baseband?

PS: If you have any recommendations what would be a good filter, to get the best out of it, I would like to hear that to.

Short answer: yes; digitally mixing this up to a carrier frequency and just using the real part will suffice. For more info, refer to the usual "complex baseband" introduction.

Your filter question is completely impossible to answer – it depends on the characteristics of your transmission and your audio channel, and you mention nothing of that.

Also, you say DBPSK is fast. I don't know what you really mean with "fast", here, but I assure you that a 1bit/symbol modulation isn't per se a very efficient one. Unless you're in a very special case, you typically stay far away from being a good modulation – it's easy to demodulate, indeed, because you don't need knowledge about the initial state and hence can vastly simplify synchronization, but that makes it an especially spectrally inefficient modulation. That means you waste bandwidth that you could use for more bit/second in favour of a simpler system.


Estimating the complexity by the calculation of the number of real multiplications in a general polynomial function For the following memory polynomial equation,

$y(n)=\displaystyle\sum_{m=0}^{M-1}\displaystyle\sum_{p =1\\p\,\textrm{odd}}^{P}a_{m_p} \, x(n-m) \, |x(n-m)|^{(p-1)}$

The total number of coefficients is $ M\frac{(P+1)}{2}$

I need to calculate:

1- the total number of complex multiplications

2- the total number of real multiplications

The direct computation should be straightforward. Since I suspect the exercice has other purposes (more advanced algorithmic complexity), there additionally are a couple of tricks to reduce the complexity of power evaluations.

First, if you have to evaluate all powers $y^0,y^1,\,\cdots,y^Q$, you can store the power-of-two powers, and use the binary expansion of a power $q$. It can be called exponentiation by squaring. For instance, to compute $y^7$, you need one multiply for $y^2 = y\times y$, one for  $y^4=y^2\times y^2$ and two more for $y \times y^2 \times y^4$, so four multiplies instead of six by direct computation. And the $y^2$ and $y^4$  can be stored for the next powers.

Second, a classical polynomial evaluation is "Horner's scheme":

$$a_0 x^n + a_1 x^{n-1} + \cdots + a_n =  \left(\left(\left( a_0x+a_1\right) x \right)\ldots a_{n-1}\right)x+ a_n$$
which can be adapted to your case.


wavwrite with 32 bits per sample seems to overflow? I am trying to generate a 32bit PCM wav file with octave, the file shall contain a single full-scale sample (0db impulse). I am getting the expected result when using 24 bits per sample:

>> wavwrite(1, 96000, 24, 'foo.wav');
>> wavread('foo.wav')
ans =  1.00000                          # as expected

However, when using 32 bits per sample, the value seems to overflow.

>> wavwrite(1, 96000, 32, 'foo.wav');
>> wavread('foo.wav')
ans = -1                                # WHAT ?!

I am expecting the second wavread('foo.wav') to return 1 not -1. I have checked the wav file with a hex editor, the sample raw value is 0x80000000, I am expecting it to be 0x7fffffff.
Is it an octave bug or am I missing something?

Using 0.99999 instead of 1 generates the expected result.

It seems like it is a bug. Octave's wavwrite() function writes 32-bit wav in 32 bits int (type 1 format) instead of normalized 32-bit floats (type 3 format). Unfortunately, 1 can't be represented in 32-bit type 1 format which is why you get this result.
The solution could be to use the audiowrite() function, or to write 0x7fffffff instead of 1.
More about wav formats.
EDIT: seems from the bug report that the audiowrite function has the same problem...


Baseband vs Passband modulation I am taking a graduate level course in Networking theory and applications and I am still a bit shaky on the difference between baseband and pass band scheming. The book I am using is Computer Networks 5th edition -- Andrew S Tanenbaum

To my knowledge, passband scheming uses the changing of frequency, amplitude, and phase to modulate a analog signal. And base band uses encoding techniques like NRZ, NRZi, and Manchester encoding.
How does this work along side, frequency multiplexing and CDMA for multiple device access?

And can someone give me a layman explanation of carrier signals? And how I can maybe think of them as offsets like in assembly languages.
Thanks for helping!

There is no unique definition, but in general baseband signals have a spectrum that starts close to DC. Most analog information sources that are considered in communications are baseband: audio, video, signals from sensors, images, etc. Also, any signal that is a sum of pulses of different shapes is a baseband signal; examples are NRZ and Manchester-encoded signals. Yet another definition is that a baseband signal has not been modulated, or its spectrum has not been shifted, for the purpose of matching it to a particular communication channel.

Passband signals were originally baseband, but were modulated, or frequency-shifted (or up-converted) for some purpose, usually to make them fit in a passband channel. Consider the case of AM radio. An audio signal is low-pass filtered to 10 kHz, resulting in a baseband signal with a spectrum from 0 Hz (or strictly, 20 Hz) to 10 kHz. Then, this signal is shifted up in frequency (upconverted) to make it fit into the 20 kHz channel assigned to the AM station. The receiver down-converts the signal back to baseband and plays it over the loudspeakers.

As you say, the upconversion process involves changing the amplitude, frequency and/or phase of a high-frequency carrier signal proportionally to the information-bearing baseband signal. The receiver uses the amplitude, frequency and/or phase of the carrier to reconstruct the baseband signal.


What's the number of data points in a Welch-based PSD? In Welch periodogram Power Spectral Density estimate, we devide the N-long signal into K segments, each of length L with overlapping D, such that N=L+(K-1)*D.

In this paper, just above eqn (7), the authors mention that the step of performing the DFT is to be done to each (windowed) segment (of length L) such that the length of the resulting DFT is L/2. Finally all DFT's will be averaged and normalized.

At the end, you get a PSD of a max. length L/2. Why is that? I would expect the length of the resulting PSD to be related to N, not to L...

And that's absolutely correct. As you mentioned, procedure can be described as follows:

Take long signal of length $N$ and slice it with some overlap and window into segments of length $L$.
For each of these segments calculate the DFT. Assuming no zero padding, you will get spectrum with $L$ points - same as segment length.
Because signals are real, then we just care about first half of spectra (squred magnitude to be correct), thus take only $L/2$ points.
Do this procedure for each segment and you will end up with lot's of vectors (dependent of overlapping and signal length $N$) of length $L/2$.
Now you need to average all the spectra. This is simply average for each frequency bin across number of segments and PSD is calculated.

Let's assume you will end up with 50 segments of length $L$, then you must average 50 spectra of length $L/2$. The longer is your signal, more segments will be averaged (better estimate), but obviously their length is not changing - you've chosen it in the beginning.


Predicting Execution Time of DSP programs I've been tasked with predicting the performance of our DSP code, for the obvious reasons (reduce time to market, add certainty in choosing chips, etc.)

Consider a system composed of many subsystems with timing results $P_{i}$ (in seconds) (filters, gains, dynamics processing etc.) with a clockrate $C$ (seconds/sample) and sampling rate $Fs$ (samples/second)

Our initial assumption is that:

$$\text{CPU utilization} = 100\% * \sum_{i}^{i=N}{P_{i}*Fs/C} \stackrel{!}{=}
\texttt{top -d 1}$$

i.e. CPU utilization estimate should equal the output of top -d 1`.

However, this doesn't hold up to real data! It's always underestimated.

We measure individual subsystems by wrapping each of them in a test harness and timing their execution. Our system runs in a single thread and on a single core.

Depending on the platform, we can get within 5 to 20%, but always underestimating top.

Why is this?

I can find nothing in the performance literature that lines up with our assumption, but the explanations are so numerous that it's almost like "we don't really know why we can't predict performance."

I am beginning to think that the solution is monitoring like a hawk during early development and scaling back the system complexity as we uncover new solutions.

However, it would be pretty awesome if we could design a system with the performance requirements ahead of time. Our requirements have largely been determined by "this is better than the old thing the vendor had, and what we had before worked with our old system."

Top is terrible tool for measuring CPU time. Set up a bunch of counters or timer around the execution of individual and collect statistics on that.
Each of your block should have a "should cost" that's based on number of arithmetic operations plus some overhead for looping and conditionals relative to the instruction set of the processor. Compare each blocks actual cost to the "should cost" and look for major discrepancies
If the total doesn't add up, it's often caused by cache misses, pipeline stalls, memory bottlenecks, interrupts, or some other service barging in. Again, proper instrumentation can shed some light on what's happening.


Besides computational costs, is filtering, decimating, then filtering different from just filtering? Background: I have a variety of BPSK which uses a raised cosine transmitted pulse shape:

$$ h(t) = {1 \over 2}\, (1+\cos(\pi t))\, \Pi(t/2) $$

After matched filtering in the receiver, the end result has ISI. Decimated to 1 sample per symbol, the pulse shape works out to [1/6, 1, 1/6].

So, I might filter and decimate this signal to 1 sample per symbol, then apply an equalizing filter to mitigate that ISI.

My question: modulo the potential consequences to computational complexity, is this any different from combining the matched filtering and equalization into a single filter?

I know filtering by two successive FIR filters is equivalent to one filter with the impulse response of the two convolved together. But with the decimation step between, I'm less sure.

Ignoring the higher computational costs, it's possible to do what you want. Let's consider the case without ISI, where the channel response $C(f)=1$. You transmit the sequence $$s(t)=\sum_{k=0}^N a_k h(t-kT),$$ where $h(t)$ is your raised cosine pulse shape and $T$ is the pulse rate. Ignoring noise, the output of the matched filter at times $t_k=kT$ is the data sequence $a_k$.

If the channel response $C(f)$ is time invariant but not flat, then the received signal (before the matched filter) is $$r(t) = s(t) \ast c(t),$$ where $c(t)$ is the channel impulse response. Now you need to design a matched filter + equalizer combination filter $g(t)$ such that $r(t) \ast g(t)$, sampled at times $t_k$, is the data sequence $a_k$.

Note that all I have done is simply to delay the decimation until after the combination filter. You say that

  I know filtering by two successive FIR filters is equivalent to one filter with the impulse response of the two convolved together.

This is true when the filters are time-invariant. A decimator is time-variant, and it can't be trivially combined with other LTI filters.


OFDM violates sampling theorem The "textbook" OFDM block diagram is as the below figure. But I wonder if it satisfies the sampling theorem. Here, we have N points $X[k]$ before IFFT, and

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

This derives from the analog signal

$$x(t)=\sum_{k=0}^{N-1}\mathrm{e}^{\displaystyle -j 2\pi f_k t}=\sum_{k=0}^{N-1}\mathrm{e}^{\displaystyle -j 2\pi k\Delta f t} $$, where $\Delta f$ is the subcarrier space.

We then sample this signal with a  sampling frequency $f_s>2f_{max}=2(N-1)\Delta f$, which can be taken as $f_s=2N\Delta f$. The time domain sample points are at

$$nt_s=n/f_s=\frac{n}{2N\Delta f}$$

Then we define

$$x(nt_s):=x[n]=\sum_{k=0}^{N-1}\mathrm{e}^{\displaystyle -j \frac{2\pi k n}{2N}} $$

and then we define a expanded

$$\tilde{x}[n]:=\begin{cases}
x[n], 0\le n\le N-1\\
0, 2N-1\ge n\ge N\\\end{cases}$$
$$

That means we should do a length $2N$ IFFT rather than a length $N$ IFFT

$$
\tilde{x}[n]=\sum_{k=0}^{2N-1}\mathrm{e}^{\displaystyle -j \frac{2\pi k n}{2N}}$$

Is the textbook wrong?

This may sound trivial but the sampling theorem can only be violated in the process of sampling an analog signal. In your question you would like to sample the signal
$$
x(t) = \sum_{k=0}^{N-1}\mathrm{e}^{\displaystyle -j 2\pi k\Delta f t}
$$
and you state correclty that, in order to not violate the sampling theorem, $f_\mathrm{s}>2(N-1)\Delta f$ must hold. In contrast, OFDM systems usually work with a sampling frequency $f_\mathrm{s}=N\Delta f$ and you may now ask where this discrepancy originates.

If we look at $x(t)$ we notice that it has a one-sided spectrum with discrete frequencies at $0, \Delta f, \ldots, (N-1)\Delta f$. That's an odd choice because the signal bandwidth is $2(N-1)\Delta f$ but only $N$ subcarriers are used. With other words: this system has potentially $2N$ subcarriers but the subcarriers $-N, -N+1, \ldots, -1$ are unused and the corresponding spectrum is wasted. We should rather define $x(t)$ like this:
$$
x(t) = \sum_{k=-N/2}^{N/2-1}\mathrm{e}^{\displaystyle -j 2\pi k\Delta f t}
$$
The absolute maximum frequency is now $\frac{N}{2}\Delta f$ and accordingly $f_\mathrm{s}>N\Delta f$.

For creating a discrete-time OFDM signal, however, $x(t)$ is quite irrelevant. You cannot violate the sampling theorem when generating a discrete-time signal, because it can be considered as already sampled. Furthermore, due to the cyclic property of the DFT, the limits of the sum in the definition of $x[n]$ are correct. Whether they are $0 \ldots N-1$, $-N/2\ldots N/2-1$ or $N+3\ldots 2N+2$ doesn't matter.

So I wouldn't say that the textbook is wrong. But if OFDM is motivated by sampling an analog signal, then $x(t)$ as defined in your answer is a poor choice.


Matlab - Creating a characteristic for Peak and Shelving filter I have designed few shelving and peak filter. And I have calculated the a, b coeficients for several gains for these filters. Now I would like to create some kind of characteristic in Matlab so I could present my design of these filters?

Something like this for shelving filter

I suggest you to use fvtool in MATLAB:

fvtool(b,a)


Generating FIR filter in Web Audio API I'm creating an HTML that applies a low-pass Windowed Sinc filter with a window size of 8192 to white, pink, or Brownian noise. There are nine radio buttons, eight of which set a cutoff frequency and one which disables the filter.

I've got the noise generator working, with all three "colours" generating noise with no problem.

I've been looking for a way to implement a Windowed Sinc filter using the Blackman window in JavaScript. I've found one site that seemed to be pretty promising,
but on that site, the event listeners are implemented and the processing itself is done (in a script called "dfilter-0.js") which is extremely large and messy. There are so many unnecessary functions and chains of functions that I can never figure out where anything is or what anything does.

My question is, could I get some help cleaning up the "dfilter-0.js" (none of the online tools I've found worked)? Or alternatively, how does one implement a low pass Windowed Sinc filter with the parameters (cutoff frequency, Blackman window size) and connect it to an AudioContext?

I've found an Audacity plugin that performs the filtering that I want, but it uses the Nyquist language. How would I implement that in JavaScript and HTML?

A FIR filter is just a "block" of coefficients, a vector, being convolved with your signal.

Hence, https://developer.mozilla.org/en-US/docs/Web/API/ConvolverNode is kind of audio filtering object you need to create and apply.


What Is Proper Paper and MATLAB Toolbox Describing and Discussing Wavelet, Ridgelet, Curvelet and Contourlet Transforms in Detail? Question 1:

Does anybody know any paper proper for getting started understanding the concept of Wavelet, Ridgelet, Curvelet and Contourlet transforms and their differences in detail?

Question 2:

does anybody know a MATLAB toolbox for implementation of these transforms?

On A panorama on multiscale geometric representations, you find a tutorial paper on those directional 2D wavelets, starting with an historical perspective (eg Gaussian & Laplacian pyramids):

and highlighting details between the main fixed and adaptive decompositions (below: a contourlet atom).

It was published in 2011 in Signal processing, special issue on multirate filter banks and multiscale representations.

It also links to the two-dimensional wavelet and multiscale toolboxes used to generate the figures.


Relative Computational Complexity for Channelizer Implementations Suppose I want to separate a broadband input signal into multiple narrow subbands using an FFT-based analysis filter bank. The filter bank uses a prototype lowpass filter and is implemented using a polyphase structure.

See Matlab's Channelizer class

Does it require fewer FLOPS (floating point operations) to do the filtering in the frequency domain as a multiplication and sum or as a convolution in the time domain? Matlab implements the filtering in the time domain.

For sufficiently long filters, it's usually more efficient to implement the individual filters as FFT-based fast convolution, indeed!

Now, whether something is faster than something else really doesn't depend on the FLOPs alone, but also on the kind of operation, the memory / bandwidth requirements and so on.

In fact, this question was raised in the context of GNU Radio's PFB channelizer; luckily, Tom Rondeau wrote a short article on his findings whether direct FIRs or FFT filters in the individual branches are beneficial in the environment of general purpose computers running his PFB implementation.

From that article:

Now, for GPU-based channelizers, this might be a bit different. Not an expert here, but I'd like to point you to Mr. Jan Krämer's pretty nice presentation from FOSDEM17's SDR Devroom on GPU-Enabled Polyphase Filterbanks("Everyday I'm shuffling") (slides here):

The talk highlights that the same number of ops, depending on your hardware, and memory layout, mean a very different throughput.


Estimating SNR of an Environment where obtaining clean speech is not possible is that possible to find SNR of an environment ? ie., i have a outdoor environment where i cannot find the clean speech ,if so i want to measure the signal power to noise power how to do that

can i do some thing like put a voice activity detection and find voice frames and find the voice frames power to non voice frames power will that be SNR of that environment ?

i have seen some tools like baudline giving segmental SNR will they follow the same method i have mentioned or how this actually happen i have also listened to tools like SOX

I think it all depends on the SNR region that you are dealing with. Using a voice activity detector (VAD) seems to be worth an attempt, but consider that VADs fail if the SNR is too low. The most simple approach might be using quantiles of the power distribution. Some lower quantiles will probably not contain power that is caused by the speech signal part while high quantiles will contain power resulting from the speech signal. Choosing the right quantiles might be a challenge, but due to the straight forward approach it might by worth a try...

Furthermore, by reinterpreting your problem to noise floor estimation, you might succeed by using some noise floor estimation algorithm (e.g. Minimum-Statistics by Rainer Martin [1, 2]) and calculating the SNR from that.

[1] Rainer Martin, “Spectral Subtraction Based on Minimum Statistics,” presented at the EUSIPCO, 1994.
[2] Rainer Martin, “Noise Power Spectral Density Estimation Based on Optimal Smoothing and Minimum Statistics,” IEEE Transactions on Speech and Audio Processing, vol. 9, no. 5, pp. 504–512, Jul. 2001.


differences between two L1 norm minization schemes I was reading and working with L1 regularized least squares, where:
$$ \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| \boldsymbol{x} \right\|}_{1} $$
is used to solve for sparse solutions in $\boldsymbol{x}$. However, I also stumbled on a different minimization for a similar case:
$$ \arg \min_{\boldsymbol{x}} {\left\| \boldsymbol{x} \right\|}_{1}     \mbox{ subject to }   {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{1}< \delta $$
so in the second case the expression is not unconstrained, but they switched both parts to L1 norm... What are the reasons to do so if both look for sparse solutions in  $\boldsymbol{x}$ ?

The first equation you have is often called the Quadratic Problem, which through the use of Duality can be shown to be equivalent to the Basis Pursuit De-Noising (BPDN) given as:
$$ \arg \min_{\boldsymbol{x}} {\left\| \boldsymbol{x} \right\|}_{1}     \mbox{ subject to }   {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}< \delta .$$
In your 2nd problem, the L2 norm is replaced with the L1 norm in the constraint. When you constrain the L2 error, you tend to end up with an error vector that is quite dense, i.e. a lot of non-zero elements, and each error is small with respect to $\delta$.
By switching the error to the L1 norm, then they are trying to make the error vector sparse (as opposed to dense). So in this formulation:
$$ \arg \min_{\boldsymbol{x}} {\left\| \boldsymbol{x} \right\|}_{1}     \mbox{ subject to }   {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{1}< \delta $$
They are trying to find:

An $\boldsymbol{x}$ that is sparse under the constraint that
$\boldsymbol{e}=A \boldsymbol{x} - \boldsymbol{y}$ is also sparse.

In this sense you might think that the problem they really want to solve is:
$$ \arg \min_{\boldsymbol{x}} {\left\| \boldsymbol{x} \right\|}_{0}     \mbox{ subject to }   {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{0}< \delta $$
but because that is a very difficult problem, they are using the convex formulation they have given.


Open Source Implementation of 3D Thinning algorithm I am looking for open source implementation of thinning algorithm in 3D.

  Thinning algorithm: Convert binary shapes to 1-pixel wide lines, without shortening it or breaking it apart. For example, the thresholded version of hand written or printed alphanumerics can be thinned for better represetation and further processing.

This particular thinning algorithm I am looking for does not disconnect an object or split an object up into multiple parts. Moreover, it does not shorten the shape.

Are there such source codes available? If 3D is not available, 2D is fine too. I can take it from there.

UPDATE:

Somebody has requested me to elaborate on what am I after. Here:  http://www.mip.sdu.dk/ipl98/how_to_use_ipl98/tutorial/node52.html There are several images in the link here on thinning and skeletonization. I can see that thinned algorithm and skeletonization algorithm produce similar yet different result. But what I am after is to have all shapes in a binary image to become thinned without shortening the shape. To clarify things further, here is what I wanted:

This is original image:

This is skeleton image:

This is thinned image:

As you can see, skeleton image and thinned image have different results. Thinned image has lesser lines and produce an outline that looks more like the thinned original. But these images are in 2D, I am looking for an algorithm to do thinned binary image in 3D.

If you would be happy with 3D skeleton and are looking for a Java implementation, there a popular image processing package in Life Science called ImageJ. In a distribution of ImageJ called Fiji, you will find a solid implementation by Ignacio Arganda-Carreras.

Check here: http://fiji.sc/wiki/index.php/Skeletonize3D. You will find a link to the source as well as links to tutorials.

Here is video describing how to use it from with the Fiji software: http://fiji.sc/wiki/index.php/Skeleton_analysis_video_tutorial


Can we detect multiple point target within a narrow-band transmit signal? As I understand matched filter is optimal linear filter for maximizing the Signal to noise ratio (SNR) in presence of additive stochastic noise. Consider a radar/sonar system with narrow-band transmit signal and I guess generally people use matched filter for detection, specifically for range resolution. As per my understanding, under this condition we cannot resolve if there are multiple target within the signal, my first question: is it correct? We can detect multiple target within the signal if we use linear chirp with reasonable bandwidth. Multiple target within the signal refers to point targets within transmission time, (the return comes from more than one target).

My second question (sorry not a technical question) : Is anyone in the community is interested to resolve targets within transmission time  with narrow-band transmit signal in the expense of high (O(N)=$N^3$, $N$ is number of samples) computational cost?

As I understand matched filter is optimal linear filter for maximizing the Signal to noise ratio (SNR) in presence of noise.

some restrictions on the noise apply (mainly: should be additive), but yes!

  Consider a radar/sonar system with narrow-band pulse

A pulse that is narrow band is not really much of a pulse, if you'd ask me, but that's more of an opinion but a scientific fact.

However: the smaller the bandwidth, the longer the duration of the pulse. Basically, for pulse radar, bandwidth $b$ is inversely proportional to pulse duration $T$.

  specifically for range resolution

low bandwidth is very bad for range resolution. A bit of fiddling with the radar equations will show that your resolution is bounded by the half distance your wave travels during $\frac 1b$ (monostatic rada). Now, as mentioned above, for the pulse radar, that leads to a resolution of $c\frac 1{2T}$.

  As per my understanding, under this condition we cannot resolve if there are multiple target within the pulse, my first question: is it correct?

Exactly. See statement above.

However! You can get something that is equivalent to a higher bandwidth by multiple observations under coherent conditions. I'd doubt you'll easily distinguish two targets by that without a secondary source of information, though.

  We can detect multiple target within the pulse if we use linear chirp with reasonable bandwidth

That's why people use chirps: The same bandwidth $b$ as a pulse of duration $\frac1b$, but you can spread the energy over a larger time $T\gg \frac 1b$. So, chirps, unlike radar pulses have sub-1/signal duration resolution, because their duration is not inherently linked to the bandwidth.

  Is anyone in the community is interested to resolve close target with narrow-band pulse in the expense of heavy computational cost?

So again, you do not want to go narrowband and pulse at the same time. That's just the worst of two worlds.

What you do is apply proper signal processing gain and knowledge on the target as the signal to get higher resolution than a single transmission of given bandwidth would have, and mostly, you combine a lot of observations.

For the "far far away" target, I'd encourage you to learn what a synthetic aperture radar does, in terms of signal-processing. There's pulse compression (like you know it from the chirp) and there's the combination of a lot of observations to mimic a really narrow beamwidth and get a resolution that is far higher than one could expect from a single observation.

For the closer target: OFDM Radar is your topic. You can use the existing OFDM waveforms, for example of vehicular IEEE 802.11p "WiFi", and exploit the fact that you know a lot about these and have very good estimators for channel state that you can apply on the reflections. Martin Braun wrote a dissertation on that (which, by the way, also does a very decent job of explaining typical boundaries in estimator's variance, accuracy and radar estimation possibilities in general).


MATLAB, error while trying to implement a Designed Filter object I am using the Digital Filter with properties in MATLAB. I have defined the filter object filt(1) as:

 digitalFilter with properties:
            Coefficients: [20x6 double]
   Specifications:
       FrequencyResponse: 'bandpass'
         ImpulseResponse: 'iir'
              SampleRate: 1000
      PassbandFrequency1: 59.9000
          PassbandRipple: 3
    StopbandAttenuation2: 80
    StopbandAttenuation1: 80
      PassbandFrequency2: 60.0800
             FilterOrder: 40
            DesignMethod: 'ellip'

Similarly, I have other filters under filt object as filt(2), filt(3) up to filt(8), each at different frequencies.

Now, for some reason, I get the error when my program tries to run the line:

Q = filter(filt(1),x);

where x is a 600000 x 1 column vector with only positive and negative real values (resembling a long sine wave).

The error that I get is:

Subscript indices must either be real positive integers or logicals.

Error in digitalFilter/filter (line 870)
    y = obj.FilterFunction(obj.Coefficients,x);

Could anyone please point out the problem?

EDIT:

The designfilt command opens up an interactive interface where you can change and tune your filter requirements. I have optimized the filters, and saved the results as a .mat file. (This whole program used to work, was part of a significant project which we completed succesfully, but after a year doesn't seem to work somehow).

I have attached the filter_8th.mat file. There are various filters here. I have already created these filters before hand because MATLAB takes a lot of time to create these.

In this .mat file, the filter objects are named as D_n1_50, D_n1_60, D_sig50, D_sig60 etc. (easily recognizable from the Workspace.)

The line of code that assigns these filter to the variable filt is, for example: filt = D_sig50. This is being selected in a loop. Which filter will be selected depends on a conditional if-else structure.

The .mat file containing "filters" and the variable x are in this dropbox link:

The issue has been resolved. This is due to the absence of the filter object fields:

                    FilterFunction:
                  FiltfiltFunction:
GetCoefficientsForAnalysisFunction:

In MATLAB 2015a these fields are absent from the filters objects. But in MATLAB 2016a (which is I am using, these fields seem to be required.) As the files were imported after creation in MATLAB 2015a, that seems to be the problem. Recreating, the fields in MATLAB 2016a solved the issue. The fields have also been created.

The code ran to check:

filt = D_sig50(1);
sf = struct(filt) ;
sf.FilterFunction

Created in MATLAB 2015a:
output:

                    FilterFunction: []
                  FiltfiltFunction: []
GetCoefficientsForAnalysisFunction: []

These fields in the filter object are empty (after loading these filters in MATLAB 2016a)

Created in MATLAB 2016a:

                    FilterFunction: @filterIIR
                  FiltfiltFunction: @filtfiltIIR
GetCoefficientsForAnalysisFunction: @getCoeffsForAnalysisIIR

The fields are now present in the filter object, and the code works now.


Convolution of Input signal and Impulse response I can solve the output if the input is

$$ x[n] = \delta[n+1] + \delta[n] + 2\delta[n-1] + \delta[n-2]$$

and the impulse response is

$$ h[n]= 3\delta[n] − 2 \delta[n−1] − \delta[n−2] + \delta[n−3]$$

I calculate the convolution between the strings [1 1 2 1]  and [3 -2 -1 1].

But I got a problem in computing the output when the input is
$$ x[n] = (\frac{1}{4})^{-n}  \left( \delta[n +1]  − \delta[n −4] \right)$$

and the impulse response is

$$h[n] = \delta[n] −\delta[n −5].$$

I took different strings but actually could not find what should be input and impulse response string and why?

In the second case, which caused you the trouble, you should first use the sifting property of the impulse to find out the input signal $x[n]$ and then proceed with convolution.

sifting property is:
$$ f[n] \delta[n-d] = f[d]\delta[n-d]$$

Apply this to simplify your input signal $x[n]$ as follows:

\begin{align}
x[n] &= (\frac{1}{4})^{-n} \left(\delta[n+1] - \delta[n-4] \right) \\
& = (\frac{1}{4})^{-n} \delta[n+1] - (\frac{1}{4})^{-n} \delta[n-4] \\
& = (\frac{1}{4})^{-(-1)} \delta[n+1] - (\frac{1}{4})^{-(4)} \delta[n-4] \\
x[n] & =  \frac{1}{4} \delta[n+1] - 256 ~ \delta[n-4] \\
\end{align}

So, your input sequence is $x[n] = \frac{1}{4} \delta[n+1] - 256 ~\delta[n-4] $

Convolving this with the impulse response of
$$h[n] = \delta[n] - \delta[n-5]$$ is easy, this time by utilizing the shifting property of the impulse which is:

$$ f[n] \star \delta[n-d] = f[n-d] $$

Where $\star$ is th convolution. Then apply distribution of convolution over addition to get your result as:

\begin{align}
y[n] &= x[n] \star h[n] \\
& = x[n] \star (\delta[n] - \delta[n-5] ) \\
&= x[n] - x[n-5] \\
&= (\frac{1}{4} \delta[n+1] - 256 ~\delta[n-4]) - (\frac{1}{4} \delta[n-4] - 256~ \delta[n-9]) \\
&= \frac{1}{4} \delta[n+1] - 256~ \delta[n-4] - \frac{1}{4} \delta[n-4] + 256~ \delta[n-9] \\
y[n] &= \frac{1}{4} \delta[n+1] - (256~ \delta[n-4] + \frac{1}{4}) \delta[n-4] + 256~ \delta[n-9] \\
\end{align}

I hope this helps.


Baby Cry recognition robust against the distance I am building a baby cry detection system using deep learning algorithm, Raspberry Pi and a microphone as the sound sensor.

When the baby starts crying then the system should detect the cry. It is working good when I placed the mic near to the baby and but when the mic is away I am not able to detect sound.

Please suggest me any techniques to detect sound from max 3-meter range.

To increase the efficiency of an electret microphone based sound capture system, you can :

use an electronic microphone pre-amplifier.
increase microphone directivity by a parabolic focus antenna.
use multiple microphones and apply adaptive noise cancelling.


What is the decay rate (time constant) of two One-Pole LPF's cascaded in series? If you have a One-Pole LPF, and set the cutoff frequency to $\frac{1}{2πt}$ where $t$ is the time to $\frac{1}{e}$ amplitude, it will provide an exponential decay in response to a step function (eg. going from 1 to 0) following the curve $y=e^{\frac{-x}{t}}$.
I am wondering how to similarly predict the decay rate of two One-Pole LPF's cascaded in series.
Ie. If LPF1 is filtering the signal at $t_1 = 2.4$, and then you run the output of LPF1 into LPF2 which is filtering the signal at $t_2 = 0.15$, in a steady state decay (eg. with a step from 1 to 0), what would be the effective "$t$" of this cascaded system?
I presume the cascade would still follow a nice predictable exponential decay just at a slightly different rate? If so, is there an equation that can predict the final "$t$" based on $t_1$ and $t_2$?
Thanks.

If both filters have the same time constant, i.e., if their individual impulse responses are
$$h(t)=e^{-\alpha t}u(t),\qquad \alpha>0\tag{1}$$
then the impulse response of the cascade is
$$h_{tot}(t)=te^{-\alpha t}u(t)\tag{2}$$
If both filter have different time constants with individual impulse responses given by
$$h_i(t)=e^{-\alpha_it}u(t),\qquad \alpha_i>0\tag{3}$$
the total impulse response is given by
$$h_{tot}(t)=\frac{e^{-\alpha_1t}-e^{-\alpha_2t}}{\alpha_2-\alpha_1}u(t)\tag{4}$$
I don't think it makes much sense to ask for an "effective" time constant of the cascaded systems. How should one define such a time constant? Of course, for large $t$ the impulse response $(2)$ decays as $e^{-\alpha t}$, and the impulse response $(4)$ decays as $e^{-\min(\alpha_1,\alpha_2)t}$, but that doesn't really say much about the systems' behavior.


When to use Fourier, Laplace and Z transforms? If we have an LTI system, with an input signal $x(t)$, impulse response $h(t)$ and output $y(t)$, I was under the assumption that if the input and impulse response were continuous in time, then you would use the FT on $x(t)$ and the Laplace transform on $h(t)$, and multiply them together to find the output $Y(f)$ in the frequency domain.

Likewise, if $x[n]$ was a discrete input signal, $h[n]$ was a discrete impulse response, and $y[n]$ a discrete output, then I though we would use the DFT on $x[n]$, z transform on $h[n]$, then multiply them together to get $Y(k)$ which is a discrete output in the frequency domain.

However I've seen people saying we would use the Fourier transform on the impulse response to get the transfer functions here so it's kinda thrown me off a little bit. Could anyone clarify this?

It's natural consequence of applying a transform to a convolution relation. The output $y(t)$ of an (continuous-time) LTI system is described by a convolution integral :

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

And when you apply a Fourier transform on this relation, it turns out to be a multiplication in the transform domain such as:

$$ Y(j\omega) = \mathscr{F}\{ y(t) \} = \mathscr{F}\{ h(t) \star x(t) \} $$
$$ Y(j\omega) = \mathscr{F}\{ h(t) \} \cdot \mathscr{F}\{ x(t) \}  = H(j\omega)\cdot X(j\omega) $$

Similary you can also apply a Laplace transform on it :
$$ Y(s) = \mathscr{L}\{ y(t) \} = \mathscr{L}\{ h(t) \star x(t) \} $$
$$ Y(s) = \mathscr{L}\{ h(t) \} \cdot \mathscr{L}\{ x(t) \}  = H(s)\cdot X(s) $$

And exactly the same happens for discrete-time LTI systems.


Fast Optimization for Long FIR Filters I need FIR filter lengths in the order of 1e4 and above to obtain reasonable accuracy in desired frequency response. The problem is that optimisation in MATLAB (e.g. fircls or the Optimization Toolbox) is very slow (because the problem size becomes a bit silly). Does any software exist out there that solve long FIR filter synthesis problems quickly-ish (or quicker than MATLAB)?

Or are there better ways of generating super-narrow pass band filters (suitable for FPGAs)?

Sampling will be at about 1 GHz, output should be about 10 kHz, so 1/100000 pass to stop ratio. The FIR filter is to perform weighted averaging.

Sampling will be at about 1 GHz, output should be about 10 kHz

So, you're decimating by a factor of $10^5=2^5\cdot5^5$.

This very much says that you'd normally go ahead and successively decimate.

In an FPGA, halfband filters are especially cheap and efficient to implement; better yet: they're extremely easy to design, and you don't have to change the shape the halfband filter for each step (though you certainly can, to some gain).

I'd probably (you'll have to do your own effort estimates; Sylvain Muneaut of the Osmocom project had a handy cascaded FIR decimation effort calculator, can't find it...) do the following, from a pure gut feeling

1 GHz -> Halfband
500 MHz -> same Halfband
250 MHz -> ¼-band
62.5 MHz -> simple $\frac15$-band filter
12.5 MHz -> simple $\frac15$-band filter
2.5 MHz -> good $\frac1{25}$-band filter
100 kHz -> excellent $\frac1{10}$-band filter

That way, the easiest to calculate filters (i.e. the ones with the least non-zero coefficients) run at the highest rate, and and the closer you get to the bandwidth you want to "protect", the flatter your passbands and the steeper your transition bands can get.

For complexity reasons, you could consider transferring the 62.5 MS/s or 12.5 MHz out of your FPGA into a PC and processing it there on a multicore CPU – CPU power is much cheaper than FPGA power these days, for memory intense tasks like long FIRs, at least.

Generally, the question is whether you're not doing something questionable when sampling a 10 kHz signal at 1 GHz; but that's due to the physical reality of the system you're observing, and it's hard to help without knowing that system.


Creating a flat spectrum filterbank I'm trying to create a filterbank whose frequency response is generally flat when summed in parallel. Are there any tips to choosing center frequencies and bandwidths for a filterbank of overlapping bandpass filters? Does a gammatone filterbank sum up to a flat response?

There exist perfect reconstruction filter banks that do what you suggest.

One example of the choice of $H_i(z)$ for $i= 0 \ldots N-1$ is the FFT filters (followed by the IFFT filters $F_i(z)$).


Can FFT be used to find peaks and pits on audio files I'm able to read wav files and its values. I need to find peaks and pits positions and their values.

First time, I tried to smooth it by $\frac{(i-1 + i + i + 1)}{3} $
formula, then searching on array as array[i-1] > array[i] & direction == 'up' --> pits style solution. But because of noise and other reasons of future calculations of project, I'm trying to find a better working area.

Since couple days, I'm researching FFT. As my understanding, FFT translates the audio files to series of sines and cosines. After FFT operation the given values is a0's and a1's for $a_0 + a_k * cos(k*x) + b_k * sin(k*x)$ which k++ and x++ as this picture:

My question is, does FFT helps to me find peaks and pits on audio? Does anybody has a experience for this kind of problems?

FFT is a tool to look at your signal content in frequency domain. So if you want to find peaks of your time domain signal, FFT cannot help you. However if you want to see what frequency component in the signal is the largest then you can use FFT.
To find peaks of your time domain signal you can load it into MATLAB and use command "max" for example if your signal is called "mySignal" in MATLAB type: max(mySignal) and it will give you the location of value of the peak.


How to generate colored Gaussian noise and adding it to a ODE system - Do I need Euler-Maruyama method? In the tutorial, when white noise process is added to ordinary differential equations (ODE), the ODE becomes a stochastic process. Then the stochastic process needs to be solved using Euler Maruyama method and not ODE.

I am having a hard time understanding how to generate and add colored noise in the form of process noise to a continous system such as the Rossler system. randn() is used to generate White Noise in Matlab and I directly added it to the equations. But I don't know if this is the correct way to add noise to a continous time system or not. This is what I did to add White Noise to the equations:

noise_x = randn(1);
noise_y = randn(1);
noise_z = randn(1);
y = [
     -x(2)-x(3) + noise_x;            % dx/dt
     x(1)+a*x(2) + noise_y;           % dy/dt
     b+x(3)*(x(1)-c) + noise_z;       % dz/dt
] ;

color noise process using Euler-Maruyama Method (EM) seems to be the way to solve the stochastic system after the noise is added. But I cannot understand what is the color of the noise from this method. Looking at the code, there is no way to know if there is any color. After the color noise is generated by calling the function function [W,t]= colornoise(a,b,T,N) how to use the output W as an additive term in the Rossler system. This is not clear to me. Can somebody please provide a small code snippet on how to add White noise and color noise such as pink, red, violet to the Rossler system.

The following are my questions:

Question 1) How to generate colored noise in general?

Question 2) How to add the colored noise in the form of process noise to the Rossler system or any continuos time system and which method to use for solving the system -- ODE or Euler Maruyama?  Is the noise added Rossler system stochastic or not?

The code below is for integrating the Rossler system using ODE45 solver.

clc;
clear all
%%%% Number of variable and initial conditions:

nbvar=3;
xini=ones(1,nbvar)/10;

%%%% Time parameters:

trans=100;
tend=500;
tstep=0.01;

b=0.2;   % (default value for chaos)

ttrans = [0:tstep:trans];
tspan = [0:tstep:tend];

x0=xini;
option = [];
[t x] = ode45(@dxdt,tspan,x0,option,b);

plot3(x(:,1),x(:,2),x(:,3));
xlabel('X','fontsize',18);
ylabel('Y','fontsize',18);
zlabel('Z','fontsize',18);
box on;

% ===================================================================
% dxdt
% ===================================================================

function y = dxdt(t,x,b)

%%% parameters

a=0.2;
%b=0.2;
c=5.7;

%%% equations

y = [
     -x(2)-x(3);            % dx/dt
     x(1)+a*x(2);           % dy/dt
     b+x(3)*(x(1)-c);       % dz/dt
] ;

Question 1) How to generate colored noise in general?

Produce white noise and filter it.

  Question 2) How to add the colored noise in the form of process noise to the Rossler system or any continuos time system

This is a little bit more complicated. The addition of noise is a simple matter of $y = x + q$ where $q$ is the noise signal. The key problem here is how much noise to add? This is important for linear cases (such as $y$ previously) because random number generators are producing numbers in some fixed interval ($[0 \ldots 1]$, $[-1 \ldots 1]$ or other) which will be entirely different than the interval the signal $x$ occupies.

In this case, assuming that the signal is stationary, then given some $x$, generate a $q$, scale $q$ so that it has the same power as $x$ and then mix them together at the desired ratio. For example $y = (1 - a) \cdot x + a \cdot q$. Here, the contribution of noise is directly proportional to $a$.

This is nice because you can now use $a$ later on in your experiments as a variable that shows how much noise was added to the signal.

But, when you are dealing with non-linear systems, this becomes even more important (and difficult). The reason for this is that, depending on the system, you cannot assume that it is stationary. Or, you can, but only at a narrow range of operation and furthermore, you are going to have to make sure that the system stays in this region.

The added difficulty is in two places:

1) You are dealing with differential equations that get integrated

Therefore small spurious excitations here and there can accumulate to large numbers quickly

2) You are dealing with non-linear systems

Therefore, combined with #1 it becomes likely that you might push the system into a region that its behaviour changes. For example, you might be interested to use the system when it produces nice smooth waveforms but inadvertently you might tip it over to the region where it starts producing spiking oscillations.

The reason I am mentioning this is because you refer to "process noise". So, you can take the Rossler model and integrate it and look at its signature attractor as an experiment or, you might want to use it as an actual model of a process, in which case, $x,y,z$ now have units, $dt$ is tied to your sampling frequency or other physical limit and also your noise now has a physical meaning too.

I don't know what these are in your case, but you have to understand them in order to be able to add the noise, especially if you have to add the noise in the derivatives (are you sure about this?).

The derivative is not the signal but rather the rate of change of the signal with respect to time. This is incredibly powerful. If you were to produce some $x$ by integrating the $\frac{dx}{dt} = b$, you would get a ramp in $x$. If your derivative now was a ramp itself with $\frac{dx}{dt} = c \cdot t + b$, your $x$ would vary exponentially.

Imagine then, how much you are impacting $x$ by specifying some random factor in its $\frac{dx}{dt}$.

If you don't have to add the noise at the derivatives, it would be slightly easier, because you simply integrate your differential equations to produce some output $z$ and then add noise to this specific $z$ as if it was our $x$ from above.

  ...and which method to use for solving the system -- ODE or Euler Maruyama? Is the noise added Rossler system stochastic or not?

If you have to add the noise to the derivatives then the system becomes stochastic and it would be better to use the Euler - Maruyama. The difference of Euler - Maruyama from the plain Euler integration is in the inclusion of a derivative term for the random process as well. Of course, if you differntiate white noise, you still get white noise, so from a practical point of view, there is no difference there because you still just have to add a random number. BUT this random number has to be appropriately scaled for its amplitude and then be multiplied by the integrating step too.

Hope this helps.


Stable digital filter with leading zeroes I'm looking construct a stable pole-only filter where the feedback coefficients start with a block of zeroes, i.e.

\begin{align}
a_0 &= 1\\
a_i &= 0, \textrm{ for}\ 1 \le i \lt k\\
a_i &\in\mathbb{R}, \textrm{ for}\ k \le i \le n
\end{align}

For $k = n$, this is a single-tap feedback delay, e.g.

$$
y_i = x_i + a_k y_{i - k}
$$

which is stable as long as $|a_k| < 1$.

However, I'm looking for a way to construct such filters with arbitrary $k$ and $n$.  Are there conditions or constraints that would let me construct this?

OK, so this is a bit of an intuitive argument so I'm not 100% sure about it, and it's totally sideways to the way I normally think about filters, but:

Let's consider the "feedback kernel" ($a_k ... a_n$ - so, not including $a_0$) as an FIR filter. If the magnitude of that kernel in the frequency domain is $< 1$ including inter-bin peaks, the original filter should be stable.

If so: to construct a filter: one could generate a random set of $a_k ... a_n$, find the maximum freq response of that feedback kernel, and scale such that the peak is less than 1, to get a stable feedback kernel.

To reason about this: let's construct an infinite sequence of sequences, $Y^k$:

$$
Y^0_i = a_0x_i\\
Y^j_i = -\sum_{l=k}^na_lY^{j - 1}_{i - l}
$$

That is, each $Y^j$ is the previous $Y^{j - 1}$ convolved with our feedback kernel (the negation doesn't affect the argument).  I also think it's true that:

$$y_i = \sum_j Y^j_i$$

Let's say $x_i$ is "time-bounded" (zero outside a finite range) - it follows inductively that every $Y^j$ is also time-bounded, because it's a convolution of the previous.

Let's also say that the feedback kernel's transfer function has magnitude $< G$ for all frequencies.  This guarantees that whatever the input, the total energy in each $Y^k$ is strictly less than $G$ times the previous one.  This implies a geometrically-decreasing upper bound on the magnitude of the sample values of $Y^j$ as well.

This geometric progression of bounds means that when we sum up $y_i$ it converges if $G < 1$, and to a value bounded by $K/(1 - G)$ for some $K$.  The time-boundedness means that it decays over time (hand-waving over the formalism here).

Does this hold up?


Speech segmentation for speaker recognition I'm trying to implement a speaker recognition system and want to make sure I'm aware of the latest trends in speech segmentation. I've read a number of very different methods at a higher level but I'm not sure which one would be best for my particular situation. I'll list what I found so far and I'd appreciate either (a) comments on the research I've done so far or (b) popular methods that I didn't mention that I should know about. Then I'll take it from there in terms of determining what will work best for my situation.

Techniques for Speech Segmentation
(that I've identified so far)

I've read that strong long-term modulation frequencies between 1-16 Hz are indicative of speech activity, but haven't been able to find any good explanations of what exactly these are. What I've been able to determine so far is that they are the temporal spectrum of specific frequency bands. (Unsupervised speech/non-speech detection for automatic speech recognition in meeting rooms by Maganti et al)
Using energy variance (high variance implies speech-like segments) independently or with other methods to build a Gaussian mixture model with two components: one for speech and one for noise, and then re-classifying windows as either speech or non-speech. (An unsupervised, sequential learning algorithm for the segmentation of speech waveforms with multiple speakers by Siu et al)
Using the Voting Experts algorithm, which utilizes iter- and intra- window entropy (a window in this sense is a 'window of windows') to determine where "chunks" are using a feature that has high inter-segment entropy, and then to classify these chunks as either speech or non-speech using 2-means clustering, where each cluster of chunks contains all chunks of one of the two classes. (Voting Experts: an unsupervised algorithm for segmenting sequences Cohen et al)

The short version of my question

Ultimately, I'd really like to know what is standard in the industry as of right now.  Then I can determine what will work best for my particular situation once I have a little more base knowledge.

I'm going to answer my own question, for those who are looking for the same information as I was. Adaptive Thresholding is the most common way to segment audio, and it gives you the most bang for the buck. While there are more accurate methods out there for segmenting, adaptive thresholding is fast and very accurate, and uses simple features such as zero-crossings, energy, and entropy. I'm going to have to look into this more myself but this should be enough to get started.


efficient triangle to sine morphing? I'm coding a synthesis external for MaxMSP in C.
It's based on Peter Blasser's bounds and bounce concept ( http://petermopar.blogspot.co.uk/2014_04_01_archive.html )- basically triangle waves with constant cursor rate (rather than frequency) and modulated bounds - eg instead of moving between -1 and 1, the wave moves between lo and hi, with lo and hi modulated.

This morning I've been thinking of adding waveshape morphing between triangle and sine (or at least sinusoidal since the moving bounds will deform the sine).

I can think of a fairly simple way to do this:

// Assuming classic triangle moving between bounds -1 to 1
// and defining position of triangle cursor @ n = pos
// shape_amt between circa 0.1 (Tri) and 1.0 (sine)
out = sin(pos * shape_amt * 0.5* pi ) / sin(shape_amt *0.5 * pi())

ie - basic waveshaping technique - scale the portion of the lookup table that's being scanned (in this case some segment of half a sine wave) to vary the shape.

For my specific case all I'll need to do from there is build in scaling from relative position between to bounds to relative position between -1 and 1.

My question really is- is this a horribly inefficient way of achieving what I'm after? I'm new to DSP coding and to coding in C so I'm sure the way that's most obvious to me isn't the best.

The obvious performance improvement that occurs to me is to create a wavetable of the sin function and look up rather than calculate those values. Obviously that way I lose resolution on the "straighter" portions of the wave, but either a high res lookup or decent interpolation feels like it ought to deal with that fairly well. Mainly I'm considering this for LFO applications so quality needn't be utterly stellar (though it'll be nice nice if it's useable into the audio range without too much fuzzy aliasing).

But beyond that, is there just a completely different approach that would be more efficient?
And if this is the best approach, are there any things I should think of when pursuing it?

Edit - Thanks for great answers below. I've ended up doing something fairly simple (and possibly crude) for the moment, which suits my purposes pretty well. I'll definitely be investigating both suggestions as improvements to this.

#define sign(a) ( ( (a) < 0 )  ?  -1   : ( (a) > 0 ) )
#define LKTBL_LNGTH 2048

void setup_lktables (t_object* x)
{   // create lookup for 1/4 sine
    int i;
    for(i=0, i< LKTBL_LNGTH, i++){
        x->sin[i] = sin(pi * i * 0.5 / LKTBL_LNGTH) ;
    }
}

double lookup_shape (t_double pos, t_double shape, t_double lo, t_double hi)
// acting on triangle wave (phase + (2*gradient*f/sr) per sample, wraps & changes dir  at <lo and >hi)
// shape comes in as restricted to (-1...-0.05, 0.05 ...1), defines the portion of lookup to use

{
    t_double midpoint, halfwidth, ph, fracph, shaped;
    t_int maxph, intph, sign;
    // get relative position between bounds for waveshaping lookup
    midpoint = lo + 0.5f * (hi - lo);
    halfwidth = midpoint - lo;
    // prepare phase values for lookups
    maxph = shape * LKTBL_LNGTH - 1;
    ph = (pos - midpoint) * maxph /  halfwidth;
    sign = sign(ph);
    ph = ph * sign;
    intph = (int)ph;
    fracph = ph - intph;
    // lookup, scale & lerp
    shaped = sign * (x->sin[intph] * (1.f - fracph) + x->sin[intph+1] * fracph) / x->sin[maxph];

    //now return  waveshaping output scaled to actual bounds
    return midpoint + shaped * halfwidth
}

This answer is about using sines instead of triangles in the bounds and bounce concept. There is no way to do it so that these defining conditions of the original concept are all satisfied:

While the boundaries stay constant the output is a pure sine wave (triangle for the original algorithm) that goes from the upper boundary to the lower boundary.
Moving the boundaries suddenly does not create a discontinuity in the output.
Moving the boundaries suddenly does not create a discontinuity in the derivative (slope) of the output.

Say the sine wave has reached the upper boundary and the boundaries are updated. Sine waves of other amplitudes cannot have zero derivative at that amplitude value, so the scheme fails.

But you can make it work if you constrain the boundary changes to times of zero crossings of the sine wave:

You just have to decide on the slope at zero crossing and your sine turns out as:  $$\text{boundary}\times\sin\left(\frac{\text{slope}\times\text{time}}{\text{boundary}}\right).$$


Real Fast Fourier Transformation (FFT) in 2D changes with changing the axes order I am using Python to calculate the real FFT for a 2D array. I found that the real FFT function does not return an array with the same size as the input array, rather, it returns an array with the same size as the first axes and half of the size of the second axes or vise verse based on the axes order that we choose.

I attached the complex FFT, real FFT with axes order (0,1), and real FFT with axes order (1,0) of an array with size(1400x144). See the difference between the second and the third image. Any idea why that happens, and what does it mean in terms of the linear property of the FFT. Thanks

A real, 1D DFT is fully described by the non-negative frequencies, since real DFTs are symmetric about the 0 frequency bin.

The function you linked is only computing the FFT along one dimension, and is only returning the non-negative frequency bins.

From the manual page:

"Notice how the final element of the fft output is the complex conjugate of the second element, for real input. For rfft, this symmetry is exploited to compute only the non-negative frequency terms."

You can use symmetry to fill in the other half of the result, if you need it.


How to design a filter with custom group delay -- why is it an IIR? Say I have an FIR filter with desired magnitude characteristics. I could imagine using an allpass IIR in series with the FIR to alter the group delay characteristic to provide both a 'magnitude equalization' and 'delay equalization' curve custom for the application.

I imagine this would be simpler in the frequency domain where the imaginary component of each slice of the FFT would adjust the group delay of a bandwidth corresponding to the nyquist frequency of the system divided by the FFT length, and then IFFT'ing to find the impulse response, windowing that impulse response to design the filter.

I see Matlab has an IIR feature for custom group delay in fdatool. Why does it have to be an IIR? Can an FIR have a custom group delay?

Thanks for any pointers!

As pointed out by Peter K., it is true that many well-known techniques for designing FIR filters actually only design linear phase filters. However, FIR filters are very well suited for delay equalization, simply because the design process is much simpler than for IIR filters. The reason for this is the fact that the design problem can be formulated in such a way that it is still linear in the unknown filter coefficients, which is usually not the case for IIR filters. This type of design problem is known as complex approximation problem, because in the case of a desired non-linear phase response a complex frequency response must be approximated by the filter, unlike a real one in the case of a linear phase filter. You can find a lot about this topic, including Matlab programs, in this thesis. Have a look at the design example #2 on page 119 ff. where an FIR magnitude and phase equalizer is designed compensating the effects of an analog anti-aliasing filter. The following figure is taken from the thesis and it shows the group delay of the analog low pass filter, the group delay of the equalizer(s) (there are two with different magnitude responses), and the total group delay of the cascade of analog and digital filter, which is supposed to be constant in the passband:


PLL and zero crossing detection Is there an algorithm which can give the frequency and phase of a three phase signal (three phase voltages) without using zero crossing detection.

I mean an algorithm which is done in open loop.

What about this: if your three signals are balanced (I mean with equal amplitudes and $2\pi/3$ shifts), you can find the amplitude by summing the squares (check that it doesn't pulse). And from there, the arguments of the three sinusoids ($\omega t+\phi+k2\pi/3$). Actually you solve this geometric problem, where the three projections are given:

Then, repeating the computation at two close instants, you get the frequency. It may be useful to maintain a linear fit on the last few $n$ phase values to improve accuracy and ease handling of phase wraparound.


Relationship between windowed fourier transform and wavelet transform I was reading on windowed fourier transform and wavelet transform, and i was thinking that the windowed fourier transform is a subset of wavelet transform. Is that true?

Wavelets commonly factor the frequency of the basis vectors into their temporal widths, some with more than one temporal location.  Windowed FFTs fix the temporal width of each basis vector to the length of the FFT, thus independent of any bin frequency, and with the basis location fixed to the center of the FFT.

I suppose one could consider the FFT to be similar to a wavelet transform with a bad set of wavelets (too many high frequency ones and with poor temporal/locality resolution, and too few low ones and with poor frequency resolution).


Relating transfer functions with step responses Relate the transfer function to its' corresponding step response.

First, I tried setting up the poles and zeros of the transfer functions. This helped a bit since I know that $G_A (s)$, $G_B(s)$ and $G_C (s)$ should all have oscillatory step responses (because of conjugate poles). The problem arose when trying to determine the differences between $G_A (s)$ and $G_B (s)$, and we thought that there has to be a better way of understanding this. Working with poles and zeros, we didn't really understand what we were doing.

So, I'm now thinking that I need to use the fact that the inverse laplace of the transfer function equals the impulse response, and the integral of the impulse response equals the step response. This is how we relate the step response with the transfer function, but we don't know how to work with it.

Edit:
For example, when I take the inverse laplace of $G_D (s)$, I get $\frac{1}{5}$, and the integral of this is $t/5$, which can be translated to the first curve above.

But when I take the inverse laplace of $G_F (s)$, I get $9te^{-3t}$, and the integral of this is $-9e^{-3t}(3t+1)$, but this is too much of a hassle, how can I relate this to one of the step responses in an easy way?

Maybe I should go back to using the poles and zeros, any tips would be greatly appreciated.

In the given example you have 3 types of systems:

the ideal integrator (D) with a step response that grows infinitely (fig. 1); that's the obvious one, as you've found out by yourself.
underdamped systems with complex conjugate poles (A, B, and C): their step response oscillates
critically damped systems with a double real pole (E and F): no oscillation in the step response

So systems A, B, and C must correspond to figures 2, 4, and 5, whereas systems E and F must correspond to figures 3 and 6. You had already arrived at that conclusion.

But we have more cues. Clearly, system C has less damping than systems A and B, which have the same damping. Smaller damping corresponds to more oscillation (i.e. higher overshoot and undershoot). So we can conclude that system C corresponds to figure 4. The difference between systems A and B can be seen by using the Initial Value Theorem:

$$f(0^+)=\lim_{s\rightarrow \infty}sF(s)\tag{1}$$

For system A we get

$$g_A(0^+)=1$$

whereas for system B we have

$$g_B(0^+)=-1$$

These are the values of the impulse responses, i.e., the derivatives of the step response. We see that the step response in figure 5 has a negative derivative at $t=0^+$, from which we can conclude that it corresponds to system B. Then we're left with figure 2 for system A.

The difference between systems E and F is the distance of the double pole from the imaginary axis. The further away from the imaginary axis, the faster will the impulse response decay, and consequently, the faster will the step response rise. So system E must correspond to figure 6, and system F corresponds to figure 3.