State Estimation by steady state Kalman Filter I am working with a discret Kalman filter on a System

$x_{k+1}=A_k x_k+B_k u_k+\omega_k$

$y_k=C_k x_k+\upsilon_k$

$E[\omega_k\omega_k^T]=Q$

$E[\upsilon_k \upsilon_k^T]=R$

I have estimated the State from the available noisy $y(k)$, which is generated from the same system state equations with Reference Trajectory of the state. Then I have tested it with a wrong initial state $x_0$ and a big initial covariance (simulation 1). I have noticed that the KF works very well, after a few steps the gain $K$ quickly converges to a very small value near zero. I think it is maybe caused by the Process noise $Q$. I have set it small because the $Q$ stands for the accuracy of the model.

Now I want to modify it to a Steady State Kalman filter. I used the steady gain from simulation 1 as constant instead of the calculation in every iteration. Then the five equations can be simplified to one equation:

$\hat{x}(k+1)=(I-KC)A \hat{x}(k)+(I-KC)B u(k)+K y(k+1)$

I want to test it with same initial state and covariance matrix as the one in simulation 1.
But the result is very different from reference trajectory, and even the result of simulation 1.
I have tested it with the covariance matrix $p_\infty$, which is solved from the discret Riccati equation. And the use the equation:

$k_\infty=\dfrac{p_\infty C'}{C p_\infty C'+R}$

This neither works. I am wondering:

How should I apply Steady State Kalman Filter and how should I set the initial state for it?
Is Steady State KF only used for scalar system?
Should I use it with LQ-controller or some others?

Eddit

I want to compare the performence between standard Kalman filter and steady state Kalman filter like this:
firstly, I set a wrong Initial State: $X_0=[5, -6, 0]^T$ but in fact it should be $X_0=[0, 0, 0]^T$. And I set the $P_0$ like a diagonal Matrix, on the main Diagonal are $25, 36, 0$. I get the steady gain $k_\infty$ and  $p_\infty$ from solving the discret Riccati equation

$P_\infty=AP_\infty A^T-(AP_\infty C^T)(CP_\infty C^T+R)^{-1}(CP_\infty A^T)+Q$

This can be solved by using Matlab Toolbox. I adjust the Koefficient $Q$ ($R$ is constant and available) until the "real" Kalman gain same as $k_\infty$ and at the same time standard Kalman filter works well. Then I think $Q$ in this time the best one. The result is good

And now I don't know, how to set the parameter in steady state Kalman filter.
I set  $X_0=[5, -6, 0]^T$ and $R$ same as ones in the standard kalman filter and use $k_\infty$ as constant Kalman gain. But it works with same parametmer not well.

I like to know, what's wrong with my design.
The plot is below:

I've just tried a toy example of doing what you're suggesting and it works fine. The plot below shows how the standard Kalman filter (red line) and the steady-state Kalman filter (green circles) outputs are the same, except for some minor differences at the beginning.

The bottom plot shows the evolution of the Kalman gain coefficients for the standard Kalman filter.

Without seeing more detail of what you're doing, I don't think I can help much.  Note that I'm using SISO example, which may not be what you're after... but I haven't ever seen reasons for MIMO models to work differently than expected.

R Code Below

#27879 and 29963

T <- 100

#H = [3 -3 1]
H <- matrix(c(1,0,0),c(1,3))

#F is just a tapped delay line
alpha <- matrix(c(0.9,0,0,0,0.9,0,0,0,0.9),3,3)
A <- matrix(c(3,1,0,-3,0,1,1,0,0),3,3)
B <- matrix(c(1,0,0),c(3,1))

sigma_v <- 1
v <- rnorm(T,0,sigma_v)
sigma_w <- 0.01
w <- rnorm(T,0,sigma_w)

x <- matrix(c(0,0,0),c(3,1))
z <- rep(0,T)

for (t in 1:T)
{
  z[t] <- H %*% x + w[t]
  x <- A %*% x +v[t]
}
par(mfrow=c(2,1))
plot(z,type="l", col="blue", lwd=5)
title('Measurements and both Filtered Measurements')

Q <- matrix(c(0,0,0,0,0,0,0,0,sigma_v^2),3,3)
R <- sigma_w^2

library("MASS") # For pseudo inverse ginv()

xkm1km1 <- matrix(rep(0,3*T+3),3,T+1)
xkkm1 <- matrix(rep(0,3*T),3,T)
K <- matrix(rep(0,3*T),3,T)
Pkm1km1 <- matrix(c(1000,0,0 ,0,1000,0, 0,0,1000),3,3)
zhat <- matrix(rep(0,T),c(T,1))

for (k in 1:T)
{
  xkkm1[,k] <- A %*% xkm1km1[,k]
  Pkkm1 <- A %*% Pkm1km1 %*% t(A) + Q
  K[,k] <- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
  xkm1km1[,k+1] <- xkkm1[,k] + K[,k] %*% (z[k] - H %*% xkkm1[,k])
  Pkm1km1 <- (matrix(c(1,0,0,0,1,0,0,0,1),3,3) - K[,k] %*% H) %*% Pkkm1
  zhat[k] <- as.numeric(H %*% xkkm1[,k])
}
lines(zhat, col="red")

xkm1km1 <- matrix(rep(0,3*T+3),3,T+1)
xkkm1 <- matrix(rep(0,3*T),3,T)
zhat <- matrix(rep(0,T),c(T,1))

for (k in 1:T)
{
  xkkm1[,k] <- A %*% xkm1km1[,k]
  K[,T] <- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
  xkm1km1[,k+1] <- xkkm1[,k] + K[,T] %*% (z[k] - H %*% xkkm1[,k])
  zhat[k] <- as.numeric(H %*% xkkm1[,k])
}
points(zhat, col="green")

plot(1:100, K[1,], type="l", ylim=c(-1,1))
title('Kalman Gain coefficients')
lines(1:100, K[2,], col="red")
lines(1:100, K[3,], col="green")


Quick question about difference between two implementations of MFCC I'm trying to learn more about voice recognition, and I was looking this library here and noticed there were two implementations of MFCC (one called mfcc and another called mfcc2). I ran both of these on the same data and got different results, and was wondering why there are two different version of MFCC that give two different results. Thanks.

At a glance - the number of filters is different, the first algorithm has it equal to 48, while the second one - to 32. But even these small things aside, MFCC computation is basically a pipeline of algorithms with 4 or 5 stages. And each stage has some input parameters varying which you can get slightly different output. Given that each stage is implemented properly these differences have negligible effect on the recognition result.


Classifying animated GIFs as animation or film Is there any known way of differentiating between frames of real footage and animated content?  It seems like there's probably range of colors more likely to appear in a histogram for real footage, or maybe a level of color variation more likely to be in animation.

The case I'm looking at is converting GIF images to h264, so there's a max of 256 colors to work with.  I'm mostly trying to figure out when I should use the animation profile so the output is screwed up as little as possible in the transcode.

It's a studied problem. Here are some relevant papers:

Detecting cartoons: a case study in automatic video-genre classification
A color-action perceptual approach to the classification of animated movies
Detecting Cartoons: Automatic Video Genre Classification


Equivalent 2D mask of moving-average I have the moving-average mask as

   mask = [1 1 1; 1 1 1; 1 1 1];

and then I compute the convolution 3 times

   imageF = conv2(conv2(conv2(originalImage, mask), mask), mask);

I want to know how can I get an equivalent mask to compute the filter with just one convolution

   imageF = conv2(originalImage, equivMask);

Convolution in linear time-invariant system is asociative. So to get the equivalent mask you just need to convolve the kernel with itself twice. This will then then give you a 7x7 kernel:

octave:1> a = [ 1 1 1 ; 1 1 1 ; 1 1 1 ]
a =

   1   1   1
   1   1   1
   1   1   1

octave:2> conv2(a,conv2(a,a))
ans =

    1    3    6    7    6    3    1
    3    9   18   21   18    9    3
    6   18   36   42   36   18    6
    7   21   42   49   42   21    7
    6   18   36   42   36   18    6
    3    9   18   21   18    9    3
    1    3    6    7    6    3    1

octave:3>

Note that the original mask is not normalised - it has a gain of 9 at DC - so with three convolutions you get an overall gain of 9^3.

Depending on what you're trying to achieve though you might just be better off with a 7x7 Gaussian.


Very basic question about how we define frequency in signal processing When talking about general periodic continuous-time signals for which $$x(t + T_0) = x(t)$$ where $T_0$ is the fundamental period we define the fundamental frequency $\omega_0$ as $\omega_0 = 2\pi/T_0$.

The way I interpret this is that $1/T_0$ is the frequency of the signal in cycles per second and there are $2\pi$ radians in one cycle, therefore the angular frequency is $2\pi/T_0$ radians per second. But are there really $2\pi$ radians in one cycle in general?

Take $x(t) = \tan(t)$ for example. Its fundamental period is $\pi$ and using the definition of fundamental frequency above, its fundamental frequency is 2 radians per second. But in the case of the tangent function, aren't there $\pi$ radians in one cycle, as opposed to $2\pi$ radians? This would give us a definition of frequency as $\omega_0 = \pi / T_0$. Or do we refer to the unit circle by convention when talking about "cycles"?

Apologies if this question sounds very basic.

For the tangent function a cycle is indeed done in $\pi$ radians, because you want

\begin{align}\tan(x) &= \tan(x+T)\\
\implies \frac{\sin(x)}{\cos(x)} &= \frac{\sin(x+T)}{\cos(x+T)}
\end{align}
$$
\implies \sin(x+T)\cos(x) - \sin(x)\cos(x+T) =0
$$
$$\implies \sin\big[\left(x + T\right)-x\big] =0\equiv \sin\left(T\right)=0$$
$$\implies T = \pi k, \quad \text{with}\quad k\in \mathbb Z$$

Smallest strictly positive $k$ (i.e. $k= 1$) gives us $T_0 = \pi$. Anyway, we know that.

Regarding the interpretations, let's use $\Omega \ [\text{rad/sec}]$ to avoid ambiguity with $\omega \ [\text{rad/sample}]$ mostly used in the discrete case.  From the definition of the angular frequency $\Omega = 2\pi F$ with the fundamental period $T_0 = 1/F$, as you have already described, we have
$$
\Omega_0 = \frac{2\pi}{T_0}\tag{1}
$$
The interpretation should go from $(1)$ and read as the number of cycles of the periodic signal in $\mathbf{2\pi}$ radians. As examples, look at the plot below,

As seen above,

$\tan(x)$ has $T_0 = \pi$ implying $\Omega_0 = 2\pi/\pi=2 \ \text{rad/sec}$ or 2 cycles in $\mathbf{2\pi}$ radians (or equivalently one cycle in $\pi$ radians)
$\sin(x)$ has $T_0 = 2\pi$ implying $\Omega_0 = 2\pi/2\pi=1 \ \text{rad/sec}$ or 1 cycle in $\mathbf{2\pi}$ radians
$\sin(2\pi x)$ has $T_0 = 1$ (see $x$-axis) implying $\Omega_0 = 2\pi/1=2\pi \ \text{rad/sec}$ or $\mathbf{2\pi\approx 6.2832}$ (see $x$-axis) cycles in $\mathbf{2\pi}$ radians.

See it this way, the ratio in equation $(1)$ is simply stating how many  $T_0$ (i.e. fundamental periods or full cycles) you have in $2\pi$. You have $T_0 = 2\pi \implies \Omega_0 = 1 \ [\text{rad/sec}]$ meaning you have $1$ fundamental period (full cycle) in $2\pi$. Similarly, when $T_0 = 1 \implies \Omega_0 = 2\pi \ [\text{rad/sec}]$ meaning you have $2\pi$ fundamental periods (full cycles) in $2\pi$. Etc. for other cases.


Which book of Proakis has dealt with continuous time signals and systems? I have heard of three books of Proakis namely Communication System Engineering, Digital Signal Processing and Digital Communications. I love the way he writes and would like to read about continuous/analog signals and systems by the same writer. Some of it has been dealt with in his book 'Communication System Engineering' namely Fourier Series and Fourier Transform. I want to read other continuous time topics like Laplace transform by the same author.

  Is there any book by Proakis available which fully deals with
  continuous time signals and systems?

No, not as far as I know. I can recommend you two books dealing with continuous-time systems, the Fourier transform, and the Laplace transform, among other topics:

Signals and Systems, Oppenheim, Willsky: a classic, used copies can be found quite cheap
The Fourier Integral and Its Applications, Papoulis: less scary for an engineer than the title might make you think; it's actually much broader and very accessible; it deals with linear systems in general, low pass and band pass systems, the Laplace transform, etc.


Is it possible to do deconvolution with two data sets that have different sampling rates? I have some terahertz spectroscopy time series data, a reference set with 2048 data points taken every 0.0521 picoseconds, and the sample data set with 544 data points taken every 0.0781 picoseconds. I'm using Matlab to take a FFT of both sets, with zero padding on the sample set up to 2048 points, and I'm supposed to do a deconvolution on those transformed sets on a specific range of frequencies, but because the sampling rate is different for the two, these correspond to different starting points in the data sets, and the lengths won't match up either. When plotting the two DFTs, they only line up properly when accounting for the different sampling rates when spacing the frequency axis.

Is there something I'm doing wrong, or how should I do the deconvolution?

Resample your data to be close to the same rate. Note that the ratio between the time samples is 1.4990403....   if you resampled the 0.0781 data by a frequency ratio of 3/2, you would result with a time error of 0.000033 ps per cycle. From your sampling rate and block length you can determine how much of an impact this would have to see if that scale is sufficient: 0.0521/.000033 =  157.9 samples to slip $2 \pi$ radians. Resampling with higher ratios to support larger data block sizes: For instance 9369/6250 will get you to 1.49904 with a resulting time error of 3.07E-7 ps per cycle, with a block length of .0521/3.07E-7 = 169,650.6 to slip $2 \pi$ radians.

I have not tried what I am about to suggest next so my thought process could be flawed, but it seems that you could feasibly do the simpler 3/2 resampling ratio and compensate for the remaining frequency error with a quadratric phase versus frequency all pass function. This is because the residual frequency error is the error at the sampling rate, and the error at all frequencies will be in proportion of their frequency to the sampling rate (so a linear function over frequency). Phase is the integral of frequency versus time, but I am thinking of the possibility of a compensating quadratic phase (vs frequency) function to reduce the residual error. Note in this case the derivative of phase versus frequency is delay, so I believe that this would be similar to emulating a phase shift with a time delay; in that it will have the correct phase delay as anticipated at one frequency only, and then will have linear error (in phase) as you deviate from that frequency.  By getting our frequencies close with the 3/2 resampling, and assuming sufficient points in the fft such that the frequency spacing is small I could see the possibility that a smaller overall error is achieved with such an approach. This needs to be pursued further mathematically or experimentally but is a thought. If you do anything along this lines sucessfully please share your results! (Or if anyone sees the flaw with this approach that makes it not feasible please illuminate and I will delete this 2nd paragraph). I did post this technique as an additional question: Compensate for residual resampling error with a quadratic phase vs frequency- possible?


Denoising before or after Geometric Transform? I am working on a image which needs denoising (bilateral filter) and a geometric transformation. I am not sure in which order these filters should be applied. Personally I think it is better to use the transformation first because otherwise I would loose details removed by the bilateral filter. But in the same moment the transformation would amplify some of the existing noise wouldn't it?
I tried it both ways and objectively it's hard to tell which one is better.

So which order is right? And how to order filters in general?

The transform will alter the spatial spectrum of the noise, at least for most areas of the image. So denoise before the transform. More specifically, denoise in the coordinate system that the noise was generated in. This assumes, of course, your denoising algorithm is tailored for the untransformed noise.


What are the Segmentation algorithms that can be used to process Satellite images (obtained from Google Earth) to segment forest and barren areas? I have tried Region based Segmentation - through Seed Growing, Watershed algorithm, and Graphcut algorithm, but none have been successful so far. I am using MATLAB.

You could try image saliency or entropy filtering. Entropy filtering calculates the entropy of gray values in the image Image within a rectangular mask. The result is written to the output image. The process is repeated in a sliding window fashion.

In the resulting image, the forest areas should appear more whitish while the barren areas appear black, because they don't contain enough information.

Of course I assume the roads would make a problem. However you could get rid of them using morphological edge finding operations.


intuitive interpretation Fourier Transform of two different Rectangular pulses From my understanding of Fourier transform, a Fourier transform of a signal in time domain will give the different frequency components of the signal in frequency domain, specified by their contributions(energies) in the original.

The Fourier transform of rectangular function is a sinc function

If I take Fourier transform of two different rectangular pulses, with width w1 and w2, with w1>w2, the sinc function of one width w1 is narrow and has large peak as compared to the one with width w2(wider and smaller peak). How do I intuitively understand this behavior?

Note that the peak height of the Fourier transform of a wider rectangular function is proportionally higher.  This is necessary to preserve total energy, as per Parseval's theorem.  However, the energy required to create the sharp edges of a rectangle in the time domain is still out there on the tails of the Sinc function in the frequency domain, roughly the same (but more wiggly).

But when you scale the vertical axis of 2 plots so that they have the same peak height, the envelope of the FT of a wider rect function looks narrower in the plots of that transform.  But that's only in relation to the peak with a scaling applied.


How to go from wav file to spectrogram back to wav file in python? Close to this question but different : Wav To Spectrogram, Back To Wav

I already posted this on StackOverflow(https://stackoverflow.com/questions/54743968/how-to-go-from-wav-file-to-spectrogram-back-to-wav-file-in-python), but I figured this question might be more suited to be here :

Here's my problem. I have some wav files. I want to create their spectrogram. Then I want to do some stuff with the spectrogram and then I want to recreate a wav file associated to the new spectrogram. I am aware that the spectrogram doesn't have the phase information, however, since I have an original file, I can (theoretically) extract the phase from the original wav file.

I want to do that in Python (if possible), do you know a good library that would allow me to do that?

You need to keep the phase information, which isn't included in a spectrogram.  "Spectrogram" is just the magnitude of the STFT output.

So this is conceptually the same as STFT ↔ ISTFT, which is supported in Python now:

STFT
ISTFT

The Examples section in the ISTFT link does exactly what you want, modifying a STFT and inverse transforming to reduce noise.


Purpose of multiplying fft output by Hann window? I'm working with some acoustic propagation models at the moment.
A feature I've noticed in 2 of these models is that the output of the model is often multiplied by a Hann window in the frequency domain immediately before reporting the results in the time-domain.

I've been struggling to understand the idea behind this, as I'm only familiar with window functions applied in the time-domain.

To give some background... these models operate by (to my understanding):

Simulating a theoretical rectangular spectra (with equal phase values for all frequencies).
The models then calculate the loss per frequency as sound propagates through the environment.
Frequency domain and time domain plots are typically returned to view the results in an intuitive fashion (dB loss per frequency).

I am assuming that to convert this back into the time domain (to show the effects of spectral dispersion), the Hann window somehow makes this signal viewable, as otherwise the flat phase values would result in non-sensical time domain signals.

Can anyone briefly describe the DSP theory behind multiplying a flat (rectangle) spectra with a Hann window of the same size?

Below I have shown some relevant matlab code.
This piece of code (as I understand) frequency shifts the signal (phs), multiplies the frequency domain by a hann window (hanwin), then converts back into the time domain fft.
The resulting tlpressd is used to plot the dB propagation loss per time in the resulting plot.
The looping term just runs through all the output ranges of the environment.

hanwin=hanning(nf+1);
hanwin=fftshift(hanwin(1:nf)/sum(hanwin(1:nf)));
phs=2*pi*cfreq*time;
if optout == 1
  for irng=1:nrout,
    pressd(irng,:)=fftshift(fft(pressd(irng,:).*hanwin')).*exp(-j*phs);
    timeout(irng,:)=time+rngout(1,irng)/(c0/1000.);
  end
% pressd = complex pressure field, row = range; col = time
tlpressd=-20*log10(max(abs(pressd),1.e-20));

Notes

FYI, the code above is from Kevin Smith's "Monterey-Miami Parabolic Equation" propagation model.  It converts the resulting green's function solutionsof the model into a one of several visual representations (time vs Frequency, or Frequency vs Depth ,etc)

To give more context (and possibly more confusion).  If I replace the variable hanwin with any time-domain signal of my choice, the end of the code block provides be with the convolution of that signal with the propagation effects the model calculated (i.e. it simulates the effects of propagation on the provided signal).  Its clear that I'm convolving the Hann window with my model output... but I have no idea what the signal theory behind this is.

Can anyone briefly describe the DSP theory behind multiplying a flat (rectangle) spectra with a Hann window of the same size?

Without any reference to the actual model and its purpose, the only thing one can say about this is that due to the convolution theorem, the described multiplication results in some form of a digital filter. That is all one can say about "modulating" the frequency domain with any sort of function.

BUT, based on this excellent exposition of the way parabolic equations can be solved via the Discrete Fourier Transform, it seems that the multiplication is performed in both the spatial and frequency domain to isolate specific parts of the signal "gracefully".

Now, the authors seem to indicate that the multiplication is done to "...avoid aliasing...", but I cannot see how the results of these operations can lead to aliasing when:

One is in full control of everything about the simulation
The first "plane" (the one that is generated at the source) would only contain harmonics form the desired spectrum and be zero elsewhere anyway.

So, I cannot see how extra harmonics are likely to creep in between integration steps to have a need to isolate them.

Therefore, a more likely reason might be to be able to "isolate" (or focus on) only a part of the domain in a graceful way, rather than setting it to zero.

Hope this helps.


Conditions for symmetric and unimodal windows in both time and frequency domains After a lecture on harmonic analysis and time/frequency methods, I reconsidered the Gaussian kernel, defined in continuous time.

It is unimodal and symmetric, and its continuous Fourier transform is another Gaussian, thus unimodal and symmetric (from Fourier transform of a Gaussian is not a Gaussian, but that's wrong!):

As a filter or a window, the Gaussian ensure a monotonic weight decay around its center. The interpretation in the frequency domain is similar: frequencies are monotonically attenuated away from a center frequency.

Such a property allows a straightforward interpretation in both the time and the frequency domain.

For continuous kernels, are there generic characterizations (necessary or sufficient conditions) under which  unimodal and symmetric windows also have a  unimodal and symmetric amplitude spectrum?

A sufficient condition is, that the distribution is spanned by zero-mean Gauss distributions with non-negative amplitude weights. The spanned space contains only symmetric and unimodal distributions and is invariant under a Fourier transform.


camera movement calculation I'm trying to work out how to estimate camera movement between 2 equirectangular photos (static scenes indoor). Both photos are taken with the same camera (Ricoh Theta S) on the tripod. Height of the camera from the floor is known and same for both photos. Pitch and yaw is known as well. There is a line of sight between both placements of the camera. Is what I'm after even possible? Cheers. Kuba.

Is what I'm after even possible?

yes.

Let's start with the simplification your camera was a perfect pinhole camera.

For starters, imagine you pick an arbitrary point object in your scene. You then get a camera matrix $\mathbf C$ that allows your to predict the 2D position $\mathbf y$ in your image of your 3D-positioned object at $\mathbf x$:

$$
\mathbf y = \mathbf {Cx}, \quad \mathbf y= \begin{pmatrix} y_1\\y_2\\1 \end{pmatrix}\in\mathbb R^3,\,\mathbf x= \begin{pmatrix} x_1\\x_2\\x_3\\1 \end{pmatrix}\in\mathbb R^4,\,\mathbf C \in \mathbb R^{3\times 4}
$$

Notice that $\mathbf C$ isn't reversible, ie. there's multiple $\mathbf x$ that explain a $\mathbf y$.

Now, translating your camera by a movement vector $\mathbf m$ is the same as shifting the scenery by $\mathbf{-m}$, which leads to a $\tilde{\mathbf{x}}= \mathbf x - \mathbf m$.

Correspondingly, we get a

$$\tilde{\mathbf{y}} = \mathbf C\left(\mathbf x - \mathbf m\right)\text.$$

Now, since $\mathbf C$ is a matrix, i.e. a linear operator, this leads to

$$\begin{align}
\tilde{\mathbf{y}} &= \mathbf C\mathbf x - \mathbf C \mathbf m\\
 &= \mathbf y -\mathbf C \mathbf m\quad \iff\\
\mathbf C\mathbf m &= \mathbf y - \tilde{\mathbf y}\text.
\end{align}$$

As mentioned above, $\mathbf C$ isn't reversible, which means that from a single photo you can't tell the 3D coordinate of a point without additional info. However, if you can identify three  objects at different points $\mathbf y_{a,b,c}$ and $\tilde{\mathbf {y}}_{a,b,c}$, you should be able to elimate the unknowns from the system of equations.

Also note that your real camera isn't a perfect pinhole camera, and that might be used to actually get a less underdefined system, so you might be able to mathematically derive a unique solution with less than three points.


How to choose the frequency grid when using DFT with numpy (FFT) I'm currently trying to learn how to calculate DFT using FFT with Python. However, I have been struggling to find how should I choose the frequency grid when plotting the final calculated DFT values out of the data samples that I have. Any helpful insights will be appreciated!

It really depends on your application and what you want to see and the actual FFT call you make.  Typically if your time series is real valued, the DFT will be conjugate symmetric with respect to negative and positive frequencies and the magnitude information is redundant.  Let assume your data consists of N samples and N is even and you used a complex DFT,  also assume you use C style indexes, i.e. X[0] is the first element in your array. X[0] is going to correspond to the "DC" bin or a frequency of zero.
X[N/2]will correspond to the Nyquist sample frequency.  Most People with ignore X[N/2+ 1] to X[N-1] (or often X[N/2] to X[N])  The thing to remember is that each bin has a conjugate valued twin.  X[0] and X[N/2] are pairs of a sort but don't follow this rule.  X[1] is paired with X[N],  X[2] with X[N-1] , X[3] with X[N-2], and so forth.  One is the positive frequency and the other the corresponding negative frequency.

If you recall your sampling theory, the Discrete Frequency Spectra is periodic, so you can think of the values X[N/2] to X[N] as either the negative frequencies of the spectrum of the cycle centered at DC or the negative frequencies of the cycle centered at the twice the sample (+) frequency. I typically will reorder the upper frequencies so that they are displayed below so the DC bin is the center point of my plot and Matlab has a function called FFTSHIFT that does that, and it also works if N is odd.

If your time series is complex, there are different possibilities.  If x is a Hilbert transpose, the negative frequencies are zero, so unless you want to verify that, it doesn't make much sense to display zeros.

If x is a signal such as an I and Q pair out of a demodulator, using a FFTSHIFT ordering is generally preferable because the frequency shifted carrier will be in the center.

The fact that the upper half of the output are actually the negative frequencies was something that confused me many moons ago.


Finding filaments in an image I am at the moment working on images such as this one:

What you see are filamentous structures / bundles. Other images coming from slightly different experiments could have more sparse / thick / bright filaments, but the order-of-magnitude and qualitative average image is like the one above.

What I would like to do is to extract as many information as possible from this image and the best thing would be to identify each single filament/bundle or a "skeletonization" of the image, but I am happy already with a vector field indicating the local orientation of the filament so that I can compute its divergence etc. etc.

So the question is, do you know of any program / conceptual tool / framework which may help me identyfyinf the filaments and or their local orientation?

I usually work in Python and any example-code is welcomed, but I would be more than satisfied with a purely mathematical / pseudocode / conceptual answer.

What I did so far:

I tried to use PCA (Principal component analysis). I select a ROI, threshold the image, create the dataset $(X, Y)$ of the pixels with high intensity and compute the covariance matrix using the intensity as a weight. Then I diagonalise it finding the "longer" eigenvecotr and assume that that is the local orientation. Works ok-ish, but I am afraid is very paramtere-dependant and may fail when filaments become too thick or too thin. An example is below. It looks good but I fear is too simplistic an approach.
I tried to implement (in Python) the procedure they describe here. Maybe I am tweaking the parameters wrong or is supposed to work only for very thin filaments, but it works very badly in my case.

I know there are pre-made tools online, but I'd rather write everything from scratch.

So summing up, the big question here is:

What is the best way to find the filaments and/or their orientation at a given point?

Thank you all!

There seem to be a number of hits in the literature from other domains that use various methods to segment the image (areas likely to belong to a filament and 'background') followed by filament identification  / extraction which you might want to have a look at anyway.

I was intrigued by the choice of PCA as a first indicator of extracting the local gradient and so this is a response with the "low hanging fruit" you can perhaps try (and in fact extremely easily) before jumping to the big guns of pattern recognition.

  ...the big question here is:
  What is the best way to find the filaments and/or their orientation at a given point?

Unless you have prior information (or assumptions) about the direction of the flow / filaments, it is impossible to derive their orientation. In other words, you can get some indication about the local orientation of (something that looks like) an edge, but not the orientation of the flow of the filament. You can assume, for example, that there is a particular intensity profile along the body of the filament which could then be used to determine "which part is the 'head' ".

The first thing to do (anyway) is to obtain the image gradient. The image gradient will give you the local slope of the intensity field (a.k.a image). All code here in Octave but only for demonstration. The operators are easily transferable across to scipy.

%Load the image
Q = imread('someImage.jpg');
%Resize it to 1/4 of its original scale, this is not essential but it helps with visualisation.
Qs = imresize(Q,0.25)
%
% Get the gradient
% [gradMag, gradDir] = imgradient(Qs);
% gradMag is the Magnitude, gradDir is the direction in degrees
% Break down gradDir to its components to process and visualise it
X = cos((gradDir/180.0).*pi);
Y = sin((gradDir/180.0).*pi);
%Let's have a look at the image with its gradient superimposed
imshow(Q2);hold on;quiver(X,Y);

The gradient of the (rescaled) image looks like this:

And if we zoom in (to the plot, not the image) around a bundle, we can see something like this:

So, predictably enough, the arrows are pointing away from the top of a ridge that signifies a filament and this can be used, to an extent, to derive its local orientation.

BUT, the bundle's intensity is not constant and the local gradient is not informed by any prior so that it "knows" which way to point.

Now, once you have moved from a scalar (Qs) to a vector field (X,Y), you might find other vector field operators useful as well, depending on what you are trying to do.

Can we do better?

Yes, we can apply a transformation to the intensity values which dampens the soft values and boosts the strong values. This will 'intensify' the ridges and basically result in a more well defined gradient later on. This is a modification of the image's dynamic range and is also known as 'expanding'.

An exponent is a good choice, it effects the following mapping:

Where we simply raise the pixel value to a power and take care of normalisation so that it does not overshoot the typical image limits.

This transforms the original (scaled) image to something like:

And conversely the local gradient now will be more well defined.

Now, if you are trying to detect broader structures which you would then infer (with some certainty) that they are formed by or belong to filaments, then you can do a couple of things. First of all, you could use Curl to identify "turbulent" (or absolutely no turbulence) areas and from there try to make sense of 'circularity' and then bundle together those areas.

The other thing you can do is something along the lines of line or circle detection through the Hough Transform. The Hough Transform is a "variant" of the Radon Transform. The latter is transforming an image from a $x,y$ plane to a $r, \theta$ plane by obtaining line integrals (the magnitude of which is $r$) at angle $\theta$ around the image. This looks like this (on Qs):

[R,xp] = radon(Qs);
imshow(R./max(max(R)));

The Radon transform converts lines (in Qs) to points. A line will produce a high value integral when $\theta$ is parallel to the bearing of the line. Conversely, arcs (and circles) are converted to arcs in the $r, \theta$ plane and if you follow intense arcs along the $\theta$ direction (corresponding to $x$ in a typical image) you will be able to identify the curled ridges that are formed by the filaments.

This is used by the Hough Transform which essentially operates on a thresholded version of the initial image (e.g. Qs) to detect lines and circles by looking for points with specific spatial relationships in the $r, \theta$ plane.

Hope this helps.


s-Domin or z-Domain - What to Use for Mixed systems I recently had to deal with a power electronic system where I had to implement a dynamic model of a power converter in order to design a suitable controller for that converter. The control will be implemented on a DSP which reads in the digitize values (voltages and currents) from some ADCs. The plant itself consists of a PWM modulator and a 4th-order LC filter. Since ADCs and control are digital systems, I can describe them using the z-transform, whereas the LC filter is an analog system which I usually describe using the Laplace-transform.

The question which now arises for me is: should I discretize the time continuous transfer functions of the filter and work out the whole system dynamics using z-transform, or do I describe the digital controller in s-domain and describe the whole system as a time continuous system?

During my studies I always worked with Laplace transform. I heard that there is z-transform too, but I never really worked with it. This was usually working fine, however, in most of the papers I find the authors are using z-transform and I don't really know when I should use one or the other for mixed systems. As far as I know, it doesn't really play a role in which domain I work as long as the frequencies I'm concerned about are "low enough" compared to the sampling frequency of the system. Is this true?

Another problem I face: when I'm working with these dynamic models in s-domain in matlab and let the program calculate transfer functions, the models which result exhibit insanely high orders (e.g. order 100 for a system which should have order 7). I know that this is numerical issues and that most of those system poles are cancelled out by zeros laying extremely close to the poles. There is a matlab command ("minreal") which eliminates some of those pole-zero combinations but not all. Could I reduce this issue by switching into z-domain? If not: how trustworthy is the "minreal" command? Do I risk loosing important information if I apply this command on every transfer function?

First choice :

Convert the Laplace transform of your process to the Z-domain using the ZOH method as it models your DAC. In your case, your DAC is a PWM.

Second choice :

Work in the Laplace domain but take in account the delay due to your PWM, which should be T/2 where T is your PWM period. To model the delay, you can use the Padé approximation.

I personnaly used the second technique. It was easier, more intuitive. But both techniques will work.

For the numerical issues :
The order of your transfer function might get large when including the controllers, I used minreal too, but you can convert your tranfer functions to ZPK, this solves many of the numerical issues you might have.

Edit : As Tim Wescott mentioned in the comments. The ZOH modelization is accurate when using a real DAC. Not so for a PWM. It's better off to work in the s-domain and try to model the delays (sampling, PWM, computation, etc.) as accurately as possible.


Power Spectrum and Power Spectral Density From signal theory we know that a very useful representation of some power signals is that of its power spectral density, whose curve represents how the total power of the signal is distributed at all frequencies. For instance, it is a couple of dirac pulses for a sine wave of frequency f0 because its power is concentrated at frequencies f0 and -f0.

But in some different situations I have seen that we often use not exactly the power spectral density representation, but the power representation.

EXAMPLE 1: Measurements with Spectrum Analyzer

If we use a spectrum analyzer, we see in abscissa the different frequencies, and in ordinate the power of the signal. For instance, this is a measurement I have done:

I have sent to the instrument a sine wave of -13dBm of power, and its total power is represented in ordinate of the peak. What does it mean from a math point of view?

I think this representation has no sense: the ordinate should be the power spectral density, not its integral in a neighborhood around the frequency considered!

Example 2: simulink

This is the representation of the power spectral density of a certain noise:

This is the representation of the power of the same noise:

What is the link between them? Power should not be function of frequency (it is the integral on all frequencies of power spectral density), power spectral density yes! So what does it mean that we see the "power spectrum"?

What is the link between them? Power should not be function of
  frequency (it is the integral on all frequencies of power spectral
  density)

As far as a Spectrum Analyzer is concerned, to display the power spectrum in dBm, what it does it to take the FFT and compute square of magnitude $P_n = 20  \times log_{10}(|X(f_n)|^2*1000/100)$ with $f_n$ being the discrete frequencies separated by resolution bandwidth (see it is set to 1kHz for this case). I multiplied by 1000 since unit is in dBm and I divided by 100 because to get to power, I have to pass voltage $X(f_n)$ through a resistance of $50\Omega$. ($\frac{V^2}{2R}$).

To derive power spectral density from power spectrum, I am sure it is just dividing power in each bin by RBW. $P_n/f_{rbw}$. Though I am not very sure whether it is the (dBm value divided by rbw) or (dBm value of linear power value divided by rbw).

  I think this representation has no sense: the ordinate should be the
  power spectral density, not its integral in a neighborhood around the
  frequency considered

Well, in this case it is the latter. I think that is how all commercial spectrum analyzers represent spectrum (if you choose units as dBm or Watts). The comment by MBaz correctly mentions this.

Note - The final plot will be slightly different because your FFT size output will be different from number of trace points used to display on the screen. In order to combine several FFT points into 1 trace point on screen, a detector logic is used.


Overview of regression methods between the data I am trying to find different approaches (classic and modern) in order to find a regression function between data (in time domain) with high-frequency noise.
Consider the result of this code (in Matlab):

fs = 10000;
t = -1:1/fs:1;
target = tripuls(t,50e-2);

noise = target + randn(1, length(target));

plot(noise,'o')
hold on
plot(target,'r')

I would like to find methods that, starting from noise find a function similar to target (for convenience, I attach here an image).

A classic solution can be to implement a low-pass filter. Did you know other (and maybe better) way to reach the expected result?

Thanks

There are multiple ways to recover the signal from noise. As you have suggested one such method, the low pass filtering which is a classical signal processing method.

The linear signal model that you have considered is $y[k] = f[k] + b + n[k]$, where
$y[k] =$ output signal, $f[k] = $ true signal, $b=$ noise offset and $n[k] = $  zero mean additive noise

If you want to look from a data point of view the methods of interest would be

Nearest Neighbor Approximation (or local average) : This method is similar to low pass filter but here you can choose a suitable window, $$\widehat{f}_{NN}[k] = \frac{1}{2L+1}\sum_{i=k-L}^{k+L} y[i] - \hat b$$ Here $\hat b$ can be obtained by an initial pilot sequence of the zero signal.
Basis Representation followed by Ridge Regression : With suitable signal class assumptions you can represent $f[k] = w_1 \phi_1[k] + \ldots + w_p \phi_p[k]$ - known as the basis expansion (In the example you have described a triangular wavelet basis would work). In this setting the goal is to learn the weight parameters $w_i, \ldots, w_p$ subject to the prior knowledge on noise, that is the mean and variance of $n[k]$ is known apriori. In case the prior information are unknown you can use the method used to estimate $b$, explained in method 1.
Alternatively, There is an equivalent formulation in convex optimization known as ridge regression which allows for regularization cost and no assumption on the noise parameters.

These methods work out mostly with a training dataset and hence mostly does non-linear (filtering in) signal processing.

Remark and Disclaimer : Please feel free to correct any mistakes in the arguments presented above, as I am an amateur in the area of modern methods of signal recovery.


MIMO channel capacity being complex? I'm considering the following simple system model:

$$\mathbf{y}=\mathbf{Hx}+\mathbf{n}\quad\textrm{with}\quad \mathbf{n} \sim \mathcal{CN}(\mathbf{0},\mathbf{I})$$

The resulting MIMO channel capacity for known and constant $\mathbf{H}$ should be:

$$C = \log_2\det\left(\mathbf{I}+\mathbf{HQH}^H\right)$$

where $\mathbf{Q}$ denotes the covariance matrix of $\mathbf{x}$.

In my opinion the channel capacity should be real valued. But if I simulate it in MATLAB and compute the capacity I obtain a complex capacity. I generate a random complex $\mathbf{H}$ since $\mathbf{H} \in \mathbb{C}^{n \times n}$ and chose a random covariance matrix $\mathbf{Q} \in \mathbb{C}^{n \times n}$ which is symmetric and positive semi-definite.

Am I doing something wrong in MATLAB or do I violate and conditions, constraints or rules?

My Matlab code is:

%% Init
N = 2;
% Generate complex random channel H
H = sqrt(1) * (randn(N,N)+1i*randn(N,N));
% Generate complex, symmetric positive semidefinite Q
Q = [1,complex(3,-7);complex(3,-7),13];
%% Check if Q is a valid covariance matrix
% Check if Q is symmetric
sym = issymmetric(Q)
% Check if Q is positive definite
psd = all(eig(0.5*(Q+Q'))>=0)
%% Compute capacity
C = log2( det(eye(N,N) + H\*Q\*H' ) )

The way you generate the covariance matrix is not right. I would either use Matlab's gallery, or generate a random matrix as follows:

A=rand(N)+j*rand(N);
Q=A*A';

To check if a matrix Q is positive semidefinite, you can look at the eigenvalues and see they are all nonnegative. You can use eig(Q) in Matlab.


Subgradient Method for K-Means Like Problem I want to solve k-means-like problem. But at first, I wonder that we can apply subgradient method to k-means problem.

We want to find the optimal clustering $S = (S_1,S_2,...,S_k)$ such that

$$ {\operatorname{arg\,min}}_S  \sum_{i=1}^{k} \sum_{\mathbf x \in S_i} \left\| \mathbf x - \boldsymbol\mu_i \right\|^2 $$

You should use something like Robust Convex Clustering:

Assume that we are given $ n $ data points and each data point
is described by a $ d $ dimensional feature vector. We collectively
represent the data by a data matrix $ X \in \mathbb{R}^{d \times n} $.
The convex clustering method clusters data points into groups via solving
a convex optimization problem:

$$ \arg \min_{P} \frac{1}{2} \left\| X - P \right\|_{F}^{2} + \alpha \sum_{i < j} {w}_{i, j} \left\| {P}_{i} - {P}_{j} \right\|_{p} $$

Where $ \left\| \cdot \right\|_{p} $ is the $ {L}_{p} $ norm, $ P \in \mathbb{R}^{d \times n} $ is the matrix consisting cluster centers and assignments, $ {P}_{i} $ is the $ i $ -th column of the matrix $ P $ and represents the centroid of cluster that the $ i $ -th data point is assigned, $ \alpha $ is a non-negative regularization parameter that controls the number of clusters. $ {w}_{i, j} $ is the weight between $ i $ -th and $ j $ -th data point specified by the user. Clusters are implicitly given by $ P $ matrix and $ {P}_{i} = {P}_{j} $ means the two data points are in the same cluster. When $ \alpha $ is zero, each data point forms a cluster since $ P $ equals to $ X $ when Eq. (1) reaches the
optimum. As $ \alpha $ increases, $ \left\| {P}_{i} − {P}_{j} \right\|_{p} $ becomes smaller for all pairs so that the whole function can reach the minimal value. Hence more data points will be grouped into a cluster.

Another great reference is Splitting Methods for Convex Clustering.

As the formulation you chose is not convex and of top of that I don't think anyone has ever calculates its Sub Gradient.

Update

I encountered another simple yet effective of K-Means - Revisiting K-Means: New Algorithms via Bayesian Nonparametrics.
It formulates the problem in as, intuitively, the Lagrange form of the K-Means with continuous parameter $ \lambda $ instead of the discrete parameter $ k $.


Understanding circular convolution Let $a$ and $b$ $\in \mathbb C^N$ and $a[k] = a[k \mod N]$ (same with b). Then the circular convolution of $a$ and $b$ is defined by

$$ (a * b)[n] = \sum_{p=0}^{N-1} a[p] b[n-p].$$

I have a problem understanding the flip operation on the convolution. Let $(Ja)[n] = a[-n]$, then

$$ J(a*b)[n] = \sum_{p=0}^{N-1} a[p]b[-n-p].$$

so we first apply the convolution and then we flip the n - is this correct? Or do we flip the whole index of $b$, meaning that

$$ J(a*b)[n] = \sum_{p=0}^{N-1} a[p]b[-n+p].$$
 On the other hand we have

$$((Ja)*(Jb))[n] = \sum_{p=0}^{N-1}a[-p]b[n+p].$$
so we first flip $a$ and $b$ and then apply the convolution. However, my book says $J(a*b)[n] = ((Ja)*(Jb))[n]$, which is not the case with my calculation, so I'm sure did something wrong. Can you help me finding the logical mistake here?

Your first version of $J(a*b)[n]$ is correct:

$$J(a*b)[n] = \sum_{p=0}^{N-1} a[p]b[-n-p]$$

because you simply have to flip $n$, and not the summation index.

The problem is in your formulation of the circular convolution of $(Ja)[n]$ and $(Jb)[n]$. It should be

$$((Ja)*(Jb))[n]=\sum_{p=0}^{N-1} a[-p]b[-n+p]$$

because you have to negate the arguments of both sequences, so $p$ becomes $-p$, and $n-p$ becomes $-n+p$. Now you can simply change the sign of $p$ because it doesn't matter in which order you compute the sum. This gives

$$((Ja)*(Jb))[n]=\sum_{p=0}^{N-1} a[p]b[-n-p]$$

which is identical to $J(a*b)[n]$ as given above, as it should be.


Normalizing loudness level of a set of songs I like to listen to instrumental music while trying to sleep. So, I've selected 20 songs from my music archive to listen to when I go to bed. The problem is that the loudness levels of the songs differ. Therefore, I have difficulty listening to some songs while the others disturb me for a specific volume level. I want to normalize the loudness of them and set it to a common value programmatically (e.g. using MATLAB).

I thought the following:

n <- number of songs
a <- an array of values initially 0 (length n)
for i = 1 to n
    partition song i into overlapping windows (e.g. window length of 1024, shift of 512)
    let s_j's involve the values in j'th window where j = 1 to #windows
    a(i) <- max_j sum(abs(s_j))
end for
m <- mean(a)
for i = 1 to n
    rescale all the amplitude values of song i by m/a(i)
end for

What should I take into consideration to normalize the songs' loudness levels? What strategy do you recommend? I can write the MATLAB code easily from a pseudocode.

You can spare yourself time by using off-the-shelf software such as http://mp3gain.sourceforge.net/


Oversampling in OFDM Why do we need an oversampling in OFDM signal? and how to do it?
If around half of the Sub-carriers are only used to send data, then do we need to perform oversampling to the signal? Is it necessary to have LPF in OFDM transceiver?

Oversampling in OFDM is usually implemented by modulating some subcarriers at the spectrum margins with zero. This zero-padding in frequency domain corresponds to oversampling in time domain. If only 50% percent of the subcarriers are used, then this corresponds to an oversampling factor of 2. Specifically, if the total number of subcarriers is $N$ and if $M$ of these subcarriers are set to zero, then the oversampling factor is $N/(N-M)$. Note, that the $M$ zero subcarriers have to be arranged in a row, cyclic shifting allowed. Normally, they're arranged symmetrically at the lowest and highest frequencies.

An OFDM transmitter might require a low pass filter in order to keep out-of-band radiation low. The cut-off frequency is usually half of the sampling frequency ($f_\mathrm{s}$).

An OFDM receiver requires a low pass filter before A/D conversion (sampling) with cut-off frequency $f_\mathrm{s}/2$. Otherwise aliasing can occur.


Properties of DTFT to Infer the Inverse DTFT of Altered Data Given a signal $ x \left[ n \right] $ and its DTFT $ X \left( {e}^{j \omega} \right) $.

Which property of the DTFT allows you to easily compute the inverse DTFT of $ \frac{4 X \left( {e}^{j \omega} \right)}{\pi} - 2 $?

The property is the Linearity of the DTFT.

Linearity means that if your input is a linear combination of signals the output will be the same linear combination of each input by itself:

$$ \operatorname{DTFT} \left( \alpha x \left[ n \right] + \beta y \left[ n \right] \right) = \alpha \operatorname{DTFT} \left( x \left[ n \right] \right) + \beta \operatorname{DTFT} \left( y \left[ n \right] \right) $$

This also holds for the Inverse DTFT.

By the way, if you want people to keep answer your questions:

Mark them as solved once someone solves them for you.
Don't use images to post the questions. Write them so the question will be independent of out of site resources.


Ideal Low-Pass Filter Output? I'm very new to digital signal processing, I've been reading a text book to try to figure out this assignment at work (avionics related).  I'm supposed to implement a two-pole low-pass filter with a fixed cutoff frequency (of 0.02Hz)

I'm not sure if my current output is correct, given my understanding of how low pass filters are supposed to work.

The top is my input signal. 3Hz to the left, 1Hz to the right.

Given my current level of understanding of digital signal processing, my expected output of a low-pass filter with a frequency cutoff of 2Hz is in the middle. 0 to the left, and the original 1Hz signal on the right.

The bottom is the output of my current implementation. I don't know if my current implementation is wrong, my understanding, or my expectations of how close to ideal the filters actually are.

This actual implementation is a single pole recursive filter, just trying to get a working example (and understanding) I can test against, before I design everything the way they want it.

The sample rate is 100Hz, so for f = 2Hz, I tried poles of both x = e^(-2*pi*2/100) and x = e^(-2*pi*2/200) (the formula I saw wanted a value from 0 to .5, so assumed it was some fraction of sample rate?)

Any ideas?

From the description I'll take you have a filter with 3 coefficients. (and the formulae are missing the symbol i). With a sample rate of 100 Hz you would require a filter of length 33 to make it cover the complete wave form of one cycle of 3 Hz signal. Just looking one tenth of a waveform a filter would have to do magic to attenuate (just) a 3Hz signal but leave 1Hz signals intact.

If you'll have octave or matlab, you might want to see the frequency response of your filter:

 `a=exp(i*pi*2/100);freqz([1,a*conj(a)],1)`

Here the 3dB cutoff is approximately at Fs/4 (one quarter of the sampling frequenzy).


Variable sample rate interpolation I have a set of data over which I would like to interpolate, with a sampling rate say about $f_0$ Hz with a significant uniform random sampling jitter, such that it's more or less:

$$f_0 - \frac{1}{4}f_0 \le f_s \le f_0 + \frac{1}{4}f_0.$$

Normally, with uniform sampling, I can perform sinc interpolation:

$$x(t) = \sum\limits_{n=-\infty}^\infty x[n] \cdot \mathrm{sinc}\left(\frac{t-nT}{T}\right), \quad \text{with}\quad T = \frac{1}{f_s}$$

Can I do something similar if my sampling rate is non-uniform?

Apart from the fact that $\rm{sinc}(x)=\sin(\pi x)/(\pi x)$ is smooth and bandlimited, the important property for interpolation is that $\rm{sinc}(k)=\delta(k)$ for $k\in\mathbb{Z}$ and the Kronecker symbol $\delta$.

For an irregular sampling grid you want a function with similar properties, i.e. smooth, (approximately) bandlimited and vanishing on the grid points other than some $k_0$, where it should evaluate to 1.

Fortunately, the sinc function is an analytic function, and for functions of that kind knowing all the zeros already defines the function up to a global factor. Let's look how that works for the sinc function. We know the zeros are at $\mathbb{Z}\setminus\{0\}$ (and fortunately there are no non-real zeros) and the scale factor is determined by $\rm{sinc}(0)=1$.

With this we can rewrite the function as a collection of infinitely many linear factors:

$$\rm{sinc}(x)=\frac{\prod_{k\in\mathbb{Z}\setminus\{0\}}(x-k)}{\prod_{k\in\mathbb{Z}\setminus\{0\}}(0-k)}=\prod_{k\in\mathbb{Z}\setminus\{0\}}\frac{(k-x)}{k}$$

This function keeps behaving well if you replace the zeros that are at the integers $k$ with slightly individually shifted zeros at $x_k$, which are in general non-integer. Since the grid is non-regular, it makes sense to define a distorted sinc-like function with a center index $k_0$ in mind, which corresponds to a position $x_{k_0}$, at which the function evaluates to 1.

So the new generalized interpolation function is

$$\rm{sinc}_{\{x_k\}}(x,k_0):=\prod_{k\in\mathbb{Z}\setminus\{k_0\}}\frac{x_k-x}{x_k-x_{k_0}}$$

for a given set $\{x_k\}$ of sample points and the center sample $x_{k_0}$. The resulting interpolation function is the sinc function for a regular grid and a slightly deformed smooth function for a non-regular grid. The bandwidth of the function is a continuous function of the grid points, that means small deviations from the regular grid will preserve band-limitedness quite well.

To evaluate this function you have to multiply infinitely many fractions, so this is not suitable for numerical evaluation. The function does also not converge very quickly, so that in practice the number of terms is still rather large. This changes if you apply a window function first so, that the interpolation is practically restricted to a small neighbourhood of the center sample. Then it suffices to evaluate a rather small number of terms to get a very good approximation.


Calculate cut-off frequency of lowpass IIR-filter Given a an IIR lowpass-filter in z-space:

$ H(z) = \frac{\sum_{i=0}^P b_{i} z^{-i}}{1+\sum_{j=1}^Q a_{j} z^{-j}} $

How to calculate it's 3dB cut-off frequency?

I about evaluating it's fourier transform (as long it exits) and take it's magnitude:

$\frac{1}{\sqrt{2}} = \frac{|H(e^{j\omega_c)}|}{|H(e^{j0})|}$

This is then transformed to calculate $\omega_c$. Is this correct? Are there other ways to do that?

This will work but only if the it's actually a low pass filter. Your formula describes any possible IIR filter so it could really be anything.


What is the difference of windowing functions for FIR filtering? What is the difference of window types such as Hann, Rectangular, Hamming in the application. Is there such a thing that the one is better than the other?

There will be an oral exam for DSP course so I want to understand this concept. I know that they have different windowing functions and also different time&frequency domain responses. Moreover, their parameters (such as main lobe width, side lobe width) are different. But I can't say which one is better for which application or I can't compare them intuitively.

So, from the discussion in the comments it's clear you know most you need to know.

The window method for FIR filter design is based on this idea:

We know the "ideal" frequency response $H(f)$  we want. Often, that's something like a rectangle in frequency design.
Well, the easiest thing to achieve that shape would simply be transforming $H$ to time domain, yielding the impulse response $h(t) = \mathcal F^{-1}\{H\}(t)$.
small problem here: If $H(f)$ is very steep, then the properties of the Fourier transform implies that $h(t)$ will be very long – in fact, if $H(f)$ is a rectangle (we call that the "brickwall" or "boxcar" filter, by the way), then $h(t)$ will be a sinc; and that is infinitely long. That's not useful for a FIR....
So we decide we need to cut off $h(t)$ after some length. Now, if we simply decide to cut off at some point, we have multiplied with a rectangular window $w_\text{rectangle}$, and we know that this is equivalent with convolving with a sinc in frequency domain. And convolving our nice ideal frequency response $H(f)$ with a sinc will totally "smear" and make it round and wobbly – not really what we want, usually.
Therefore, instead of using the rectangular window, we use a "rounder" window.
Thus, we realize we need to find a trade-off between being "sharp" in frequency and not "too smoothly running out" in time domain (because that means our filter remains long).

Now, the filters you mention have different properties, and were all designed to fulfill some roles.

For example, the Hamming window has a very nicely suppressed first sidelobe in frequency domain – which means that after convolving our ideal $H$ with it, the transition from passband to stop band will show the same nice behaviour, and signals close to the edge of the passband are already well suppressed. That might be useful when you want to extract a single channel from many channels that lie close together in spectrum.

Hann, on the other hand, doesn't suppress the first sidelobe as much - but in exchange, the power of the sidelobes drop faster the further you get from the mainlobe. This helps you if you want to design e.g. an anti-aliasing filter, where all the frequency that you don't suppress ends up aliased in the passband.

Now, for your oral exam: Maybe head over to wikipedia's list of window functions, and look at the frequency domain plots of them. When your examiner asks you what the up- and downsides for a given window for a given application are, you'd start by explaining that you need to convolve the ideal frequency response with the spectrum of the window, and then work your way towards an answer.

If your examiner is like any of the examiners that I was the transcript writer for (or any examiner that orally examined me when I was still a student), then that ability to work towards a solution based on an understanding of windowing and an analysis of what the  application calls for is far more interesting than the ability to have an answer to the question "What window will you use for designing an anti-aliasing filter?" instantly ready.


What is the best codec for streaming of audio in this case? I am trying to stream audio from device. I can't use the PCM audio because it is not compact enough, and I don't have so much bandwidth.

What is the option to substitute wave PCM with some codec? I am interested in good quality for speech and sampling frequency up to 16 kHz. Is ADPCM a good choice? Where can I find info about that?

Opus would be a good choice. It claims to be the superior solution for Internet (i.e. lossy transmission with inconsistent latency) music and speech streaming applications. Specifically, it provides tight control over encoding latency when that's important. It's based on technology from two previous codecs: SILK, developed by Skype for speech transmission, and CELT, a codec developed for low-latency streaming audio.

One nice aspect of Opus is that it is an open standard, so a reference implementation is available that has source code for both floating-point and fixed-point computations. It also claims to have support for multiple architectures, including x86, ARM, and MIPS, which might be useful.

As far as quality goes, there are a number of detailed comparisons available between Opus and other audio codecs like AAC or MP3. Here is a plot that seems to summarize it well:

Basically, the makers of Opus claim superior audio quality to its competitors over a wide range of bitrates. Especially in the pretty-low-bitrate range of 16-32 kbps, it seems to have no peer at the present time. This page has a number of live tests that you can run in your browser, showing its performance with both speech and audio, its robustness to packet loss, and its ability to change bitrates (and therefore audio quality) on the fly.


Understanding Voss-McCartney pink noise generation algorithm I'm implementing the Voss-McCartney pink noise generation algorithm.

If you follow the link above, you can read:

  from James McCartney  2 Sep 1999 21:00:30 -0600:

    The top end of the spectrum wasn't as good.  The cascade of sin(x)/x
    shapes that I predicted in my other post was quite obvious.
    Ripple was only about 2dB up to Fs/8 and 4dB up to Fs/5.  The response
    was about 5dB down at Fs/4 (one of the sin(x)/x nulls), and there was
    a deep null at Fs/2.
    (These figures are a bit rough.  More averaging would have helped.)

  You can improve the top octave somewhat by adding a white noise
  generator at the same amplitude as the others. Which fills in the
  diagram as follows:

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
x x x x x x x x x x x x x x x x
 x   x   x   x   x   x   x   x
   x       x       x       x
       x               x
               x

  It'll still be bumpy up there, but the nulls won't be as deep.

If I understand it well, this algorithm generates pink noise by adding random (white?) noise sources at different frequencies1

However, I don't fully understand the explanation given in the quote above for the extra white noise generator on the "top row". Can someone clarify how/why it improves the algorithm? Does that make it a good algorithm for pink noise generation for audio applications? Especially, shouldn't I discard the first samples until all the "rows" were mixed into the signal (in the ASCII art quoted above, that would mean discarding 15 first samples)?

1 I'm not sure of the wording here. Do not hesitate to correct me if I'm wrong


Matlab filter implementation I have a filter response in the $z$-domain.
$$H_0(z) = 2 + 6z^{−1} + z^{−2} + 5z^{−3} + z^{−5}$$
How do I implement this filter on MATLAB?
I can do this if I have the coefficients but I do not know how to get them.

Your filter is a 5-tap delay line FIR filter. The FIR filter's coefficients are [2 6 1 5 1]. If x is your filter's input sequence and you execute the following MATLAB command:
h = [2 6 1 5 1];

you can then implement the filter (to produce a y output sequence) using the following MATLAB command:
y = filter(h, 1, x);

You can see your filter's frequency magnitude and phase responses by executing the following MATLAB command:
freqz(h, 1, 256)


Impulse Response of real coefficient, LTI System I'm trying to obtain the impulse response $h[n]$ of a system whose frequency response is $H(e^{j\omega})=R(\omega)e^{-25j\omega}$. I believed that $h[n]=h[n-25]$, would be the correct answer, however I was told it would be something along the lines of $h[n]= h[50-n], 0<=n<=25$.
Can somebody explain why?

It must be added to the problem that $R(\omega)$ is a real-valued, possibly bipolar function. In that case, its inverse discrete-time Fourier transform must be even:
$$r[n]=r[-n]\tag{1}$$
From the given relation between $H(e^{j\omega})$ and $R(\omega)$ it is clear that
$$h[n]=r[n-25]\tag{2}$$
must hold. I'm sure that you'll manage to combine $(1)$ and $(2)$ to come up with the desired result.


LPC produces unexpected residual I have written a vocal pitch shifting algorithm which uses LPC to separate the signal into its spectral envelope (predicted by the coefficients) and the pitch information (the residual error).

My understanding (and intuition) was that the residual signal would be unintelligble, but if I listen to the residual created by algorithm it is easily identifiable as a very low volume, slightly distorted version of the original signal.

Is this correct, or should I be getting what I originally expected?

Cheers,
Ben

It sounds like you're doing things correctly.

Think of LPC as a compression method. A lossy method at that, if you discarded the residual and replaced it with a synthetic pulse train (which might consider using for pitch shifting).

In a perfect world where the LPC coefficients were a perfect match to the formant filtering applied by the throat/mouth/nose, the residual would be exactly the glottal impulse.

However, since the LPC is just an approximation, the residual contains additional information.


Are there amplitude and phase characteristics of harmonic signal components? Nonlinear amplifiers are characterized by their AM/AM and AM/PM. That's what I thought. But now I have a question: The amplitude and phase distortion characteristics, do they only belong to the fundamental signal component?

I mean, to fully characterize my amplifier, are there AM/AM and AM/PM characteristics for the harmonic components as well?

If the amplifier is memoryless (i.e. there is only AM/AM conversion and no AM/PM conversion) and the input signal is bandpass (centered around some carrier frequency $F$), for input signal:

$$
x(t) = a(t)\cos(2\pi Ft + \alpha(t))
$$
the output signal
$$
y(t) = v(x(t)) = v(a\cos(2\pi F t)) = v(a \cos \theta)
$$

can be expanded using a Fourier series:
$$
v(a\cos\theta) = v_0(a) + v_1(a)\cos(\theta) + v_2(a)\cos(2\theta) \ldots
$$

If the amplifier is memoryless, the coefficients $v_m(a)$ are real and determine the harmonic amplitudes, and there is only AM/AM conversion and no phase shift on the harmonics.

The coefficients are known as the "Chebyshev Transform" of the nonlinearity $v(x)$ and are given by:
$$
v_m(a) = \frac{2}{\pi} \int_0^\pi v(a \cos\theta)\cos(m\theta)d\theta
$$
If the amplifier has some memory, but "not much", it is known as a "quasi memoryless" system, and the coefficients are complex (giving AM/PM conversion).

To characterize the amplifier, you would need to measure $v_m(a)$ for each harmonic $m$. However, since usually all you need is the first harmonic because others will get filtered out later, you can measure $v_1(a)$ from standard power out/power in curves. (The same ones used to measure 1dB compression and 3rd order intercept). The power gain measured in this curve is just $(v_1(a)/a)^2$.

Source:
appendix C of:
https://smartech.gatech.edu/bitstream/handle/1853/5327/ku_hyunchul_200312_phd.pdf?sequence=1&isAllowed=y


How to cover all spectrum in wavelet when it which is behalved? As I know , wavelet decomposition behave with halve.
I mean that a way like this HH HL LH LL.

But I wonder, how to corver the all spectrum when it was halved?
I think if it want to all cover spectrum then it might be needed infinity halve.

How can cover all the spectrum?

It is covered - its just the signal that is the residual part of what you didn't filter yet.


Are the RLS filter and Kalman filter gradient methods? I would like to extend my previous question
What is difference between LMS and gradient-descent adaptation?
with this other question.

I found out, that RLS and Kalman filter learning seems to be somehow similar. My question is: Can be those algorithms called gradient descent methods? If not, how is this kind of algorithms called?

Recursive least squares (RLS) filters don't use gradient descent. As their name suggests, they use a least-squares fit to determine the optimum coefficients at each time step. Via clever formulation of the filter structure, one can use the calculations done from time step $n$ to recursively calculate the updated coefficients for time step $n+1$ without having to do the full least-squares fit again.


Decomposing a DFT into multiple FFT calls I'm using a good fast FFT implementation (vDSP) that will only work on power of 2 blocks of audio data.  Now I have a problem where I would like to be able to apply the calculations to non powers of 2 audio.

Some degree of reading around suggests that I can do this by breaking the DFT into smaller power of 2 blocks.

Can anyone explain how to do this, if indeed it is possible at all? I have been unable to find any explanations of how to do it ...

The simplest thing to do if you have access to only a power-of-two FFT is to zero-pad your data going into the $N=2^p$ FFT and interpolate what comes out.  i think that's what MATLAB does for fft() for non-power-of-two DFT.


Is the Canny algorithm enough for creating a feature descriptor image for input to an SVM? I retrieve contours from images by using the Canny algorithm.

Is this enough preprocessing to have a descriptor image for input into an SVM to find similarities?

Do I need other features like elongation, perimeter, area?

I am asking about this because I was inspired by this example. I have given my image in greyscale first, in the Canny algorithm style second, and in both cases my confusion matrix had plenty of 0s like precision, recall, f1-score, support measure, etc.

Yes and no.

Yes, the binary image output from the Canny edge detector is a feature set.

As the Feature (computer vision) wiki page describes:

  ... an edge can be represented as a boolean variable in each image point that describes whether an edge is present at that point. Alternatively, we can instead use a representation which provides a certainty measure instead of a boolean statement of the edge's existence and combine this with information about the orientation of the edge.

What you have is the former; a boolean variable for every pixel in the image, true where there is an edge. A Hough Transform would be an example of a method to get the latter; edges represented by measure and orientation.

If you take the same image and convert it to edges, the feature set will be the same. If you take two different images and convert them to edges, the feature set will likely be different. So yes, you can model an image using the Canny edge detector alone.

Now, no, the binary image output from the Canny edge detector is probably not a very useful feature set.

Consider a single training image of an apple and three test images: an apple, an orange, and a banana. By shape, what we are focusing on by focusing on the Canny edge detector, we would expect the test apple to be most closely related to the training apple, the orange second, and the banana last.

In a candy land of image analysis, where the images are all perfectly composed, their edges could line up such that the apples edges overlap exactly, the orange's edge overlaps a lot with the apple, and the banana's edges only overlap for the few pixels they intersect. However, in reality you are probably comparing images of different apples with different backgrounds and different composition. Imagine the test apple is shifted one pixel down and to the right, now it hardly overlaps with the training apple at all; the orange will be a better match. In this example, the binary image output from the Canny edge detector is not a robust feature descriptor.

Think of the Canny edge detector algorithm as a transformation of the image. It can be a useful step in getting to a powerful set of features which describe your image, but it itself is not usually a powerful feature descriptor.

So, what would be a powerful feature descriptor?

..well you didn't ask that ;)


Ratio of expected values of squares of errors in power quantities, in dB An amplitude quantity $a$ is estimated by two estimators $\hat a_1$ and $\hat a_2.$ The error of estimator $\hat a_2$ is compared to that of $\hat a_1$. One possible comparison is the ratio of expected values of absolute error:

$$\frac{\operatorname{E}[|a - \hat a_2|]}{\operatorname{E}[|a - \hat a_1|]}.$$

This is a ratio of amplitude quantitities and can be expressed in decibels as:

$$20 \log_{10}\left(\frac{\operatorname{E}[|a - \hat a_2|]}{\operatorname{E}[|a - \hat a_1|]}\right) \text{ dB.}$$

Another possible comparison is the ratio of expected values of square error:

$$\frac{\operatorname{E}[(a - \hat a_2)^2]}{\operatorname{E}[(a - \hat a_1)^2]}.$$

This is a ratio of power quantities and can be expressed in decibels as:

$$10 \log_{10}\left(\frac{\operatorname{E}[(a - \hat a_2)^2]}{\operatorname{E}[(a - \hat a_1)^2]}\right) \text{ dB.}$$

The squares of the estimators can be used to estimate $a^2.$ The estimators $\hat a_1^2$ and $\hat a_2^2$ can be compared by the ratio of expected values of square error:

$$\frac{\operatorname{E}[(a^2 - \hat a_2^2)^2]}{\operatorname{E}[(a^2 - \hat a_1^2)^2]}.$$

This is a ratio of squares of power quantities. How should it be expressed in decibels?

$$10 \log_{10}\left(\frac{\operatorname{E}[(a^2 - \hat a_2^2)^2]}{\operatorname{E}[(a^2 - \hat a_1^2)^2]}\right) \text{ dB,}$$

or perhaps:

$$5 \log_{10}\left(\frac{\operatorname{E}[(a^2 - \hat a_2^2)^2]}{\operatorname{E}[(a^2 - \hat a_1^2)^2]}\right) \text{ dB} = 10 \log_{10}\left(\frac{\sqrt{\operatorname{E}[(a^2 - \hat a_2^2)^2]}}{\sqrt{\operatorname{E}[(a^2 - \hat a_1^2)^2]}}\right) \text{ dB},$$

where the square roots are power quantities, or is there no commonly agreed convention?

Your question raises many very present concerns about measuring differences, and optimizing, in signal/image processing. Measuring can help:

to compare different outcomes from two different processes,
to optimize processing per se (for denoising, restoration).

Most use cost functions related to  $p$-norms (when $p\ge 1$) or quasi-norms (when $0<p\le 1$): $$\mathcal{C}_p(a,\hat{a})=\left(\sum_k|a_k-\hat{a}_k|^p\right)^{1/p}\,.$$ Those functions are symmetric, and are consistent with the physical units of $a$. The behavior of $\mathcal{C}_p$ depends on power $p$: classical norm inequalities tells you that $\mathcal{C}_p(a,\hat{a})\ge \mathcal{C}_q(a,\hat{a})$ when $p\le q$. Those (quasi-)norms are $1$-homogeneous. Using different $p$ or $q$ can be useful, especially in a sparsity context, see for instance references in: Euclid in a Taxicab: Sparse Blind Deconvolution with Smoothed \ell_1/\ell_2 Regularization. Here is a rendition of classical penalties:

Thus, quasi-norm or norm ratios are used to compare estimates, in a unitless fashion. The $\ell_p$ signal-to-noise ratios are typically of the shape:

$$\frac{\mathcal{C}_p(a,0)}{\mathcal{C}_p(a,\hat{a})}\,.$$

They are unit-less (ratios of same units). To better express orders of magnitude, one usually take logarithms. The logarithms can also been interpreted as unit-less measure. Interestingly, either $\log$ or $\ell_p/\ell_q$ ratios (like $\ell_1/\ell_2$) have been used a non-convex penalties in signal restoration.

Your functions with power $\alpha$ can be written as:

$$B\log_b \frac{\mathcal{C}_p(a^\alpha,\hat{a}_2^\alpha)}{\mathcal{C}_p(a^\alpha,\hat{a}_1^\alpha)}\,.$$

Now, those logarithms of ratios are unit-less, and are usually employed in a relative way, for comparison only. With $B>0$, $b>1$, $t \to B \log_b (t)$ is a monotonic increasing function. So, for value comparison only, we don't really care about the values for $b$ and $B$. However, for consistency in dB scale (as Olli Niemitalo already answered properly), signal-to-noise ratios tend to be computed in $10\log_{10}$ of energy measures, so homogeneously:

$ 20\log_{10}(|a|)$
$ 10\log_{10}(|a|^2)$
$ 5\log_{10}(|a|^4)$

make sense. However, it is not easy to compare them objectively, but they can be plotted on similar graph.


Study Signal Processing I'd Like to ask two questions :

What is the difference between studying Signal processing (both Deterministic and statistical) in Department of Electrical Engineering versus Department of Mathematics and statistics

What's All the mathematical topics (both time/n-domain and f/S f/Z domain ) which is a prerequisites of studying Signal Processing (Both continuous and discrete) or reading a book like "signals and systems" by "Alan Oppenheim" ??

Thanks All

i remember this question being asked on USENET at comp.dsp, but i don't remember all of my answer.

there are several mathematic disciplines that inform signal processing and formalize some of the concepts.

let's assume you have calculus of one and many variables and differential equations (any scientist or engineer should have those).

well, Complex Variables and Analysis of Complex functions to start with.

Linear Algebra: Matrices and Determinants

Fourier Series and Transforms

Probability, Random Variables, and Random Processes

Numerical Techniques

Metric Spaces and Functional Analysis

all those are useful math topics, sufficiently useful that i would recommend taking math classes in these if you end up in grad school in EE.  there are many signal processing topics that make reference to one or more of the above math topics.


How do I recover the signal from an ECG image In my project, I have to digitize an ECG image taken with a normal camera (jpeg). For example, I have the following camera captured image:

and I want to get something like this:-

and then the digitized data (x,y points), like in this video on digitization of ECG

I have no idea how to do it, so I searched and consulted several research papers. The general approach of the algorithms is like:-

change to gray level image
delete gridlines
add missing points
convert 2D image to 1D image

I am stuck with the second point, i.e., deleting the gridlines. I looked up some more references to do this and found histogram analysis might be helpful.

Can you please guide me on how to do this (I'm using MATLAB 2010)? Any help would be appreciated.

Sorry, I use Mathematica, but it should be really easy to implement the idea into Matlab. I give the code anyway, so when my description is not detailed enough, you can get the rest from the code.

Basic idea is: you look at your image column-wise. Process every column of pixels separately. Note, that in the plot, I inverted the gray-values. So black is 1 and white is 0.
If you plot the (inverted) brightness-pixel-values you have basically only two situations. The first one is, when your column is not on a vertical grid line. The the plot looks like

The second situation is, where you are directly on a vertical grid line. Then the grid line influences the brightness of the whole column

But what you see is, that your dark EEG seems to be always the maximum. Therefore, the very complex algorithm is: Go through every column and take the position of the most black pixel.

img = ColorConvert[
   ImagePad[
    Import["http://i.stack.imgur.com/500Kg.jpg"], {{0, -20}, {0, 0}}],
    "Grayscale"];

Image[
 Transpose[Function[With[{m = Min[#]},
     Map[Function[{v}, If[v == m, 1, 0]], #]]] /@
   Transpose[ImageData[img, "Real"]]
  ]
 ]

Note, that I croped a bit of the right side of you image, since it was completely white there. The result is

Now you can join the points or interpolate them in every way you like and you get your EEG


How to calculate the similarity of two electric waves I have two electric signals from different sources and I think that are similar.
I have the FFT of the signals and I have their spectrums of frequency  in separate arrays.
How should I proceed to find the percentual of similarity the two signals?

You could measure the Cross-Correlation between the two signals.

Here's an almost identical thread, with a filling answer:

How do I implement cross-correlation to prove two audio files are similar?

Some general info on cross-correlation from wiki:
https://en.wikipedia.org/wiki/Cross-correlation

Hope that helps.


Why must I take abs of filter frequency response? I was playing with filters in SciPy and noticed something that I'm not that familiar with:

I'm using the code:

import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt

def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y

# Filter requirements.
order = 6
fs = 30.0       # sample rate, Hz
cutoff = 3.667  # desired cutoff frequency of the filter, Hz

# Get the filter coefficients so we can check its frequency response.
b, a = butter_lowpass(cutoff, fs, order)

# Plot the frequency response.
w, h = freqz(b, a, worN=8000)

taken from: https://stackoverflow.com/a/25192640/4959635

And doing some plots of the given freqz results:

plt.plot(w,h)
plt.show()

plt.plot(w,abs(h))
plt.show()

So why must I take abs to get "the plot that I want"? And what's the first one with negative amplitudes?

Frequency response has two parts: amplitude response and the phase response. Both of these are represented as a complex signal when you get the response from freqz. In order to plot the amplitude response you need to use abs. Otherwise I doubt it only shows you the real part which I think is what you see in the first figure. Note that when dealing with the amplitude response, in most cases you will do a further step and convert it to $\textrm{dB}$ scale.


Connection between filter equation in frequency domain and difference equation in time domain? I know that the frequency response of a filter is described by equation:

H(w) = 1-1.176*e^(-jw) + e^(-j2w) , where w is the angular frequency.

Then in time domain I have a filter that is described by difference equation:

y(k) = x(k)-1.176*x(k-1)+x(k-2)

Now I'm asked to explain the connection between these equations and why the coefficients for the both equations are the same [1 -1.176 1]?

I've been reading books and trying to find information around the internet but without any solutions for this. Could someone explain me the connection of these equations? Why are the coefficients same?

I have an idea that these two equations might represent the same filter but in different domains (time and frequency). Does convolution have something to do with this?

Hint:

The transfer function $H(\omega)$ is equal to the ratio of output and input spectra:

$$H(\omega)=\frac{Y(\omega)}{X(\omega)}$$

So your first equation can be rewritten as

$$Y(\omega)=X(\omega) - 1.176X(\omega)e^{-j\omega} + X(\omega)e^{-2j\omega}$$

What do you get if you transform this equation to the time-domain?


What is the general formula for radon back projection for a javascript implementation? I'm aiming to implement a simulation of computed tomography back projection in javascript/HTML5 canvas.

Trying to figure the correct approach for doing a back projections and I have been studying the presentation below:

http://www.dtic.upf.edu/~afrangi/ibi/reconstruction_color_2.pdf. On pp. 9-13 of the above pdf presentation a step by step example of back projection is shown.

This example involves a 2x2 matrix (with 4 projections and a total of 10 line integrals) and computes the line integrals of each projection where the value of each corresponding line integral is added to each matrix cell, respectively. Subsequently the total number of line integrals is subtracted from the value of each matrix cell. Finally the resulting value in each matrix cell is divided by the number of projections minus one.

Does this approach hold in general for larger, say a 100x100 matrix with larger cell values, say around 100?

To implement projection the simplest way is to rotate your image then sum over a row or column. The simplest way to implement back projection is to take a line of your sinogram which is a projection from certain view (angle) then repeat that line to form an image then rotate that image regarding the projection view, then sum all of backprojected views.
There is a unique formula for continuous case but for discrete case there are many different implementation, this one is the simplest.

Good luck.


Is Spectral Coherence (MSC) dependent to the initial amplitude of the comparing signals? Given the signals

\begin{align}
x(t) &= A_1\sin(\omega t)\\
y(t) &= A_2\sin(\omega t)
\end{align}

Do $A_1$ and $A_2$ play a role in the $C_{xy}$ coherence estimator? MATLAB's mscohere.m shows no differences for $A_1=A_2$ and $A_1\neq A_2$ at all.

From the definition of the (magnitude-squared) coherence

$$C_{xy}(f)=\frac{|G_{xy}(f)|^2}{G_{xx}(f)G_{yy}(f)}$$

with the cross-spectral density $G_{xy}(f)$, and the power spectra $G_{xx}(f)$ and $G_{yy}(f)$, respectively, it is clear that scaling of $x(t)$ or $y(t)$ does not change the value of $C_{xy}(f)$, because the scaling constants appear in the numerator as well as in the denominator and cancel out.

The coherence is a normalized quantity (and hence scaling-invariant), and it always satisfies $0\le C_{xy}(f)\le 1$.


Could 2 synchronised RX chains be used to receive a signal too wide for just one of them? Say the goal is to demodulate a "wide" signal, e.g. 300 MHz wide. An SDR is available with a maximum RX bandwidth of 150 MHz (ADC sampling freq. limitation), but it comes with 2 RX channels sharing the same clock source.

The proposed approach is connecting the same antenna to both RX chains. Ideally, they would take the same clock source, but one ADC would trigger on a rising edge and the other on the falling edge of the clock, thus effectively doubling the ADC sampling rate and allowing downstream demodulation logic to correctly demodulate a 300 MHz signal.

If this is practically doable, what are the constraints (w.r.t. SDR analogue RF frontend capabilities, compensation for loss of signal strength due to splitting it into 2 RX inputs etc.) and if not, why not?

In lack of more background, I can only give a general answer:

What you describe is simply an interleaved ADC; that is a technique employed in very high-rate ADCs. And high-rate definitely is more in the range of hundreds of MHz, not a mere handful! For your proposed rates, cheap ADCs can directly sample sufficiently fast. In fact, in the context of SDR, you'll find that most ADCs used are actually doing higher rates on two simultaneously sampled channels in order to allow IQ signals to be meaningfully digitized.

So, yes, that is possible, but it's incredibly hard to do well-synchronized and balanced enough. Just get a single sufficiently fast ADC. In fact, I'm not even sure what SDR-centric ADCs you're looking at that these rates become a concern at all.


How to proceed for BOW creation From sift features I extracted sift feature to construct a bag-of-words model. I did the following to extract the feature vector:

Extracted features from the key frames.
After the feature extraction I tried using Principle Component Analysis to reduce the size of my features.

Query: What further pre-processing steps should I take before creating a codebook out of the set of features extracted ?

There is now support for the bag-of-words model in the Computer Vision System Toolbox for MATLAB.


How do I compare 2 audio files, returning a percentage of difference Is there a way to compare 2 audio files, whereby a percentage of difference will be returned?
I am using C# VS08, .net framework 3.5.

Edited: difference in the sound (like for eg. Audio 1: "HELP", Audio 2:"HELP ME PLEASE", it will return the difference in percentage between this 2 audio.)

I am currently saving recorded files in wav format, and its bit rate and quality will also be the same as what will be recorded and compared against. I wish to ask if there is a way to compare the sound (not just the binaries) in 2 audio files of the same properties and same format (WAV), and show the percentage difference between the 2 sounds in the 2 audio files.

Help will be greatly appreciated.

Could you please define the "percentage of difference"?

For example, please tell us which value do you expect to get:

Between two versions of the same file encoded with different settings?
Between two versions of the same file, except that one is slightly faster than the other?
Between two recordings of the same audio input, but through different microphone/recording equipment?
Between two recordings of the same person saying the same word?
Between two recordings of the same person saying the same word, with a different prosody (rhythm/melody)?
Between two recordings of different persons saying the same word?
Between two recordings of the same person saying the same sentence except for a few words?
Between a recording of a person and that of a cow?
Between a recording of an airplane and music?

What do you want to measure? The speech contents (words)? The melody, rhythm? The overall audio similarity? Encoding/recording equipment differences?


How can I lowpass filter by only reducing peak data? I have a 2D image, which I want to lowpass filter, with these constraints / quality metrics:

I can not "add" light to the image, so each pixel in the result
should be <= the corresponding pixel in the input.
The lowpass cutoff frequency should be a parameter, to experiment with
Applying this filter repeatedly should not change the result in a significant way.
The time it takes to run this algorithm (5 minutes for a 5MPix image seems reasonable)
Minimizing the amount of light that is filtered away.

Below are some approaches I have tried, together with their shortcomings:

Gaussian filter like normal, then pull the result down to comply with constraint 1. This complies very well with the first 3 points, but reduces a lot more light than necessary.
Fitting "upwards" parabolas through the "low" points and "downwards" parabolas between them to smooth out. This works great in 1D, but applying it first horizontally, then vertically produces bad results in 2D. It takes a lot longer, but not too long for my application. However, repeatably applying this filter will drastically change the result. If the (1D) input is a perfect "downward" parabola (which should not be filtered at all), it will be replaced by 2 "upward" parabolas sitting at the start/end.
Using some other form of 2D "basic" functions and linear solving to find the optimal parameters. This is an idea only currently, not implemented / tested yet.

My domain of experience in signal processing is almost exclusively image processing, so I hope to find alternatives to this problem with the input of experts active in other areas of signal processing.

update 2011/08/18

Based on the current reactions, I decided to make things a bit more clear by adding graphs of a typical input and the results of the 3 approaches I described originally + the suggestions I received so far.
For easy comparison, I used just 1D filtering in these examples.

Input data:

Gaussian filter + bring it down to comply with requirement (1).
You can see that bringing it down results in unnecessary light reductions on the right hand side.

Parabolas
As far as I am concerned, this is pretty much excellent, sadly it does not translate perfectly into 2D by applying first horizontal, then vertical. In this case, you also see I can evaluate the fitted parabolas in floating point resolution, which is a small benefit, but not absolutely required.

Grayscale erosion
Based on the suggestion from rwong, I tried grayscale erosion. I used a structuring element with the same parabolic shape as my "fitted" parabolas. The result is almost exactly the same, so this looks promising. However, there are still a few problems:
 1. My structuring element was not "big enough" (although it was already 801 pixels wide)
 1. I only have "upwards" parabolas, no "downward parabolas to smooth the transition from one parabola to the next.

Median filtering
Only included for completeness, it is not really what I want.

raw data
I pasted the raw input data + the various python commands onto pastebin, so you can experiment with the same data too.
http://pastebin.com/ASnJ9M0p

There is indeed a 2D version for your attempt #2 - it is similar in theory, but it cannot be decomposed into two 1D operations. Please read about "2D grayscale morphological filtering". It is faster than curve fitting.

http://opencv.willowgarage.com/documentation/image_filtering.html
http://en.wikipedia.org/wiki/Dilation_(morphology)

Median filtering might also be useful if you are trying to remove speckles. A more advanced form of median filtering is "ordinal filtering".

In all cases, requirement #1 can be met trivially by taking the pixelwise minimum between the output and the input. It is an important quality criteria, but it will not limit the choice of algorithms.

Gaussian filtering (and a number of other useful filters) can be decomposed (first from 2D to 1D operations, then via Fourier transform), but there are many other useful image processing techniques that are not decomposable, which make them slow but does not diminish their usefulness.


How to Combine Low and High Frequencies of Two Images in MATLAB? I would like to combine two images A and B in the following way:

1) I want to take a Fourier transform of both of them

2) For image A I want to apply a weighted filter, which gives more emphasis for low frequencies

3) For image B I want to apply a weighted filter, which gives more emphasis for high frequencies

4) I want to combine these frequencies and take the inverse Fourier transform

Can someone give me any guidelines where I should start (which functions etc.) to do this in Matlab? =) I'm kinda learning about Fourier transform and I want to play around with images. I was hoping if someone could give an example of code how this could be done etc.

Please note that I'm new at this stuff and I'm not very familiar with all the terminology yet. I'm an amateur trying to learn about Fourier transform by doing an experiment with it :)

Thank you for any help! =)

P.S. I would appreciate if someone could give me a minimal code snippet showing me what I need to do =)

The easiest way doing so would be Laplacian Pyramid.

Yet, it can be done just by using simple Addition operator.
Just add the High Frequency of one image to the Low Frequency of the other.
Keep in mind few things:

Dimensions must be the same.
Otherwise, interpolate to the same dimensions.
It is better to use HPF which is built from the same LPF used.
To the least, they must have the same Cut Off Frequency.
Use floating point arithmetic for this procedure.