Spectrogram with square or non-square magnitude of STFT: power vs. magnitude As seen in this question and answer, to do a spectrogram, it's common to plot either:

the square magnitude $|\text{STFT}(\text{frame}, \text{bin})|^2$ ("power spectrum")
the magnitude $|\text{STFT}(\text{frame}, \text{bin})|$

In an audio context, when to choose the one or the other? For which reason is one more useful than the other?

In the context of noise reduction, there are both Magnitude Spectrum Subtraction or Power Spectrum Subtraction, as seen around page 5-8 in http://dsp-book.narod.ru/304.pdf. When to use the one or the other?

Also, what do most audio editors display? spectrum or power spectrum?

In the image you gave, the right-hand scale is in decibels (dB). So, essentially, the square turns into an affine scaling in the logarithm domain, which essentially yields the same image, at least the same relative dynamic range. So, in that case, it does not matter.

Aside, $x \to x^\alpha$ transformations with $x\ge 0$ (a positive signal, a magnitude spectrum, a spectrogram) are customary when signals have different orders of magnitude that cannot be caught by human eyes or ears: $\alpha > 1$ puts emphasis on the largest components,   $\alpha < 1$ relatively enhances smaller values, which is useful when you want to detect low-amplitude signals.
Depending on whether you want to capture audio fingerprints, or perform subtle audio forensics, one or the other should be chosen. Such $\alpha$-power transformation are also useful in your other spectral subtraction question, What is the name of this very simple spectral subtraction technique?


Smallest sampling frequency to fully reconstruct a signal Given the spectrum of analog signal $x_a(t)$ which is imaginary and band-limited, find the lowest sampling frequency to be able to reconstruct $x_a(t)$ from samples $x[n]$.

My attempt: Bandwidth of the signal is $\Delta \Omega = \Omega_2 - \Omega_1$
Such the sampling frequency should be $\Omega_T = 2\Delta \Omega$

Your understanding is not correct for complex sampling:

For the specific case of a complex bandpass (analytic) signal whose spectrum is by definition zero for the negative frequencies such as:

$$ H(\Omega) =  \begin{cases}  0 &, \text{ by def. for } \Omega < 0 \\ A_\Omega &, \text{ nonzero for } \Omega_1 < \Omega < \Omega_2 \\ \end{cases}$$

then the necessary minimum sampling rate which would avoid aliasing (spectral overlap) is given by $$ \Omega_s > ( \Omega_2 - \Omega_1 ) $$ Note that if the signal were real with a symmetric bandwidth then the minimum sampling rate would be twice that of the complex case. Also then the allowed range of valid sampling rates would be found differently.

The reconstrcution filter that would yield the original complex signal $x(t)$ back will be a complex bandpass filter, not a real one. You can obtain the complex bandpass filter $H_+(\Omega)$ from the real bandpass filter $H_r(\Omega)$ using a Hilbert transformer such as:

$$ h_+(t) = h_r(t) + j \mathcal{H} \{h_r(t) \} $$

Note that for a given sampling rate you will be getting $F_s$ complex samples per second which is equivalent to $2 \times  F_s$ real samples per second. Hence the apparent advantage of complex sampling is actually not realized , as the total number of samples per second will be the same in both cases.


Regularized Least Squares by Laplacian Operator - Optimal Value of the Regularization Factor (Lagrangian Multiplier) Consider the cost function

$$f(X,\lambda) = \|AX-b\|_2^2 + \alpha \|LX\|_2^2$$

$A:$Measurement matrix($R_{m\times n}$,$m \ll n$), $b:$observation vector($R_m$), $L:$Laplacian operator($R_{n \times n}$), $X:$vector form of an image($R_n$)

(The above cost function can be read, non-mathematically as minimize

$AX=b+error$ subject to $LX=0$)

$$\frac{\partial f}{dX} = 0 \implies 2X^T(A^TA) - 2X^T(A^Tb)+\lambda(L^TL) = 0.........(1)$$

$$\frac{\partial f}{d\lambda} = 0 \implies \|LX\|_2^2=0.........(2)$$

From $(1)$ $X_s = (A^TA+\lambda L^TL)^{-1}A^Tb.........(3)$

$(3)$ in $(2) \implies \|LX_s\|^2_2 = 0..........(4)$

It is very clear from $(3),$ $\lambda$ plays a very crucial role in solving for $X$ because changing the value of $\lambda$ gives a different estimate for $X$.

using $(4),$ is it not possible to auto compute the optimal value of $\lambda$?

If it is possible to auto compute optimal $\lambda,$ why is it in many image de-noising problems we specify the optimal lambda $\lambda$ as an input to the algorithm?

It it is not possible to auto compute $\lambda,$ why is it so?

Thanks for reading my question with patience..

Let me rewrite the problem:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| L x \right\|}_{2}^{2} $$

In order ot say something intelligent about the optimal value of $ \lambda $ one must set an optimization criteria.
Since you don't have one there is no well defined answer.

What might assist you is reformulating the problem in an equivalent form:

$$\begin{align*}
\arg \min_{x} \quad & \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} \\
\text{subject to} \quad & {\left\| L x \right\|}_{2}^{2} \leq \alpha &
\end{align*}$$

Both forms are equivalent in the sense that for any $ \alpha $ of the 2nd form there exists $ \lambda $ such that both solutions are the same.

Now if you can say something about the smoothness level of the image which will let you set $ \alpha $ intelligently than you'll be able to extract $ \lambda $ to match it.

Pay attention that the connection isn't trivial and depends on the data itself among other things (See my answer to Significance of λ in Basis Pursuit).


How to change the measured curve to be close to the target curve I am new to Signal processing and want to implement for the measured curve from my experiment.
I have a measured curve (like skewed gaussian).
I like to change that measured curve to be as close as possible to the target curve (gaussian) (statistically and visually). Visually means, the target curve's width is abit longer than the measured curve when we the two curves are overlapped. Which signal processing should I use?
I like to attach a picture, but it is not allowed for the new person.
Thanks
thomas

It is unclear what is your application from your question.

If you aim at "gaussianizing" your data -- for example to make better use of statistical modeling, regression, or machine learning methods which work better under a gaussian assumption -- you can use the Box-Cox transform. It consist in applying the following transformation to your data:

$$y = \frac{x^\lambda - 1}{\lambda}$$

$\lambda$ is chosen as to maximize the log-likelihood of a gaussian model given the transformed data - this is a maximization problem usually solved using Levenberg–Marquardt.

I would also strongly suggest posting this to the stats SE.


Is DCT (Discrete Cosine Transform) of Type-2 lossless or lossy? As I know, in general, the DCT is lossless. But I'm not exactly sure about Type-2 of DCT.
Is it lossless or lossy?

Mostly yes, but it depends on the context. Let us elaborate.

DCT-II is one of the many forms of Discrete Cosine Transforms, and probably the  most widely used one, as it is (somehow) present in JPEG or MP3 formats. "Lossy" often refer to the compression standard which uses it, because the  main loss results from quantization (and generally not the transform by  itself). By itself, the transform is invertible, even more orthogonal or orthonormal, so in theory you have no loss (it is bijective).

As @Fat32 points out, its coefficients sometimes involve cosines, hence irrational numbers (except at specific values) that are not easily  represented in finite arithmetic (float or int), and can induce round-off errors.

But:

DCT-II of size $1$ is $1$, DCT-II of size 2 (unnormalized) is: \begin{smallmatrix} 1 & 1 \\1 & -1\end{smallmatrix} so DCT-II can be lossless (in special cases though),
For integer data (like images), it seems possible, with enough bit-depth computations (related to the initial integer range and the DCT size), to remain practically lossless, as you know that the  transform/inverse transform result should be integer,
Very  accurate integer, or dyadic approximations exist, meant for integer hardware: for instance the binDCT.

My final answer is thus: DCT-II is mostly lossless, if you take care of it.


Good inhomogenity measures in images? Can you advise me on good inhomogenity measurements in an image?
I am trying to find the borders where homogenous zones touch inhomogenous ones in a image.

I am going to work with sick lung ct images like this:

I am looking for different ways to measure inhomogenity in the texture in the image.

Thank you

The variance is a good inhomogenity measure.
For example, take 20*20 patches and calculate the empirical variance of the patch using these 400 samples.
For very simple purpose, this could work.

Edit : After viewing your example images, try to just play around with simple tools like manual thresolding and edge detection.

Edit 2: Added Example with different patch sizes (40, 10, 2)


Design a linear-phase FIR filter approximating the magnitude of a given IIR filter I have a biquad IIR filter from which I want to get a linear-phase FIR. I saw that related question but the OP does not care much about phase.

I have troubles getting the exact same magnitude curve when extracting an FIR from my IIR.

So far, my method has been:

Apply my IIR filter on a 1024 samples buffer which contains nothing but a Dirac as first sample.
Shift the impulse response obtained in step 1 and make it symmetric (to get phase linearity)

If I stop after step 1, I get the exact same magnitude curve for my IIR and my FIR but I also get the exact same phase curve (which is non linear, and hence not interesting).
If I stop after step 2, I get an linear phase but not the exact same magnitude.

On the following picture, IIR's and FIR's (after step 2) transfer functions are (respectively) in red and blue:

What have I been doing wrong?

What you do in step 1 is simply truncate the infinite impulse response to approximate it by an FIR filter. If you use sufficiently many filter taps, the approximation becomes arbitrarily accurate. This means that the resulting FIR filter approximates the magnitude and the phase characteristic of the original IIR filter. So with this approach the phase will never become linear.

Making the impulse response symmetric to obtain phase linearity, as you do in step 2, will of course change the magnitude response.

What you should do is use the magnitude of the IIR filter as a desired response in a (linear-phase) FIR filter design routine. In that case you will get an FIR filter with an exactly linear phase and with a certain magnitude approximation error. That magnitude error can be made sufficiently small by choosing an appropriate filter order. The simplest approach is probably to use a least squares approximation, which just involves solving a system of linear equations.

Example: I use a peaking EQ filter as the IIR prototype. The coefficients are (b are the numerator coefficients, a are the denominator coefficients):

b = [1.2223e+00, 0, 7.7775e-01];
a = [1.1250e+00, 0, 8.7502e-01];

You can use the magnitude of the IIR filter's frequency response and combine it with a linear phase to obtain the desired response for the FIR filter design routine (N is the filter length). The code is Matlab/Octave syntax:

[H,w] = freqz(b,a,256);
N = 61;
D = abs(H).*exp(-1i*w*(N-1)/2);

You can use a least squares FIR filter design routine called lslevin.m, which you can find here.

h = lslevin(N,w,D,ones(length(w),1));
Hh = freqz(h,1,256);

The figure below shows the magnitudes of the two frequency responses (IIR and FIR):


Conditions for finite bandwidth I'd like to know which conditions must be satisfied in order to get a signal with a finite bandwidth.

Thank you for your time.

Let us  stick to traditional functions, considering the almost everywhere equivalence class.  From an actual signal perspective, I am leaving distributions aside.

According to the Paley-Wiener theorem, functions with compact frequency support (corresponding to Bernstein spaces) can be extended
from the real line to the complex plane $\mathbb{C}$ as entire functions
of exponential type; entire functions are holomorphic at all finite points of $\mathbb{C}$ ,  exponential type means that $|f(z)|\le C e^{a|z|}$. They can be called $B$-functions.

So they inherit from all their properties: no bounded support in time, infinite differentiability, Taylor expansion everywhere, controlled decay at infinity...


How to Apply Statistical Algorithms of Signal Processing to Regulate Variation of a Curve? Below I am posting 2 graphs. I want to regulate the curvature of first graph using some statistical methods such as use of standard deviations, and modulate my graph to look like second one.

I am not asking here for any matlab command, but for some statistical method or a signal processing algorithm I can implement using MATLAB.

I just want transformation of shape as that of image 2, I do not want to achieve exact values of image 2.

Thanks & Regards, Nick.

.

You can regulate its second derivative which is the curvature.

Something like:

$$ \hat{x} = \arg \min_{x} \frac{1}{2} \left\| x - y \right\|_{2}^{2} + \frac{\lambda}{2} \left\| D D x \right\|_{2}^{2}  $$

Where $ y $ is samples you have and $ D $ is the Derivative Operator (In Matrix Form). By applying it twice we're regulating the Second Derivative.
You can find a statistical explanation for this (Prior of Gaussian Distribution of the Second Derivative) or just call is Tikhonov Regularization.

This has a closed form solution:

$$ \hat{x} = \left( I + \lambda {D}^{T} {D}^{T} D D \right)^{-1} y $$


Possible to get transfer function coefficients from window? I am hoping to use scipy.signals.filtfilt() to smooth some signals in Python, and wanted to build the filter based on a window like a hanning window or whatever.  E.g.:
import scipy.signal.windows as windows
window = windows.hann(filter_width)

But standard filters don't just take in windows, they take in numerator and denominator transfer function coefficient arrays a and b:
data_smoothed = scipy.signal.filtfilt(b, a, data_noisy)

Is there a way to calculate the transfer function coefficients a and b from a window? I like filtfilt() more than straight-up convolution with the window because it has a lot of useful features baked in.

What you are describing is an FIR filter, such that all the denominator coefficients are zero, save the basis, a[0]=1.  So you could do something like:
data_smoothed = scipy.signal.filtfilt(window, 1, data_noisy)

There is a notable point.  The DC gain of the filter is equal to the sum of the coefficients for FIR filters.  Your window is likely normalized to 1, so the sum is probably higher than 1, which means your filter will have gain at low frequencies.  You would want to divide all the coefficients by the sum of the window to keep the gain to unity.


Estimating the PSD (Periodogram): How do you find the Expected Value The periodogram $\widehat{P_x}(e^{jω})$ is an estimate of PSD $P_x(e^{jω})$ and for mean-square convergence -

$\lim_{N \rightarrow \infty}E\left\{\widehat{P_x}(e^{jω})\right\}=P_x(e^{jω})$

My question:

If I calculate the expected value, which is according to definition an average, this way

$E\left\{\widehat{P_x}(e^{jω})\right\} = \frac{1}{2π}\int_{<2π>} \widehat{P_x}(e^{jω})dω$

Then, $E\left\{\widehat{P_x}(e^{jω})\right\}$ should have a single value. If this was true, then it would mean that $P_x(e^{jω})$ has the same value for any $ω$ and this clearly is not true. So, where am I mistaking?

TLDR;

The mistake is here:

  If I calculate the expected value, which is according to
  definition an average, this way

  $E\left\{\widehat{P_x}(e^{jω})\right\} = \frac{1}{2π}\int_{<2π>}
> \widehat{P_x}(e^{jω})dω$

The mistake is in your interpretation of expectation. It is not an average over all frequencies.

For each fixed $\omega$, the expectation is taken over many realizations of the random signal $x$.

Gory details

Suppose you are given a complex valued random signal $x(n)$ for $0\leq n\leq N-1$. I would emphasize that this is a random signal. If you repeat your experiment you'll obtain a different set of $N$ values each time (called "realizations" of $x$).

Recall how the periodogram $\hat P_x(e^{j\omega})$ is computed. You first compute the sample autocorrelation at each lag $m$:
$$
\hat r_x(m) = \frac{1}{N-m}\sum_{n=0}^{N-m-1}x(n)x^*(n+m) \;\;\;\;\;\;\;\;\;\;\;\;(1)
$$
and then get the periodogram estimate by taking the DTFT:
$$
\hat P_x(e^{j\omega}) = \sum_{m=-(N-1)}^{N-1} \hat r_x(m) e^{j \omega m}.
$$
Since we started with a random $x$ we would like to know how this estimate performs on average i.e. average over the randomness in $x$:. So we take the expectation:
\begin{eqnarray}
\hat P_x(e^{j\omega}) &=& \mathbf E\left[\sum_{m=-(N-1)}^{N-1} \hat r_x(m) e^{j \omega m}\right]\\
&=& \sum_{m=-(N-1)}^{N-1} \mathbf E[ \hat r_x(m)] e^{j \omega m}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(2)\\
&=& \sum_{m=-(N-1)}^{N-1} \mathbf E\left[ \frac{1}{N-m}\sum_{n=0}^{N-m-1}x(n)x^*(n+m) \right ] e^{j \omega m} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(3)\\
&=& \sum_{m=-(N-1)}^{N-1} \frac{1}{N-m}\sum_{n=0}^{N-m-1}\underbrace{\mathbf E[x(n)x^*(n+m)]}_{\mbox{autocorrelation function of }x} \; e^{j \omega m}\\
&=& \sum_{m=-(N-1)}^{N-1} r_x(m) e^{j\omega m}  \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(4) \\
&=& P_x(e^{j\omega})    \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(5)
\end{eqnarray}

where (3) follows by substituting (1) in (2), (4) is by definition of autocorrelation $r_x$ of the random sequence $x$ and (5) is by definition of $P_x$.

Once again, note that the expectation was taken to "average out" the randomness introduced by $x$. It is not an average over frequencies because the frequencies are not random here.


how can i perform a registration to 2D images using ''rigid " transformation in matlab? i need to perform registration to 2D image using"rigid transformation" in  matlab .i need this for measuring the similarity between my images .any one know how this can be done using matlab?

an example on my images :

regards

A rigid, or isometrie transformation, is a transformation in which the distances between the points that are transformed, are preserved. This means that an object formed by these points will have the same shape and size after transformation. Examples are rotations and translations.

A 2D column vector can be rotated by $\theta$ as follows:

$\left(\begin{array}{c}x' \\ y'\end{array}\right) = \left ( \begin{array}{cc} \cos \theta & - \sin \theta \\ \sin \theta & \cos \theta \end{array} \right ) \left(\begin{array}{c}x \\ y\end{array}\right)$

In your case you will probably be seeking a combination of rotation and translation between two images meaning that the relation between the points of these images can be defined as follows:

$\left(\begin{array}{c}x' \\ y' \\ 1\end{array}\right) = \left ( \begin{array}{ccc} \cos \theta & - \sin \theta & t_x\\ \sin \theta & \cos \theta & t_y \\ 0 & 0 & 1 \end{array} \right ) \left(\begin{array}{c}x \\ y \\1\end{array}\right)$

Or short:

$\textbf{x}' = \textbf{A}\textbf{x}$

This particular transformation has 3 degrees of freedom, one for rotation and two for translation. This means that you can solve this matrix for 2 point correspondences $\{ \textbf{x}' \leftrightarrow \textbf{x} \}$. That means in your case that you need two points in the original image and two matching points in the distorted image to ascertain the translation and rotation between them.

Ultimate reference: Multiple View Geometry

EDIT

The deformation between these images is obviously not an isometrie. For one the branches appear and disappear between them. I played a little with a slightly more involved form of feature based registration based on the assumption that the distortion is a similarity transform with 4 degrees of freedom, and the features will make up for the difference in edges. In short:

SIFT feature detection and matching
RANSAC estimation to detect inliers
Direct Linear Transform (DLT) variant for estimation of scale and rotation
Run of the mill cross correlation to detect the translation

Results, not brilliant, but here goes

All code on Github, demo file is so_13_02_2012.m.


Find similar Music using FFT Spectrums I have been experimenting some weeks to find a way to match/find similar Songs in a Library containing different Genres of Music.

My first try was to Detect Features like Tempo or how much Bass there is in Songs to form groups, but I didnt get far with this approach (Volume Changes Based Beat Detection) since on about 20% of the Songs Beat dont have to be counted always, sometimes 1/2 or 1/3 of them and I couldnt implement that.

After some weeks of failed trying I got a new idea which is described later in this Post. Simply put it works by taking Spectrum Samples of Files, making something like an "Average Spectrum" of Files to compare them.
The Idea behind was that for example Hardstyle has a lot more Bass than average Rock Music, I also verified this by looking at some Spectrums in Audacity.

File 1: Take full File FFT Spectrums (2048 Sample Size atm,Amplitudes Log. scaled)
Sum all Spectrum Arrays, take Averages of each Bin
Do the same some other Files, store all Results
Make List of FFT Values Differencies between File 1 and other Files
Make Average of Differencies between File 1 and File X
Sort Ascending by these Averages
Songs with the Lowest "Difference Value" are considered to be Similar.

Can some of you who have good knowledge tell me if this would be the right/good way to implement my Idea?

What you are trying to do has been tried over and over by hundreds of researchers and there is quite a large body of work about this. Check the proceedings of the ISMIR conference. Even if it is not up to date, read Elias Pampalk's thesis : http://www.ofai.at/~elias.pampalk/publications/pampalk06thesis.pdf

To quickly orient you on the right track:

Music can be similar according to many dimensions: a) timbre/texture/genre ; b) rhythmic pattern ; c) melody/chord progression... and more! From your message it's not clear what you want to measure!

If you are interested in a) the features you might want to look at are MFCC (Mel Frequency Cepstrum Coefficients), since they somehow capture the way human hearing works (frequency warping, log scale), since they are decorrelated (making modeling easier), and since they have lower dimensionality (13 coefficients vs 2048).
If you are interested in b), look at a feature called "Fluctuation Patterns" (Pampalk, in short autocorrelation of the signal in the 0.1 .. 10 Hz range over a few bands) ; or Whitman's "Penny" features (FFT of the MFCC along the time axis).
If you are interested in c), look at chromagrams. Start with Ellis' chromagram code (http://labrosa.ee.columbia.edu/matlab/chroma-ansyn/) then move up to Mauch's implementation if you need something more robust (http://isophonics.net/nnls-chroma).

That's for the features. Now you'll have to think of a better way to compare your songs once they have been represented as a sequence of those features. Computing pairwise differences between sequences is not very smart - eg: comparing a song and the same song offset by some silence will yield a difference while it is exactly the same! You'd rather compare the distribution of those features ; for example compute the mean / standard deviation of the features over song A and the mean / standard deviation of the features over song B and then take a probabilistic distance (KL, Bhattacharyya over those).

Last point, but which will matter later: computing the distance between a song and the rest of the corpus to find the nearest matches is quite inefficient. When dealing with large collections, techniques like LSH or Ball trees allow such nearest neighbors queries to be performed without explicit comparison with the whole corpus.

As an aside, tempo detection is an entirely different matter. If you want to look into it, best performance / accessibility paper on the topic is Ellis' Beat Tracking by Dynamic Programming. http://www.ee.columbia.edu/~dpwe/pubs/Ellis07-beattrack.pdf . It's incredibly simple but is close to state of the art algorithms.


Help with my first (simple) Z-transform I need to transform this Laplace function to the z-domain:

From the answer I received:

$s=(1-z^{−1})/T$

Then substitution into my Laplace function would give:

$t(z) = 2R/(m*(1-z^{−1})/T + 2R)$

Is that it?

Then $z$ is the input signal's current sample? Is $z^{-1}$ is the input signal delayed by one sample?

ie. Is it:

output_of_function = 2*R / (m* (1 - input_signal_1)/T + 2*R)

What basic mechanism might this function expect to perform? Eg. a LPF? Or what?

Does the fact that you're subtracting $z^{-1}$ from 1 imply that the input must be between -1 and 1?

Thanks. (I edited this in reply to the answer given from Matt L which helped.)

First of all, it's important to understand that there is no single best way to transform a continuous-time system to a discrete-time system. The method you're using is called backward Euler method, and it is defined by the mapping

$$s\leftarrow\frac{1-z^{-1}}{T}\tag{1}$$

Note that in $(1)$ you scale by $1/T$, where $T$ is the sampling interval (i.e., $1/T$ is the sampling frequency). Eq. $(1)$ just means that continuous-time differentiation (multiplication by $s$) is approximated by a first-order backward difference ($z^{-1}$ corresponds to a delay by $1$ sample). This and other methods for converting continuous-time systems to discrete-time systems are discussed in this answer.

Here is an example where I chose an arbitrary frequency range of interest $f\in[0,10]$, and an appropriate sampling frequency $1/T=20$. The figure below shows the frequency response of the continuous-time system (in blue), and the frequency response of the discrete-time system (in green), which was transformed from the continuous-time system according to $(1)$.


Why is the ROC of Laplace transform independent of imaginary part of s? An integral is defined as converging if it yields a finite value upon application of limits of integration. It is divergent otherwise.

Now sticking to the mathematical notation of Laplace transform, we have for a causal function $x(t) = u(t)$:
$$
X(s) = \int_0^\infty x(t)e^{-st}dt = \int_0^\infty e^{-st}dt
$$

$$
X(s) = -\frac{1}{s}e^{-st}\Biggr|_{0}^{\infty} = \frac{1}{s} = \frac{1}{(\sigma + jw)}
$$

By multiplying numerator and denominator with $(\sigma - jw)$, we get:
$$
X(s) = \frac{(\sigma - jw)}{(\sigma^2 +\omega^2)}
$$

For laplace transform not to exist, the denominator must become 0. Hence in this contrived example, both $\sigma$ and $w$ must be 0.

Conversely, if $\sigma = 0$ and $w \neq 0$, $X(s)$ exists as $\frac{-j}{\omega}$ with a magnitude of $\frac{1}{\omega}$.

Unless my basics are messed up, why in the literature is $jw$ disregarded as influencing $X(s)$? In other words, why is the ROC dependent only on $\sigma$?

I think your misunderstanding is that you manipulate an expression, which is only valid for $\text{Re}\{s\}>0$, in order to show for which values of $s$ it might be valid.

Note that the integral

$$\mathcal{L}\{u(t)\}=\int_0^{\infty}e^{-st}dt$$

only converges for $\text{Re}\{s\}>0$, so the result $1/s$ does not make any sense for other values of $s$.

Only the real part of $s$ in the term $e^{-st}$ inside the integral can provide the necessary damping such that the integrand decays sufficiently fast for the integral to converge. That's why the region of convergence only depends on the real part of $s$.


Merging pattern for FFT solutions I have the generic problem of analysing a discretly sampled signal, lets
say at 100Hz, and would like to compute features (ML application) of this signal at different scales:

An instantiation of this problem is computing the FFT on real-signal with
N samples:

Original signal: X[1],...,X[N]

Suppose I have already computed the FFT of the following
subarrays:

X[N-k+1],...,X[N]

X[N-2*k+1],...,X[N-k]

X[N-3*k+1],...,X[N-2*k]

X[N-4*k+1],...,X[N-3*k]

which is each of size k. Then I would like to compute the FFT
of the subarray

X[N-4*k+1],...,X[N]

(which is of size 4*k), WITHOUT reconstructing the original array X[.] and
running an ordinary FFT on the array of size 4k. Roughly I would like to do O(N) primitive operations in the best case if that is at all possible.

Is there some recursive property I could use?

PS: Note that if the goal of computing the FFT is substituted with e.g. finding the mean or the variance there exist easy recursive properties that allow merging of the results in constant time.

You can compute the DFT of the total signal from the DFTs of its parts but this is not efficient. What you would need to do is interpolate each DFT to the length of the total array (which is equivalent to zero-padding the sub-blocks to length $N$ before computing the individual FFTs), and then combine the interpolated DFTs with appropriate complex factors. This is computationally less efficient than going back to the time domain and computing the FFT of the whole sequence.

There is no recursion involved (unlike your examples of recursive mean and variance computation), simply because the data blocks are not related to each other. If you were to have a sliding window with a fixed length and you wanted to recompute the DFT after shifting the window by one sample, then you would indeed get a recursive formula for the DFT. This is known as the sliding DFT.


How do I set window length and shift length for computing STFT on Scipy? Does nperseg=32, noverlap=16 mean the window length is 32mS with a shift of 16mS?
Refer here for more details

for this function:

nperseg will set the number of samples to output from your STFT.
noverlap will specify how many old samples you want to take to calculate your future STFT.

Let's say you want a stft of size 256 with an overlap of 50% you will set noverlap = 128.
128 samples will then be taken from the old segment that was used to calculate the previous STFT...
For more information about STFT and the interest of overlap I invite you to watch this discussion.
STFT: why overlapping the window?


Reconstruction of Audio Signal from its Absolute Spectrogram I have the absolute Spectrogram of an audio signals.
I lost the phase data of the Spectogram because of various processing applied on the original spectrogram of the signal.

I'm trying to reconstruct the audio signal in a meaningful (Audibly) manner from teh absolute value only of the Spectrogram.
The obvious inverse won't work (The DFT inverse of the absolute, Since the Phase is significant).

The Spectrogram is a result of fusion of few audio signals as I'm trying to create a smooth transition between audio signals.

Anyone has experience with the problem?  Anyone has experience with this procedure?
Could anyone refer me to a code, article, etc...

Thanks.

One thing commonly done (for example in the source separation community) is to use the phase data of the original signal (before transformation where applied to it) - the result is much better than null or random phase, and not so far from algorithms aiming at reconstructing the phase information from scratch.

A classic reconstruction algorithm is Griﬃn&Lim's, described in the paper "Signal estimation from modiﬁed short-time Fourier transform". This is an iterative algorithm, each iteration requires a full STFT / inverse STFT, which makes it quite costly.

This problem is indeed an active area of research, a search for STFT + reconstruction + magnitude will yield plenty of papers aiming at improving on Griffin&Lim in terms of signal quality and/or computational efficiency.


Biquad Treble Shelf filter gain being weird Trying to compute Biquad Coefficients for use in a codec chip TLV320AIC3100.

Bass Shelf(boost and attenuation) work as expected. Treble Shelf attenuation works as expected but Treble Shelf boost is being weird in the sense that the amplitude difference between non treble frequencies and treble frequencies is as expected (Say 10dB) but instead of the overall gain being +10dB at treble frequencies, it ends up being +3dB(ish) and non treble frequencies end up being -7dB (instead of 0dB).

Used the TI BiQuad Calculator tool and also the tool here

Don't know where things are going wrong. Any guidance is much appreciated. Thank you.

Used the formulas from p3 of this PDF

// Coefficients from the Calculator
N0  2.336268076
N1  -2.773236549
N2  1.01094234
D1  -0.680740702
D2  0.254714569

N0_Real 2.336268076
N1_Real -1.386618274
N2_Real 1.01094234

Scaling Factor  2.336268076 (Max of the above 3 values)
Range   32767 (Fixed for the Chip - 16 bits)

Digital Coefficients
N0  32767
N1  -19448
N2  14179
D1  11153
D2  -8346

Some information:

Sampling Frequency: 44100
Treble Shelf Center Frequency: 7kHz (to get required boost/attenuation at 10kHz)

It's your scaling. Don't divide the coefficients by 2.336268076  (which is exactly the gain that's missing.

Instead use proper fixed point implementation techniques (scaling, shifting, "Q" notation, etc.) to deal with numbers larger than 1. Or much easier: use floating point, if you can.


Can someone explain why absolute is needed in power or energy signal? I know the power signal and energy signal's equation.

$$
P=\displaystyle\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{T/2}|x(t)|^2dt
$$

$$
E=\displaystyle\lim_{T\to\infty}\int_{-T/2}^{T/2}|x(t)|^2dt
$$

I am wondering whether $(x(t))^2$, instead of $|x(t)|^2$, is fine or not.

I think $(x(t))^2$ and $|x(t)|^2$ are same, aren't they?

Actually, I'm just guessing because of imaginary number, but I didn't learn or study that part, so I cannot make sure. Please let me understand. Thank you very much for reading my post.

Power and energy are real-valued quantities, so for a complex-valued signal $x(t)$ you need to use $|x(t)|^2$. For real-valued signals, you don't need the magnitude because $x^2(t)=|x(t)|^2$ holds. Note that this equality is not satisfied for general complex-valued signals.


How do I calculate peak amplitude of the signal components after zero padding and FFT? I am learning about DFT and trying to apply it to some audio processing. I am new to DSP but experienced in programming and have some background in math and physics. The FFT algorithm I use (lomontFFT) doesn't have 1/N scaling in front of the DFT sum (ie. it's defined just like on wikipedia DFT page). I will use simple example to describe my issue.

Let's say I have have 1024 samples of the pure 100Hz sine signal with peak amplitude of 1, and let the sample rate be 10240Hz. So that means that frequency resolution is 10240Hz/1024 = 10Hz. So when I perform FFT I get peak in 10th frequency bin which is correct. (since the signal is real I am only looking at the first half of the spectrum). Now, the magnitude of this peak in 10th bin is 512. How should I interpret this value? Is this magnitude somehow related to the energy of the 100Hz sine wave component in the signal?

Now, how can I get a spectrum to show amplitudes? I'd like to see what components are contained in the signal and what are their amplitudes. So in this example I should have one component at the 10th frequency bin with the value of 1 which is equal to the amplitude. I know I can divide the current value I got by 512 which seems to be the N/2. Is this correct way to calculate it?

If so, what about zero padding? If I zero-pad the same signal with additional 1024 zeros now I have total of 2048 samples. This means that frequency resolution is now 5Hz and my sine component is in the 20th bin. The magnitude is again 512, however now if I divide it by N/2=2048/2=1024 I get 0.5 instead of 1 for the peak amplitude. It seems like I should not count zero padded samples in this scaling? But what about if I just get the signal from somewhere and I don't know if signal is zero padded?

Also it seems that zero padding changes magnitudes in other bins. So in the first case without zero padding I get magnitude of 512 in the 10th bin and all other magnitudes are nearly zero (eg. values of the order 10^-10). But in the second case with zero padding, beside the magnitude value of 512 in the 20th bin I also get large magnitudes in bins around it (eg. 19th bin has magnitude of 334 and 21th bin has magnitude of 318 etc...). I guess these are side-lobes of spectral leaking but why are they showing up if I only have one sine frequency of 100Hz and frequency resolution is 5Hz so it falls exactly in the 20th bin?

I wish I could show some plots but I am programming this in C# console application so I can't provide visuals.

Note: I am using "frequency resolution" term here to denote "frequency bin spacing", I am aware that zero-padding doesn't increases frequency resolution since it doesn't add any new information to the signal.

It looks to me like your frequency resolution is sufficient and you are getting the correct value of the peak magnitude in your FFT but your scaling factor is incorrect.

You are right in dividing the magnitude of your FFT peak by N (1024 in your case). This number however, divides the power across the positive and negative sides of the frequency spectrum.

If you only want to look at one side of the FFT, you should multiply your peak magnitude by 2/N to get the correct amplitude that you are expecting. See this answer for another example.

e.g.

512 * 2 / 1024 = 1 (the answer that you were expecting)


Removing cyclic noise from signal I'm using a sensor with output that looks like the following figure.

The cyclic noise is common to ask measurements made by the sensor. Is there a way to remove the noise without negatively effecting higher frequencies.

We have attempted using a Fourier transform and then zeroing the frequencies around 6 Hz in Matlab, but that was ineffective.

You must use notch-stop or band-stop filters.

There are different methods to implement them, But you could implement them easily using fdatool in MATLAB. After you designed the filter, filter your signal using filter function.


Why does it mean that the process/signal is not stationary when its variance varied with time? Why does it mean that the process/signal is not stationary when its variance varied with time? that is,

$VAR[X(t)]= \alpha \times t$,$t$ is time,and $\alpha$ is a constant,then $X(t)$ is not the WSS process

From the Wiki page - A stationary random process is a stochastic process whose unconditional joint probability distribution does not change when shifted in time. Consequently, parameters such as mean and variance, if they are present, also do not change over time.

Additional References are: Link1 and Link2.

The implication from the definition is that the mean and variance of random process do not vary over time. Since by your definition of the variance - it varies over time, then your process is non-stationary.

A couple of notes:

Just because the mean and variance don't vary over time does not imply that the process is stationary. The definition of a stationary random process has more implications.
A random process where the mean and auto-covariance are invariant to time shifts is often called a "Wide Sense Stationary" process. There are no conditions on the higher order moments (skewness, kurtosis, etc). This is a less stringent version of stationarity that is often more applicable for real-world signals.


Short lowpass filter that also extrapolates lines I'm looking for a short lowpass filter that also extrapolates lines,
[0 1 2 3] --> 4,
so I need (correct me) $H([\ 0\ \ 0.25\ \ 0.5\ ]) \approx [ 1\ \ 2\ \ 0 ]$ .
How close can one get to this with 2 to 4 FIR and 1 IIR coefficient ?

Added: for a 3-term FIR, I want

Input     -> output
 1  1  1  -> 1  # const -> const
-1  0  1  -> 2  # extrapolate lines
 1 -1  1  -> 0  # low-pass

Solving these 3 equations gives the FIR [1.25 .5 -.75] .
I can do longer FIRs by least squares + fiddling. But how do I do FIR x 1-term IIR ?
Is there a general approach to filters that extrapolate lines ?

The most straight-forward solution is a length $N$ FIR filter with impulse response

$$h[n]=\begin{cases}\frac{N}{N-1},&\quad n=0\\
-\frac{1}{N-1},&\quad n=N-1\\0,&\quad\text{otherwise}\end{cases}$$

For $N=3$, the output sample computed from three input samples x = [1 3 5] is $\frac32\cdot 5 - \frac12\cdot 1 =7$, as it should be. Note that this filter is very cheap, it has only two non-zero taps regardless of the value $N$.


Matrix form of 2D-DFT for a vectorized image I want to apply 2d DFT to a N by N image. However, image is vectorized such that it is NxN by 1. How can i find the matrix form of 2d DFT such that resulting vector from multiplication is the vectorized NxN by 1 2d DFT of my N by N image?

In matrix form, we can write the 2-D DFT of an Image $\mathbf X$ via $\hat{\mathbf{X}} = \mathbf F \mathbf X \mathbf F^{\rm T}$.

This expression can be vectorized with the help of the Kronecker product. Using [*], we obtain $$\hat{{\mathbf{x}}} = {\rm vec}\{\hat{\mathbf{X}}\} = (\mathbf F \otimes \mathbf F) \cdot \mathbf{x}.$$

This shows that the matrix you need is $\mathbf F \otimes \mathbf F$.


Channeling two named PIPEs with a GNU Radio flowgraph | String transmitted on one end does not reach the other end 1)
I created two named PIPEs:

$ mkfifo pipe_in
$ mkfifo pipe_out

2)
Then connected pipe_in to pipe_out with the following GNU radio flowgraph:

File Source config:

Repeat: No
Vec Length: 1 (default setting)
Add begin tag: pmt.PMT_NIL (default setting)

File Sink config:

Unbuffered: Off
Vec Length: 1 (default setting)
Append file: Overwrite

3)
Finaly, I start running the flowgraph, then write a string to pipe_in hoping to see it reach pipe_out by running the following two commands in separate terminals:

TERM2 $ tail -f pipe_out

and

TERM1 $ echo 'hello world!' > pipe_in

Unfortunately the string never reaches pipe_out. Can someone shed some light please?

That's how tail works, it waits until the end to print the last n lines. If you use cat instead, it should work.

I hope it helps.


Phase-shift filter design If a time-domain signal has sharp corners, its frequency spectrum will contain high-frequency components. Truncating the spectrum results in Gibbs' phenomenon. So if you're trying to design an FIR, you really want the target frequency response to be nice and smooth so that windowing the impulse response down to a finite length doesn't distort the frequency response too much.

Currently I'm contemplating trying to design a very strange filter: One that has unit gain at all frequencies, but non-zero phase. I'm wondering whether or not a similar phenomenon occurs: If the filter has unit gain at all frequencies, then what does truncating the impulse response do to the phase alignment?

This would be an allpass filter. Except for the trivial case of unity and integer-sample delays, these can't be done as FIR filters and in general an IIR filter is required. However, they are easy to make. The zeroes of an allpass are simply the inverse of the poles (and vice versa). If you have the poles in polynomial form, you can simply flip them to get the zero polynomials. For example a second order allpass looks like this
$$H(z)=\frac{a_{2}\cdot z^{0} + a_{1}\cdot z^{-1} + a_{0}\cdot z^{-2}}{a_{0}\cdot z^{0} + a_{1}\cdot z^{-1} + a_{2}\cdot z^{-2}}$$
Strict allpass filter have $\left \|H(e^{j\cdot \omega })\right \|=1 $ for all frequencies. You can certainly design approximation using FIR filters if you only need this property for a limited frequency range and if the magnitude doesn't have to be exactly unity.


How to select the specific frame with object I am detecting the object from the live camera through feature detection with svm , and it read every frame from camera while predicting which affect its speed , i just want that it should select the frame which contain the object and ignore other frames which have no object like empty street or standing car's , it should only detect the moving object

For example , If the object came into camera in 6th frame , it contain into the camera till many frames until it goes out from camera's range , so it should not recount the same object and ignore that frames.

Explanation :

I am detecting the vehicle from video , i want to ignore the empty frames , but how to ignore them ? i only want to check the frames which contain object like vehicle , but if the vehicle is passing from video it take approximately lets assume 5 sec , than it mean same object take 10 frames , so the program count it as 10 vehicles , one from each frame , i want to count it as 1 , because its the one (SAME) vehicle which use 10 frames

My video is already in Background subtraction form , i used SURF/BOW with SVM

I explore two techniques :

1- Entropy ( Frame subtraction )
2- Keyframe extraction

I'm seeing two issues here:

Frames with no activity: You can check the gradient between successive frames and only do object detection when the gradient information passes a certain threshold.
Counting the same object: The vehicle location information can be useful for this task. If the gradient change occurs in the location of a detected object in the previous frame then it is the same vehicle. If necessary, you can also use object color and object size parameters to account for varying speeds.


MFCC, uncertain of my results and algorithm I've started to implement an MFCC algorithm on an embedded DSP, using the CMSIS DSP library. However the results I'm getting don't appear consistent. This is my algorithm so far:

Gather 32ms PCM sample from microphone, 512 samples long. 16000hz
Apply a Hamming window.
Compute the 512 FFT using this.
Compute the magnitudes to get power spectrum.
Generate a set of 30 filter banks - using 32 points on the mel spectrum
Convert my power spectrum magnitudes from 16bit ints to floats.
Get the dot product for each filter with the power spectrum.
Log each dot product.
This is the uncertain part. The DSP library I'm using only has DCT IV, and this only supports sizes 128 and up. My set of logs is 32. So I've opted to use a Real FFT function, the results I get from this vary quite a lot. If I speak the same word in the same tone.

Could this be what is hampering my result? Real FFT provides an output of both real, complex and conjugate data (described in link above) - which would I use for the MFCC? - Should I even expect my MFCCs to be similar or is it expected to vary a lot? My results tend to vary between -20.0 and 20.0.

Thanks in advance.

For the computation of MFCC's you usually use a DCT-II. Just implement it yourself if you can't find a proper implementation anywhere. You only have relatively few log spectral values, and probably even fewer DCT coefficients that you will use, so I hope that computational complexity will be no big issue. This document is a nice overview of all the steps necessary for obtaining MFCC's.


Python: Least Squares Support Vector Machine (LS-SVM) I'm looking for a Python package for a LS-SVM or a way to tune a normal SVM from scikit-learn to a Least-Squares Support Vector Machine for a classification problem.

The goal of a SVM is to maximize the margin while softly penalizing points that lie on the wrong side of the margin boundary. The function that is used is a Quadratic Programming (OP) problem.

A LS-SVM which defines a least squares cost function and replaces the inequality contraints with equality constraints and is a Linear Programming (LP) problem.

Wikipedia: LS-VSM

Does someone ever tried something like that and can help me?

There is a package called FukuML.

In their description (Version 0.4.1) they write:

Support Vector Machine

Primal Hard Margin Support Vector Machine Binary Classification Learning Algorithm
Dual Hard Margin Support Vector Machine Binary Classification Learning Algorithm
Polynomial Kernel Support Vector Machine Binary Classification Learning Algorithm
Gaussian Kernel Support Vector Machine Binary Classification Learning Algorithm
Soft Polynomial Kernel Support Vector Machine Binary Classification Learning Algorithm
Soft Gaussian Kernel Support Vector Machine Binary Classification Learning Algorithm
Polynomial Kernel Support Vector Machine Multi Classification Learning Algorithm
Gaussian Kernel Support Vector Machine Multi Classification Learning Algorithm
Soft Polynomial Kernel Support Vector Machine Multi Classification Learning Algorithm
Soft Gaussian Kernel Support Vector Machine Multi Classification Learning Algorithm
Probabilistic Support Vector Machine Learning Algorithm
Least Squares Support Vector Machine Binary Classification Learning Algorithm
Least Squares Support Vector Machine Multi Classification Learning Algorithm
Support Vector Regression Learning Algorithm

As you can see, they have the feature you're after.


Understanding of Random Process/Random Variable At a simpler level to my previous question, I wanted to confirm my understanding on Random Process based on Random Variables using an example.

So, I took this example:
If we consider a dice, which has six sides denoted by fi (i can be 1,2,3,4,5,6). So, assigning outcomes f1=10,f2=20,f3=30,f4=40,f5=50,f6=60 thereby constructs a Random Variable X.

Now, if we assign the outcomes: f1=15,f2=25,f3=35,f4=47,f5=80,f6=90 at a different time instant t2 and at time instant t3, we assign the outcomes: f1=43,f2=32,f3=16,f4=89,f5=99,f6=56.

So, if I plot this variation on a Graph, I am getting something like this below:

Now does this set of sample functions (Ensemble) constitute a Random Process?

Let me explain it in another way. Consider you have 6 different function of time.
You only throw your dice once and regarding the outcome you choose on of six functions and the chose one is one outcome of your random process.

Now expand this concept to all possible functions of time instead of only 6 functions and instead of dice throwing consider another form of random event which could append a probability to all (infinite) possible functions. The set of all these possible functions is your ensemble.


Is there a closed form for the average energy of a QAM constellation? I'm trying to come up with a formula for the average energy of the signals in a M-QAM constellation.
If $b$ is the number of bits represented by each point in the constellation, such that $M=2^b$. For $b$ even, I know the formula is:

$$ \bar{E} =  \frac{(M-1)d^2}{6} $$
where $d$ is the distance between two points in the constellation.

This can be easily verified using the average of all QAM constellation points generated by MATLAB, whose default value for $d$ is $2$:

b=4;
M=2^b;
sum(abs(qammod([0:M-1],M,0,'gray')).^2)/M

Is there a formula for the average energy in a QAM constellation with $b$ odd?

From $\bar{E}=\frac{(M-1)d^2}{6}$ as the average energy in $M$-QAM, $M$ a perfect square and a square constellation, we have that the total energy is
$$E_{\operatorname{tot}} = M\frac{(M-1)d^2}{6}$$
Now, since the $M$-QAM signal is the sum of two orthogonal $\sqrt{M}$-PAM
signals of equal energy, then, writing $N$ instead of $\sqrt{M}$,
we have each of the orthogonal signals carrying a total energy of
$$\displaystyle \frac{1}{2}\frac{M(M-1)d^2}{6}
= \frac{N^2(N^2-1)d^2}{12}.$$
Since there are $M = N^2$ constellation points, the average energy of
a $N$-PAM signal with spacing $d$ is
$\displaystyle \frac{(N^2-1)d^2}{12}$.

Now, if $b$ is odd, a $2^b$-QAM signal constellation is
a rectangular grid with $2^{(b-1)/2}$ points in one direction and
$2^{(b+1)/2}$ points in the other.  Thus, the $2^b$ points contribute an
average of
$$\frac{(2^{b-1}-1)d^2}{12} +  \frac{(2^{b+1}-1)d^2}{12}
= \frac{(2.5M-2)d^2}{12}$$

Note that it is possible to use constellations other than rectangular when
$b$ is odd. For example, the $32$-CROSS constellation is a (rotated) version
of a $36$-AM square constellation with the corner points deleted. Rotation
does not change the energy and so we can get the answer for $32$-CROSS by
finding the total energy for a $36$-QAM square constellation, subtracting
off the energy of the $4$ corner points, and dividing what is left by $32$.


Seemingly no division operation is done of the diagram of the transfer function Given the below system.
$$  y\left[n   \right] =2 x\left[n \right] -x\left[n -1\right] +0.5 y\left[n -2  \right]  $$
We'll find out the transfer function of it.
So, easiely using z-transform,
$$  Y(z) =2 X(z) -X(z) z^{-1} +0.5 Y(z) z^{-2}  $$
$$ Y(z) \left\{ 1-0.5 z^{-2}  \right\} =X(z) \left\{ 2-z^{-1}  \right\}   $$
$$ H\left( z  \right) = \frac{Y(z) }{X(z) }= \frac{2-z^{-1} }{1-0.5 z^{-2} }   $$
gained it.
Currently I can't get the below diagram.
I can't see component which operates division so it is weird that the diagram is applied to the transfer function.
About the blue domain of below, currently I can't get how the 1 of the denominator of the transfer function is attained.

What I've been missing or what should I study for next?

The diagram implements the time domain difference equation and NOT the frequency domain transfer function.
The difference equation has no divisions and therefore there are no divisions in the block diagram. The Z-transform is useful to analyze the frequency domain behavior of a system but you don't implement it directly. (at least not at this point in your learning).


How do we actually obtain samples to form discrete time samples? I'm in my second week of an online DSP class and I'm still kinda lost on some of the concepts of what we are actually doing.  I'm hoping you guys can help clear up my understanding.  Please forgive me if what I'm asking makes no sense.

I'm reading about continuous time sampling, which I don't think I fully understand yet, and discrete time sampling.  But I guess to understand those concepts I really need to understand sampling.

I think sampling is gathering discrete data points from a continuous time signal?  But how?  How is this physically done?  What are we using to get the actual data to convert it into a discrete time sample?  Is there a good example I can look at to understand this more?

Right now the text book just shows math formulas, pictures of graphs of sines with various frequencies and amplitudes, and with dots representing discrete time samples vs continuous time samples - but no application of what is actually happening?  For example, is an instrument making a consistent sound, that is being captured through some microphone and producing the signal, which is then passed through some system that calculates data points?

Is a continuous time signal just a signal being actively captured and a discrete signal just data points over a certain period of time?  How and when do I apply the mathematics we are learning about (as we haven't needed to really do so yet)?

Assume you are speaking, those vibrations (sound) are converted into electrical signals(continuous time signal) by Mike. After that, a Sample & hold circuit is placed, which stores the converted electrical signal inside a capacitor in the form of charge (now u get discrete time signal).

Every signal is  sampled at some frequency called sampling frequency. Let the sampling frequency in this case be 8KHz. It means, at time t=0, you sample (read) the electrical signal amplitude (by using sample & hold circuit) and hold it(store it) in a capacitor. You wait till some time and when time t=0.00125 seconds(1/8KHz, remember our sampling frequency is 8KHz), you again sample the electrical signal and store it . You repeat this process every 0.00125 seconds. Now you have got discrete version of your continuous time electrical signal (in fact speech signal).

In the time between every samples(0.00125 seconds), the charge stored in the capacitor is converted to digital binary values using special circuits called ADCs. These converted digital values and discrete signals in binary form(1s and 0s which may be stored or for further processing in a computer).

If you are interested, read about how sample and hold circuit works. You need some electronics background for that.


How to perform 2D image rotation by FFT based image translation using matlab? I am studying rotation using 3 translations as described in this scientific paper. According to the authors , it is possible to decompose a rotation into 3 shear translations which will be implemented in the frequency domain by FFT according to time-frequency shifting theorems.

For example the first translation is obtained as in equation (13) by doing the following operation for every pixel considering that the horizental coordinate x is variable and the vertical coordinate y is fixed:

Gx(x,y)=FFT-1(exp(-2*piau*y)*FFT(f(x,y)) where f is the value of the pixel x,y , a is a constant for the translation and FFT-1 is the inverse fourrier transform and u is x mapped into the frequency domain.

To code this transformation, one should think as follow:
We have to compute  Gx by changing the range of U from -M/(2*M) to (M/2*M) and by iterating y from 1 to N.

Surprisingly, the suggested implementation doesn't follow the procedure above:
Can somebody develop further why the implementation is as follow.

close all;
clear all;
I = double(imread('Put image here'))/255;
%imshow(I)
I=rgb2gray(I);
facteur=4;
[M, N, s] = size (I);

temp=ones(2*M,2*N);
temp(floor(M/2):floor(M/2)+M-1,floor(N/2):floor(N/2)+N-1)=I(:,:);
I=temp; clear temp;
imshow(I)
[M N s]=size(I)
 teta=pi/8;
a=tan(teta/2)
b=-sin(teta)

 Nx = ifftshift(-fix(M/2):ceil(M/2)-1);
 Ny = ifftshift(-fix(N/2):ceil(N/2)-1);
 Ix=zeros(length(Nx),length(Ny));
 for k=1:M
   Ix(k,:)=(ifft(fft(I(k,:)).*exp(-2*i*pi*(k-floor(M/2))*Ny*a/N)));
 end
 %imshow(Ix);

The following is the completed code that yields the expected result:

close all; clear all; clc;

I = double(imread('Cameraman.tif'))/255;
figure,imshow(I)

[M, N, s] = size (I);

temp=ones(2*M,2*N);
temp(floor(M/2):floor(M/2)+M-1,floor(N/2):floor(N/2)+N-1)=I;
I=temp; clear temp;

figure,imshow(I)
[M N s]=size(I)
teta = pi/3;
a = tan(teta/2)
b = -sin(teta)

 Nx = ifftshift(-fix(M/2):ceil(M/2)-1);
 Ny = ifftshift(-fix(N/2):ceil(N/2)-1);
 Ix=zeros(length(Nx),length(Ny));
 for k=1:M
   Ix(k,:)=real(ifft(fft(I(k,:)).*exp(-2*i*pi*(k-floor(M/2))*Ny*a/N)));
 end

 Iy = zeros(size(Ix));
 for k=1:M
   Iy(:,k)=real(ifft(fft(Ix(:,k)).*exp(2*i*pi*(k-floor(M/2))*Ny*b/N)'));
 end

 If = zeros(size(Iy));
 for k=1:M
   If(k,:)=real(ifft(fft(Iy(k,:)).*exp(-2*i*pi*(k-floor(M/2))*Ny*a/N)));
 end

 figure,imshow(Ix);
 figure,imshow(Iy);
 figure,imshow(If);

So this works as expected from a 2D rotation about z-axis. The paper argues that this is a fast and accurate method, but I doubt that, as a simple 2D rotation could easily be implemented using raw trigonometric computations and a bilinear interpolation (or bicubic if a sharper result is desired) other than forward and inverse FFTs.

I can say that the method is slow (on matlab, it can be further vectorized however). On the other hand the perceptual quality of the resulting rotated image is fairly good except some little ringing artefacts towards the edges.

I suggest you the ch-5 of Principles of Computer Graphics 2E by FOLEY for further understanding of the mathematics used in these algorithms and applications.


What is the best approach to classify an image containing either text or image? I am doing a research about image-based User Interface (UI) Analysis, and while UI image-based analysis is not more complicated than normal image analysis, I was able to do a segmentation process to separate text, icon and picture in the UI.

However, my next step is to perform content classification on the segmented groups.  The segmentation contains:

text
image or icon

So: what could be the best approach to classify an image containing either text or image?

The segment with text (I assume with a background) should contains text clearly visible from the background, then the text segment should have a lot of peaks in the transitions from text to background. The segment icon normally shouldn't have many peaks because in general there are few variations in an image. So if you make a stft of a part of a segment the text should have more harmonics than the icon.
However if the icon is similar to the text or if it has many peaks this idea can't works. You can also find the color density, a single color text on a single color background has the color density concentrated in two points. This is only an idea and it must to be improved.


Creating a matrix of Gaussian wavelets at dyadic scales I need to create a diagonal matrix containing the Fourier coefficients of the Gaussian wavelet function, but I'm unsure of what to do.

Currently, I'm using function ConstructHaarWaveletTransformationMatrix to generate the Haar wavelet matrix and am taking the rows at dyadic scales (2, 4, 8, and 16) as the transform:

M = 256
H = ConstructHaarWaveletTransformationMatrix(M);
fi = conj(dftmtx(M))/M;
H = fi*H;
H = H(4,:);
H = diag(H);

etc.

How do I repeat this for Gaussian wavelets? Is there a built-in MATLAB function which will do this for me?

For reference, I'm implementing the algorithm in section 4 of the paper in Compressed Sensing for Wideband Cognitive Radios.

I need to do edge detection on a 1D signal, which is constructed like so:

edges = [50, 120, 170, 200, 220, 224, 256];
levels = [24,3,30,0,36,0,0];
idxs = zeros(1,max(edges));
idxs(edges(1:end-1)+1) = 1;
psd = levels(cumsum(idxs)+1);

When I just wavdec and the differentiation matrix:

[y1,y2] = wavedec(psd, 0, 'haar');
z0 = gamma*y1.';

I get exactly what I want.

The point, however, is that I want to recreate this edge information compressively (that is, by solving a linear program).

The fact that using the matrix above creates gibberish is my current sticking point. I want to replace it with something else and see if it still works.

In this case the problem was that I was working in the frequency domain, but implementing an algorithm in the time domain.

I've removed the pre-multiplications of the Fourier matrix in the time domain and everything works fine.


Multiple biquads vs higher order filtering I've heard that you can stack biquads on top of each other to get the same (or similar) effect as higher order filters.

For instance a biquads is order 2, so putting something through 3 biquads is like an order 6 filter.

I've heard that doing this is more numerically stable than an actual order 6 filter.

I was wondering though, besides increased numerical robustness, are there any differences between doing this compared to just using a higher order filter?

Does it impact processing speed or frequency response or anything else?

For some context, this for audio synthesis for electronic music, on a modern desktop computer, using 32 bit floating point math.

A cascade of 2nd order biquadratic (biquad) sections is a commonly used approach for  implementing higher order FIR/IIR filters in fixed point systems (Once you determine that you need an IIR filter, this would be a topology to consider). A big consideration for the biquad structure is stability: higher order IIR filters are prone to instability in fixed point systems due to coefficient rounding combined with long feedback delay paths (causing poles to go outside of the unit circle). The absence of any feedback between the individual cascaded sections (decoupling the poles and zeros between sections) results in a more stable structure that will be less sensitive to coefficient quantization.
I do not know of other advantages, but the numerical robustness is a very serious consideration especially  in higher order structures (and fixed point systems).


Fourier Transform with both Time Delay and Frequency Shift I know that the Fourier transform of a function with time delay can be written as: $$\mathscr{F}\big\{x(t-t_0)\big\}=X(f)e^{-j2\pi f t_0}$$

The Fourier transform of a function with frequency shift can also be written as: $$\mathscr{F}\Big\{x(t)e^{j2\pi f_0 t}\Big\}=X(f-f_0)$$

So what if we have both shift and delay at the time domain, what will be the result in the frequency domain? E.g.:

$$\mathscr{F}\Big\{x(t-t_0)e^{j2\pi f_0 (t-t_0)}\Big\}$$

Will the result be: $$X(f-f_0)e^{-j 2 \pi  f (t-t_0)}$$

Also what will be the result of:

$$\mathscr{F}\Big\{x(t-t_0)e^{j2\pi f_0 t)}\Big\}$$

Is there an order to apply these properties?

If you are ever unsure, just go back to the definition and work out the Fourier Transform property for the specific situation:

$$\begin{align*}\mathscr{F}\left\{x\left(t-t_0\right)e^{j2\pi f_0\left(t-t_0\right)}\right\} &= \int_{-\infty}^\infty x\left(t-t_0\right)e^{j2\pi f_0\left(t-t_0\right)} e^{-j2\pi f t}dt\\
\\
&= \int_{-\infty}^\infty x\left(\tau\right)e^{j2\pi f_0\tau} e^{-j2\pi f \left(\tau+t_0\right)}d\tau \\
\\
&= e^{-j2\pi ft_0}\int_{-\infty}^\infty x\left(\tau\right)e^{j2\pi f_0\tau} e^{-j2\pi f \tau}d\tau \\
\\
&= e^{-j2\pi ft_0}\int_{-\infty}^\infty x\left(\tau\right) e^{-j2\pi (f-f_0) \tau}d\tau \\
\\
&= e^{-j2\pi ft_0} X(f-f_0)\\
\\
\end{align*}$$

$$\begin{align*}\mathscr{F}\left\{x\left(t-t_0\right)e^{j2\pi f_0 t}\right\} &= \int_{-\infty}^\infty x\left(t-t_0\right)e^{j2\pi f_0t} e^{-j2\pi f t}dt\\
\\
&= \int_{-\infty}^\infty x\left(\tau\right)e^{j2\pi f_0(\tau+t_0)} e^{-j2\pi f \left(\tau+t_0\right)}d\tau \\
\\
&= e^{-j2\pi (f-f_0)t_0}\int_{-\infty}^\infty x\left(\tau\right)e^{j2\pi f_0\tau} e^{-j2\pi f \tau}d\tau \\
\\
&= e^{-j2\pi (f-f_0)t_0}\int_{-\infty}^\infty x\left(\tau\right) e^{-j2\pi (f-f_0) \tau}d\tau \\
\\
&= e^{-j2\pi (f-f_0)t_0} X(f-f_0)\\
\\
\end{align*}$$


Frequency resolution I have made a stft_func in MATLAB which takes window length i.e t_short  as a parameter for window length. Now, this function is giving me a very blurred image with $f_s=1000\mathtt{Hz}$, radar target $10\mathtt{GHz}$ (fm) and $T_{short}=0.5 \mathtt{s}$
How to find frequency resolution from spectrogram.?

Frequency resolution of your spectrogram is connected to time window length and sampling frequency by following relation:

$$\Delta f = \dfrac{f_s}{N}$$

where $f_s$ is the sampling frequency, and $N$ is length of your time window in samples. $N$ can be easily expressed in seconds, in your case length of window is $T=0.5 \mathtt{s} $ and $f_s = 1000 \mathtt{Hz} $. Therefore window length in samples is: $$N=T\cdot f_s = 0.5 \cdot 1000=500$$
Which is giving you frequency resolution to be:

$$\Delta f = \dfrac{1000}{500}= 2 \mathtt{Hz}$$

You might also notice that substitution of $N$ with respect to $T$ and $f_s$ yields:

$$\Delta f = \dfrac{f_s}{T\cdot f_s}= \dfrac{1}{T} $$

where $T$ is obviously length of your window.


voice pattern matching algorithm My idea is to find a word in wav/mp3 files. It is just one word.
I have a database with a few hundred songs.

I am looking for pattern matching (this is not an exact pattern) it can be an word i choose but only one.

Can you recommend a pre-processing algorithm. you run the algorithm on the files?
Can you recommend a real-time algorithm?

I have found this online

If you are really interested in voice recognition. Check out Hidden markovs model. It is used by several softwares to recognize your words


Computer Vision / Speech Recognition I am not sure if the question is off topic but I need some input before I can finalize anything.

Actually to do my higher studies (in research) I wanted to make sure that I land my career in Digital Image / Signal Processing only. I have a strong background in Computer Application development. Now the Signal Processing course comes under Electronics And Communication Engineering / domain while Computer Vision & Speech Recognition subjects comes under Computer Science.

I am not sure if these subjects are different or doing Computer Vision / Speech Recognition is nothing but DSP / DIP. Is MATLAB being used for Computer Vision / Speech Recognition software?

Should I opt for Electronics And Communication Engineering course or one from Computer Science? Kindly guide me appropriately

These fall under signal processing and mostly you will find electrical engineers who are working in domain of speech recognition. The thing with computer vision is that it overlaps between both field.But as you mentioned you only have option of Electronics/Communications, Electronics is not much related to these field ,well in communications, we all use statistical signal processing which is also useful in speech processing and computer vision.

If you can get into Electrical Engineering(signal processing) that would be much better but out of these two I would suggest go for Computer Science if you plan a future in these field.


A DSP simple book/reference? I have some basic maths skill up to derivatives.
I ask you if you know any simple/concise book that can introduce me to DSP world. Thank you.

I can recommend the following two books (which are free!):

Introduction to Signal Processing by S.J. Orfanidis
Signal Processing for Communications by P. Prandoni and M. Vetterli

Also take a look at the answers to this question.


Extracting the DC component of a signal, that proportional to the displacement about the y-axis I’m not sure that I clearly enough reflected the essence of a small problem in the title of the topic, but I’ll try to reveal its essence in the question itself.

https://www.wolframalpha.com/input/...

We have two signals.

If the signal is shifted more to the positive region (up along the $y-axis$), then we "mentally associate the $+1$ sign with it, if it is negative (down the $y-axis$), then we" mentally associate the $-1$ sign with it".

If both the positive and negative regions are approximately the same, then the DC component should be approximately equal to zero.

How is it possible to extract the DC component of the signal proportional to the displacement about the y-axis:

1. Use mathematical functions? If so, which ones?

2. Use filters? If so, which filter to choose to solve this problem?

I am an amateur in signal processing, I began to study this theme not so long ago.

DC component is the zero frequency component. If you look at the DFT expression, this would mean summing up the input. So in mathematical terms this would be the sum (or the running sum) based on your application.

In terms of filtering again this would be an averaging filter. Like moving average.


what is the possible solution for the Intra-class variations poblems? my poblem is that i have a database for human dorsal hand veins ,but with intra-class variations , so this problem have a large effect on my results when i need to do matching between my samples for recognition purpose(testing phase after finishing trainning phase ) .

therefore , i need to know what are the possible solutions that may be suitable here.
are there any function in matlab can overcome this problem?

examples of my images(5 images for the same person):

class 1:

class 2:

class 3:

and the code that i wrote to obtain the above images was as below:

    in_dir = 'E:\master_matlab\HandVein_DataSet';
    out_dir = 'E:\master_matlab\final resultant images';
    for px = 1 : 50,
        for hx = 1 : 5,
            fname = sprintf('%04dhv%d.bmp', px, hx);
            current_image = imread(fullfile(in_dir, fname));

            % do processing

    level=graythresh(current_image);
    x=im2bw(current_image,level);
    % figure,imshow(x);
    nnew=current_image.*cast(x,class(current_image));
    medimage=medfilt2(nnew,[5 5]);%median Filtering image
    % figure,imshow(medimage);
    z=adapthisteq(medimage);
    % figure,imshow(z);
    H = fspecial('unsharp');
    y=imfilter(z,H);
    % figure,imshow(y);
    mIm=imfilter(y,fspecial('average',15),'symmetric');%%%%%%here mask is 15*15 and method is symmetric(not replicate)
    sIm=y-mIm;
    bw=im2bw(sIm,0);                          % Convert to binary image
    ALT_img=imcomplement(bw);                 % Complement binary image
    % figure,imshow(ALT_img);
%%%%%% morphological on binary images.

open_img =bwmorph(ALT_img,'open'  ,Inf);
% figure,imshow(open_img);title('open_img');
major_img =bwmorph(open_img,'majority' ,Inf);
ske_img = bwmorph(major_img,'skel',100);

%removal outer boundary
without_border=ske_img.*x;
% figure,imshow(without_border);

   %%%%spurs removal can be added  or removed
spur_img=bwmorph(without_border,'spur',1);%8 spurs
% figure,imshow(spur_img);

% centering the resultant spur_image
measurements = regionprops(im2double(spur_img), 'Centroid');
    [rows columns] = size(spur_img);
    rowsToShift = round(rows/2- measurements.Centroid(2));
    columnsToShift = round(columns/2 - measurements.Centroid(1));
    shiftedImage = circshift(spur_img, [rowsToShift columnsToShift]);

    %%%%%cropping to remove the below white lines

rect=[26 ,37 ,255 ,138];
ssdc=imcrop(shiftedImage, rect);
q12=imresize((ssdc),1.25);
% figure,imshow(q12);
%
        % save output
        imwrite(q12, fullfile(out_dir, fname));
    end
end

regards for all

It boils down to the features.  You have shown some images of the same class. You also need to look at images of the other class(s), and find a set of features that maximize inter-class variance and minimize intra-class variance. What features those might be depends on your data. You need to look at the types of variation that exists within each class, and try to find features, transformations, or distance measures that account for it. It is very hard to tell you anything more specific based on the information you have provided.

Edit:
One thing you can do is post some images of a different class.

Also, are the images aligned? Maybe before you compute the features you should translate all your images so that the white pixels coincide as much as possible.  Then you can define a distance measure between two images as the percentage of white pixels that co-locate, and then use a k-nearest-neighbor classifier.

Some other random thoughts:
Shape context features may be suitable for your problem.

You also might want to search the literature on recognizing hand-written characters.


Fourier transform and anti-trasform--identity missing I have a very silly doubt:

If we define the power spectral density:

S(f)=$\frac{1}{2\pi}\int exp(-i\tau2\pi f)r(\tau)d\tau$ (1)

where $r(\tau)$ is the correlation coefficient.

If we do the Fourier anti-transform, we obtain
$r(\tau)=\int exp(i\tau2\pi f)S(f)df$ (2)

Now my doubt is: if I substitute in the second equation the first equation, it seems to me that I don't find the identity
$r(\tau)=r(\tau)$

I hope you can help me, maybe I am missing something

Different ways of showing it, depends where you start. Are you willing to accept that the Fourier transform of $\delta(t)$ is $1$ and vice versa (i.e., $\int_{-\infty}^\infty {\color{red}1} \cdot {\rm e}^{\jmath 2\pi f t} = \delta(t)$)? If so, it's easy:

$$\begin{align}
r(\tau) & = \int_{-\infty}^\infty {\rm e}^{\jmath 2\pi f \tau} S(f) {\rm d}f \\
& = \int_{-\infty}^\infty {\rm e}^{\jmath 2\pi f \tau} \int_{-\infty}^\infty {\rm e}^{-\jmath 2\pi f t} r(t) {\rm d}t {\rm d}f \\
& = \int_{-\infty}^\infty \int_{-\infty}^\infty {\rm e}^{\jmath 2\pi f (\tau-t)} r(t) {\rm d}t {\rm d}f \\
& = \int_{-\infty}^\infty r(t) \int_{-\infty}^\infty {\color{red}1} \cdot {\rm e}^{\jmath 2\pi f (\tau-t)} {\rm d}f {\rm d} t \\
& = \int_{-\infty}^\infty r(t) \delta(\tau-t) {\rm d} t \\
& = \int_{-\infty}^\infty r(\tau) \delta(\tau-t) {\rm d} t \\
& = r(\tau) \int_{-\infty}^\infty \delta(\tau-t) {\rm d} t \\
& = r(\tau)
\end{align}$$

Step 1: Insert (2) into (1). Note that the inner integration variable is a new one, different from $\tau$. I call it $t$.
Step 2: Pulling the first exp inside the integral, combining the exps.
Step 3: Changing integration order (PSD and ACF are absolutely integrable), pulling out what does not depend on the inner integration variable $f$.
Step 4: Using the fact that the inverse Fourier of a constant is a delta (think $\tau-t$ as one variable here, then it's the inverse Fourier of 1).
Step 5: Using the sifting property of the delta.
Step 6: Moving out the constant term
Step 7: Area under the delta is one.

Of course, step 4 is the critical one. If you don't buy it, you need a different, more fundamental/mathematical approach. This is more the engineering point of view I'm presenting here.

Regarding your reply with $\tau$ vs. $t$: What you say is not true. See, we're computing $r(\tau)$ via $\int_{-\infty}^\infty {\rm e}^{\jmath 2\pi f\tau} S(f) {\rm d}f$, which means it's an integral over frequency and the integration kernel depends on $\tau$. The function $S(f)$ does not depend on $\tau$ of course. Now, you are replacing $S(f)$ by the inverse Fourier transform of the autocorrelation, which is an integral over time. But it's a different time variable, which I called $t$. It must be, since it if were $\tau$ it would mean that $S(f)$ somehow depends on $\tau$.

Another way to think about it: The variable $\tau$ is the independent one. Don't forget that all our integrals are definite ones (from $-\infty$ to $\infty$). We sometimes drop that for laziness, but I added them now to be more clear. Now, the integration variables on the right-hand side are the ones we integrate over, hence the result cannot depend on it. Imagine something like $x(\tau) = \int_{-\infty}^\infty \int_{-\infty}^\infty X(f) g(\tau) {\rm d}f {\rm d}\tau$. This does not make sense as the right-hand side integrates over $f$ and $\tau$ (the result is a number) whereas the right-hand side depends on $\tau$. This is why we need a new time variable.


up sample and down sample Let's say that I have a sampled signal x[n], it is being, in this exact order, up sampled by 2, down sampled by4, up sampled by 4 and down sampled by 2 to produce y[n].

It seems to me that it should be pretty self evident that since we up sampled the signal by 2 and down sampled it by 2, then up sampled it by 4 and down sampled by 4, I should just get the original x[n] back.

Am I right?

So the real question is, can the various up/down sampling pieces be readily swapped?

The downsampling by a factor of 4 can introduce aliasing (because you upsampled only by a factor 2 in the first stage). If this is the case, it cannot be undone by the following upsampling operation. So in general this system will not produce at its output the original input signal $x[n]$.


How to improve the reliability of KF estimation? Sorry for boring you my dear friends.

Rcently I'm taking advatage of the state estimates X, outcoming from a Kalman Filter, with the their statistical information, i.e. Matrix of covariance P, to implement other studies.

Theories and practical studies have shown that the quality of KF state estimates is severely depending on the 'correctness' of system modeling, while unfortunately the true situation is unknown. Generally, the KF is either too optimistic or too pessimistic about the state accuracy, i.e. the state covariance matrix P.

The question is therefore, is there any method to reflect the variation of state estimates or the reality on the state covariance matrx P?

Thank you in advance for taking a look!

For very long times, ive used the KF as a parametrized estimator, training the noise matrices, fitting them standard conditions. Just like a "neural net" on the covariance noises. And you have killer results ;).

And the method?... I've applied any standard nonlinear, gradient descent optimization algorithm, such as Newton Raphson, Levenberg Marquardt...

Sweet, simple, unexpected, and excelent results :). Though unpublishable xd.

If you have data, we can try them...


What does the term phase actually mean? I am new to seismic data processing and I really have no understanding of the term 'phase'. So, could anybody give me a simple explaination of the term?
Thank you

I think wikipedia can more than adequately answer your question.

http://en.wikipedia.org/wiki/Phase_(waves)

However, in brief, the phase term describes the relationship between a waveform and a fixed reference point in time. For example the sinusoid $sin(\omega t)$ is zero at $t = 0$ whereas $sin(\omega t - \phi)$ is zero at $t = \phi$. $\phi$ could be referred to here as the phase offset.


Sound Source Location I have four microphones which I want to use to locate the direction of an audio source. Is there any Python module that implements an algorithm that, given the recordings of the four microphones, will give me a vector pointing to/angles of the audio source?
Thanks in advance

I wish I could say that there are 2 different approaches to determining the direction of a source relative to some phase center of an array but that wouldn’t be correct.

Firstly I don’t think there is a Python module that does it generally because there are no algorithms that are completely general.

Having said that there are not 2 general classes of algorithms, I’m going to ignore what I just said.  There is beamforming and there is cross correlation.

Beamforming seeks to use coherent gain to determine direction.

Cross Correlation seeks to determine the relative delays between the sensors and then pick the best solution.

Both categories are based on the principle that a propagating wave can be described by a wave front which is a surface in 3 or 2 dimensions of stationary phase.  In a homogeneous 3 d free space, a point source the wavefront  consists of expanding spheres where the amplitude drops off as $1/r$. Far enough away from the point source, the wavefront is sufficiently flat to assume a plane wave.

So in a real problem you have to make some assumptions about the propagating wavefronts.  We can also have non radiating field components near the source.

The most typically standard approach is to assume plane waves from point sources impinging on your sensors.  There are other kinds of assumptions that can be made like near point sources having wavefronts that have some curvature but you only have 4 sensors.

In beamforming we apply a set of delays to each received signal based on the assumption that it came from some direction and stack them up and sum them.  The set of delays from the correct direction would produce the greatest coherent gain. This means that you have to form a lot of beams for a lot of proposed directions.  Being off a little bit means that the estimated direction would be off a little as well. Beamforming is expensive in real time. Old fashioned RADAR solved the problem by scanning with a rotating parabolic antenna. In SONAR, rotating parabolas are mechanically undesirable and forming multiple beams is the better way to go.

The cross correlation method goes after the delays but you have to do all pairs and then you need to figure out which direction produced those delays.

There is something called split aperture beamforming which is a combination of the two.

When you have more than one source the basic ideas mostly hold but there are issues that are best left for a book or 2 or 10.

If your sensors are too close together, there are problems.

If your source is distributed some adjustments will help.

One good aspect of these problems is the Bessel Fourier relationship because like Fourier series in time,  an arbitrary wavefront at a temporal frequency can be described by a sum of plane waves.

There is a very very large literature on direction finding. I sometimes wonder if there are more papers than there are actual arrays but that is getting less true (not actually true) as time goes on.


Wavetables and Partials I am reading Direct Synthesis versus Wavetable Synthesis. The quote below explains code I have read but I don't get why it's a division by the integer multiple, when multiples increase in frequency. Doesn't that mean we need more samples (higher sample rate mapping to more samples for same accuracy) to deal with the higher frequencies, resulting in a larger wavetable.?

  If you play a loop of a wavetable that is exactly one period at the fundamental frequency then the harmonic partials will fit within that loop. That is because the period of the harmonic partials are the fundamental period divided by N.

Harmonic partials have frequencies $f_i=if_1$, $i=2,3\ldots,$ where $f_1$ is the fundamental frequency. So the frequencies are not doubled, but they are integer multiples of the fundamental frequency. But if you talk about periods, then you have $T_i=1/f_i$, i.e. $T_i=T_1/i$, $i=2,3,\ldots$. So the period of the $i^{th}$ harmonic partial is the fundamental period divided by $i$.


What data should I use to test an FFT implementation, and what accuracy should I expect? I'm involved with an effort to implement an FFT algorithm, and am curious what the recommended advice is for the input test data to use -- and why! -- and what accuracy to expect.

On test inputs, I've found a little guidance in old Usenet posts that I'll post as an answer, but it's just one person's suggestions without a lot of justification -- I haven't found anything that looks like a solid answer.

On accuracy, Wikipedia says that the error should be O(e log N), but what's a reasonable expectation in absolute terms?

Edit to add: The actual tests are in a form where I have stored arrays of input data and pre-computed "reference" output data to compare to, so I don't necessarily need something with a closed-form solution.

If you want to verify an FFT algorithm for correctness, in the sense that it performs the desired function that has the known properties of the discrete Fourier transform, then you can use the approach proposed in:

  Ergün, Funda. (1995, June). Testing multivariate linear functions: Overcoming the generator bottleneck. In Proc. Twenty-Seventh Ann. ACM Symp. Theory of Computing. (p. 407–416).

The above paper is referenced by the makers of FFTW as their method of choice for verifying that a particular FFT implementation does what it should. The proposed technique distills the function into three main components that are verified with separate tests:

Linearity: The DFT (along with its other cousin transforms in the Fourier family) is a linear operator, so for all values of $a_1, a_2, x_1[n], x_2[n]$, the following equation must hold:

$$
\mathrm{FFT}(a_1 x_1[n] + a_2 x_2[n]) = a_1\mathrm{FFT}(x_1[n]) + a_2\mathrm{FFT}(x_2[n])
$$

DFT of the unit impulse: A time-domain signal equal to the Kronecker delta function is applied to the input of the FFT algorithm and the output is checked against the known DFT of the unit impulse function (it transforms to a constant value in all output bins). If the FFT algorithm provides an IFFT, it can be tested in reverse to show that it yields the unit impulse function again.
Time shift: Two sets of data are applied to the input of the FFT algorithm; the only difference between the two in the time domain is a constant time shift. Based on the known properties of the DFT, this should effect a known linear phase shift between the two signals' frequency domain representations, where the slope of the phase shift is proportional to the time shift.

The authors of the paper assert that these tests are sufficient to validate the correctness of an FFT implementation. I've not used this technique in the past, but it does seem to make sense, and I would trust FFTW's authors (who have produced a great piece of free software) as credible authorities on good approaches to the validation problem.


From IQ signal to FM modulated carrier, how its done? I am trying to comprehend GMSK modulation/demodulation. IQ samples are being phase keying is such a way that after multiplying them by complex carrier and sum both I and Q together, there is FM modulation. I don't understand how its being done. I guess it must be trigonometry but can someone explain me how phase changes in IQ result in frequency changes in HF signal?
Thanks in advance

The solution to this is realizing the relationship between frequency and phase:
$$f(t) = \frac{d\phi(t)}{dt}$$
Where $f(t)$ is the instantaneous frequency as a function of time, and $\phi(t)$ is the instantaneous phase as a function of time. So by such definition, a linear change in phase represents a step change in frequency. This is exactly how we can use vector modulators as frequency translators (and equivalently single-sideband modulators). A rotating phasor on the IQ plane, if rotating at a constant rate represents either a positive or negative frequency depending on the direction of rotation. So with FSK using two frequency symbols we get positive and negative rotating phasors. The frequency realized is the rate of that rotation (cycles/sec).
Starting with MSK, which is a subset of FSK, and specifically the minimum frequency shift between the two FSK tones such that the two signals are orthogonal, for this we rotate the phase 90° over one symbol duration, such that the direction of rotation determines if a one or zero was transmitted. The reason for 90° is this results in the minimum spacing of the two frequency tones as first described (hence MSK). If we rotated further within the same time of a symbol period this would result in a higher frequency and thus increase the spacing between the tones (as in generalized FSK). Below shows the phase diagram of rotation versus symbol "1" or "0". If we continued to send a long series of either symbol, it is quite clear that we are transmitting either a positive or negative frequency (relative to a central carrier frequency, so basically one of two possible tones: FSK).

When we transmit a series of symbols the frequency shift keying that is occurring is clear from the relationship between phase and frequency:

The sharp transition in frequency as a step is what creates a high spectral occupancy (sharp edges create high frequencies, consider all the spectral components of a square wave). By rounding these edges we minimize the spectral occupancy, and Gaussian Minimum Shift Keying specifically replaces the straight phase slopes with a Gaussian filtered transition; the phase specifically transitions as an integrated Gaussian pulse (thus the frequency follows a Gaussian transition- see diagram below). So the phase starts out slowly at the beginning of the symbol, reaches maximum rate mid way and then terminates slowly at the end where a possible otherwise abrupt direction change could occur. You can make a very simple (no multipliers so can be done in a micro-controller) GMSK modulator using a Gaussian Filter on the Frequency Control Word input to an NCO which is otherwise connected to the data stream (and similarly with a VCO in the analog world). (This is actually what I was getting at with this question: Simplest All Digital GMSK Modulator which can be combined with this solution for the simple no-multiplier Gaussian filter: Gaussian FIR filter with no multipliers?).  The GMSK phase and frequency diagrams below further and most immediately shows how this is a form of FM modulation.

GMSK is also commonly implemented as partial response signals where we overlap consecutive pulses with intentional inter-symbol interference which can later be completely undone (with Viterbi decoding for example), and thus further improves spectral efficiency since we can send more data in the same amount of time (at the expense of receiver complexity).  The strategic use of the wonderful Laurent Decompositions significantly simplifies those receivers— but that is another topic for another day.
Further FSK, MSK, GMSK are all constant envelope modulations so have the noted advantage that we can run the transmitter PA into saturation which offers significant power efficiency and without regard to immunity against intention jammers can hard limit the input instead of using an AGC (both of particular interest to low cost battery operated devices)


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

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

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

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


Comparison between two Digital Filters Suppose you have two designs of IIR filters with different specifications which are meant to be implemented on a DSP. On what parameters will one compare the two filters so that we can choose which filter is better (like Speed, Memory or Efficiency)? And how do these parameters vary filter to filter?

For example-
I have 2 filters, A and B:
Filter A is a 4th order Low Pass Elliptic filter and filter B is a 70th order Butterworth Low Pass Filter and both can do the required job satisfactorily.
Which filter should be implmented on the DSP and why? (Not taking into account the difficulties one might face in implementing a different algorithm)

The comments actually clearly put forward the complexity of a serious answer. But I'ld like to go simple and put mine for your example filters.

As you've stated, if two IIR filters of $70$th order and $4$th order (of whatever type) can satisfactorily do the job for you, then you should select the $4$th order filter because it will take much less processing power and processing time to compute its output.

That being said, those two filters (the elliptic and the Butterworth) have different frequency domain characteristics. They both have their own magnitude and nonlinear phase responses.

The Butterworth filter will have a monotonic frequency response in both of its pass and stop bands therefore in the design process it will be an overkill to specify a maximum allowed error ripple (from the ideal response) which will be satisfied at the cutoff freqencies but otherwise will be oversatisfied on the rest of the band due to monotonicity of the magnitude curve. That's why the Butterworth filter requires a much larger order for the same given filter specs.

The elliptic filter is known to have equiripple oscillatory characteristics on both the pass and stop bands which distributes the error energy accross all of the bands, therefore reducing the maximum (peak) of the allowed error ripple for a given filter order or minimizes the order of the filter for a given allowed peak ripple and transition width.

Since you don't want to consider the effects such as coefficient quantization or implementation architectures for high order recursive filters, then there isn't much remains to be considered for those two specific filters.


Identifying license plate number of a thief's car from a video Motivation:

I recently had my laptop stolen, with lots of pictures of my family, research work and other very important stuff. Of course the theft has been reported, but I think people here are more competent in image processing compared to the local police station :)

Video Recordings

The stealing was captured by a security camera, original .dav files can be downloaded from here:

ch10.dav
ch13.dav

Converted .avi files can be downloaded (or played in browser) from here:

ch10.avi
ch13.avi

Description

In the file named ch10 (full name: HCVR_ch10_main_20160501143700_20160501143947) the thief's car enters the frame from the right at 0:48 . In the file named ch13 (full name: HCVR_ch13_main_20160501143811_20160501144205), in the first 2 minutes, the thief's car enters the frame from the right, makes U-turn, steals my stuff,  makes a U-turn again and go.

What has been tried?

I tried cropping a few frames around the license plate, stretch so that the plate will be on same position and size on every frame, and averaging. The result wasn't better from the original frames :(

How could I get the number?

Averaging may not be the best option, due to motion compensation in video frames: one number (say $8$) may be extrapolated for another (say $3$) on a different JPEG block. It might be interesting to dig into the compressed video format, check motion vectors, but this requires a little bit of hacking.

Some software which can help:

Amped FIVE (used by the local CSI here, I learned it directly from them), with a live demo
Focus Magic, with a review here.

Unfortunately, while your videos are pretty good, the problem sounds quite difficult. For alight side, I'd like to share the following picture, that tells about issues in image forensics:


Implementing the DWT I have been given the task to implement the 5/3 CDF transform for image compression.

Given that the impulse response for the low and high pass are:

$h_1 = [-0.5, 1 ,-0.5]$ (High Pass)
$h_2 = [-0.125, 0.25, 0.75, 0.25, -0.125]$ (Low Pass)

Say I want to start by computing the horizontal 1D DWT for an image. Now given that the pixels of a row obtained are $x[n]$, we can compute the convolution $y[n] = x[n]*h[n]$. At the beginning of the convolution process we have partial overlap, then full overlap and then partial overlap again. But due to the last part of the partial overlap we obtain extra values.

For the high pass filter I end up with 2 more and for the low pass 4 more. This means that the values obtained for the row pixels are now slightly larger. What am I do about this with respect to programming? Do I eliminate the last few values so $y[n]$ will match $x[n]$ in size?

Another very important question, I have implemented this regardless but I face an even greater problem. While I am fairly new to python, the algorithm I implemented is very very slow. Is there are faster way to implement this rather than using convolution? I will refrain from uploading the code here, because I am not sure if my problem is applicable to this forum or stack overflow. I shall upload it upon request.

Treating boundaries correctly is non trivial for wavelets and filter banks in general. Symmetric/anti-symmetric filters (such as with biorthogonal wavelets as the 5/3 or the 9/7) help a little. Yet, with dyadic sub-sampling, without playing too heavily on polyphase matrices and symmetries, I would suggest the following procedure:

choose a level $L$ for the iterated decomposition
find an image size that can be divided by $2^L$
extend the original image to that size (symmetry/anti-symmetry, or zero-padding but not the best)
extend the smallest filter with zero-padding to match the longest filter size, so both filters have the same padded length.

Cleverer/sparser exist, but they are more involved. Now the number of coefficients ought to be similar. As for the Python, this is more "language-related", and not adapted to SE.DSP. IN 1D, the lifting scheme can improve the computation burden by a factor of 2 in terms of basioc operations.