Error in using Kalman Filter for 2D Position Estimation in Python This is my first question on DSP Stack exchange, so I apologise if it is poorly worded. I have some positioning data from a vehicle (GPX Format, collected through Strava) and want to use a Kalman filter as a first step to track the position of the car and determine the accuracy of the GPS measurement. However, the output estimate of the Kalman filter from my implementation seems completely wrong and absolutely does not match the position information and plot from the data. Can someone help me figure out what is wrong with this implementation and provide some methods to fix this problem?
Appreciate all the help in advance. Thanks!

As a reference, I used the code given in the below link and correspondingly modified it based on my requirements: https://stackoverflow.com/questions/13901997/kalman-2d-filter-in-pythonn.
Background: I only have a basic understanding of the working of the Kalman filter and am a new user to Python, but for this implementation, I've considered a constant velocity model with states as Position and Velocity, time step is assumed to be 1 (Considering GPS updates at 1Hz), the measurement matrix only considers the position information and the actual measurement gives the corresponding longitude and latitude values.
The test GPX file was obtained from the following link: https://github.com/stevenvandorpe/testdata/blob/master/gps_coordinates/gpx/my_run_001.gpx

My implementation in Python:

import gpxpy
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

with open('test3.gpx') as fh:
    gpx_file = gpxpy.parse(fh)
segment = gpx_file.tracks[0].segments[0]
coords = pd.DataFrame([
    {'lat': p.latitude,
     'lon': p.longitude,
     'ele': p.elevation,
     'time': p.time} for p in segment.points])
coords.head(3)
plt.plot(coords.lon[::36], coords.lat[::36], 'ro')
plt.show()

def kalman_xy(x, P, measurement, R,
              Q = np.array(np.eye(4))):

    return kalman(x, P, measurement, R, Q,
                  F=np.array([[1.0, 0.0, 1.0, 0.0],
                              [0.0, 1.0, 0.0, 1.0],
                              [0.0, 0.0, 1.0, 0.0],
                              [0.0, 0.0, 0.0, 1.0]]),
                  H=np.array([[1.0, 0.0, 0.0, 0.0],
                              [0.0, 1.0, 0.0, 0.0]]))

def kalman(x, P, measurement, R, Q, F, H):

    y = np.array(measurement).T - np.dot(H,x)
    S = H.dot(P).dot(H.T) + R  # residual convariance
    K = np.dot((P.dot(H.T)), np.linalg.pinv(S))
    x = x + K.dot(y)
    I = np.array(np.eye(F.shape[0]))  # identity matrix
    P = np.dot((I - np.dot(K,H)),P)

    # PREDICT x, P
    x = np.dot(F,x)
    P = F.dot(P).dot(F.T) + Q

    return x, P

def demo_kalman_xy():
    x = np.array([[100.0, 0.0, 0.0, 100.0]]).T
    P = np.array(np.eye(4))*1000 # initial uncertainty
    plt.plot(coords.lon[::36], coords.lat[::36], 'ro')
    result = []
    R = 0.01**2
    for meas in zip(coords.lon, coords.lat):
        x, P = kalman_xy(x, P, meas, R)
        result.append((x[:2]).tolist())
    kalman_x, kalman_y = zip(*result)
    plt.plot(kalman_x, kalman_y, 'g-')
    plt.show()

demo_kalman_xy()

The main reason why your Kalman filter is not working is because you are not converting lat and lon values to kms. In the code below, I defined a new function called lat_lon_posx_posy which converts lat and lon values to px and py values in mts. You will need to make the following changes to your code.

Include the following function

import utm

def lat_log_posx_posy(coords):

     px, py = [], []
     for i in range(len(coords.lat)):
         dx = utm.from_latlon(coords.lat[i], coords.lon[i])
         px.append(dx[0])
         py.append(dx[1])
     return px, py

Change demo_kalman_xy() to the following:

def demo_kalman_xy():

    px, py = lat_log_posx_posy(coords)
    plt.plot(px[::18], py[::18], 'ro')
    plt.show()

    x = np.array([px[0], py[0], 0.01, 0.01]).T
    P = np.array(np.eye(4))*1000 # initial uncertainty
    result = []
    R = 0.01**2
    for meas in zip(px, py):
        x, P = kalman_xy(x, P, meas, R)
        result.append((x[:2]).tolist())
    kalman_x, kalman_y = zip(*result)
    plt.plot(px[::18], py[::18], 'ro')
    plt.plot(kalman_x, kalman_y, 'g-')
    plt.show()

Final Plot:


why z-transform differentiation is needed? Find the z-transform and sketch pole-zero plot and ROC for

$$x(n)=|n|\left(\frac{1}{2}\right)^{|n|}$$

Right now, i can get the following,

$$x(n)=\begin{cases} n(\frac{1}{2})^n, & n \geq 0 \\ -n(\frac{1}{2})^{-n}, & n \lt 0 \end{cases}$$

Thus,
$$n\left(\frac{1}{2} \right)^n u(n) -n\left(\frac{1}{2} \right)^{-n}u(-n-1)$$

And,

$$\left(\frac{1}{2} \right)^n u(n) \underleftrightarrow{z} \frac{1}{1-\frac{1}{2} z^{-1}} ,|z|>\frac{1}{2}$$
$$\left(-2 \right)^n u(-n-1) \underleftrightarrow{z} \frac{1}{1-2 z^{-1}}, |z|<2$$

From here, i know my poles is $\frac{1}{2}$ and $2$. But, the solution that i have still continue working on to differentiate the equtions in z-domain and obtained the following,

$$\frac{\frac{1}{2} z^{-1}}{\left( 1-\frac{1}{2} z^{-1} \right)^{2}}+\frac{2z^{-1}}{\left( 1-2z^{-1} \right)^2} ,\frac{1}{2}<|z|<2$$

Does the differentitaion important and why does it need so?

You just wrote down the $\mathcal{Z}$-transform of $\left(\frac12\right)^nu[n]$, but you need the $\mathcal{Z}$-transform of $n\left(\frac12\right)^nu[n]$ (note the factor $n$).

In order to find that $\mathcal{Z}$-transform you can use the differentiation property (see this table):

$$\mathcal{Z}\big\{nx[n]\big\}=-z\frac{dX(z)}{dz}\tag{1}$$


OFDM DA conversion problem I'm working on a sound card based OFDM system and ran into some problems:

How can I perform  DAC on the real and imaginary parts of the IFFT? I think they are real numbers instead of 0 and 1. Thanks.

Imaginary part of a complex number is real. For example

$$z = 1 - 2j$$
$$\mathsf{Re}(z) = 1$$
$$\mathsf{Im}(z) = -2$$

The same will happen with the imaginary part of your FFT. In the end, you get two real signals. Then you modulate them by two sinusoids which are 90 degrees out of phase, and add the modulated signals at the output stage.


Sparse Channel Estimation in OFDM System Compressive sensing methods are used in channel estimation when the channel is thread as sparse where there are zero points in the signal representation.
I am in the beginning of the sparse signal processing and trying to understand how to model the sparse channel in the equation below.

Y(k)=H(k)X(k)   === > Y(k) received signal in frequency domain

My questions are:

1-What is the presentation of sparse channel h in this equation?
2-How to define the dictionary matrix for sparse channel?
3-For sparse estimation do we need to insert pilot symbols into the transmitted data x(n)?
4-Do you have a reference which explains the sparse channel properties in OFDM systems?

I would appreciate if someone please help on this topic.

Many thanks!

Are you still there, I've seen this question right now.
First let me correct your question,

Y(k)=H(k)X(k) === > Y(k) received signal in frequency domain

That's true, but it's not necessary that signal in frequency domain in order to be used in compressed sensing, it might be in time domain also.

Now, regarding your questions,
the representation of h channel can be sparse too but it's not always, it can be also in other forms .. taking an example of multi-path environment, the received signal can be sparse, and h channel can be either sparse or not, that's not important.

how to design the dictionary,  that's the goal of the compresses sensing, it's how can you design it,  you can check online, you will get many dictionary.

Estimation is different about using compressed sensing. you may use pilot or not, that's different topic.

reference for compresses sensing, is An Introduction to Compressive Sensing
Collection Editors:
Richard Baraniuk
Mark A. Davenport
Marco F. Duarte
Chinmay Hegde

Thanks


Fourier series calculation I have tried to solve, but do not know if the answer is correct or not.

A person has a periodic voltage input to a circuit. The input repeats itself every 0.02 seconds i.e. the fundamental period is 0.02 seconds.
The person measures the input with a voltmeter and finds that the input may be approximated by:

Vin(t) = V0*exp(a*t) ; 0 sec <= t < 0.02 sec ;  V0= 2 Volts a = -100sec^-1

Graphically the function may be represented by:

A)Write a mathematical expression containing an integral that may be used for calculation of the Fourier-coefficients of Vin(t).

B)Evaluate the integral analytically for the DC-term (often referred to as c0 ). Looking at the graph, how can you check if your answer is reasonable?

I know I need to use these formulas:

I've tried to calculate a0 and here I got approx 2,
according to this matlab code:

v0 = 2;
a = (-100^-1)
T = 0.02
a_0 = (1/T)*int((v0*exp(a*t)),t,0,T)

a0 calculated to 1.9998

I stuck with the an and bn, because i'm not sure if that is correct or not...

an = ((2/T)*int((v0*exp(a*t))*cos(n*w0*t),t,0,T))

an calculated to (20000*exp(-1/5000)*(exp(1/5000) - 1))/(100000000*pi^2 + 1)

bn = (2/T)*int((v0*exp(a*t))*sin(n*w0*t),t,0,T)

bn calculated to (200000000*pi*exp(-1/5000)*(exp(1/5000) - 1))/(100000000*pi^2 + 1)

Hint: this is a lot easier if you use the Euler formula https://en.wikipedia.org/wiki/Euler%27s_formula

Then the integral simply becomes
$$X(n)=\frac{1}{L}\int_{0}^{2L}e^{(-a+i\omega)\cdot t}, a_{n}=real[X(n)], b_{n} = imag[X(n)]$$

which is pretty straight forward to solve


Why plotting the fft result of a signal from -f to f gives different result compared to 0 to f? When I plot the fft result of an audio signal on the positive side of frequency axis I get the following figure, I understand the signal on the right side of the Nyquist point is just a mirror image of the original signal, the sharp noise is somewhere around 1700Hz:

Now when I plot the fft of the same signal from negative frequency to positive, that noisy component is now sitting somewhere around 2300Hz and when I check the fft of the audio signal, the sharp noise is indeed at 2300Hz and not 1700Hz, the audio content also is no longer below 1KHz and it's now instead just below 4KHz. What exactly is going on here ? Why the top plot appears to be giving me the wrong information ?

I have tried to find the reason online but I don't think I am googling the right words, hence can't find anything.

Can you advise:

What is the sample interval
What did you use to calculate the spectrum.

Edit 2:
The top one looks right, plotting from 0 Hz to 2*Nyquist. Not sure how
you calculated the spectrum, but this may be the convention
with that package.
If you prefer to present the spetrum from -ve Nyq to +ve Nyq you can't
just plot the same results with new axes.
You will need to recalculate
either with a function which uses that convention, or else do the DFT
yourself, or else you could just move the results above nyquist
and use these for the negative part.
Edit 3:
Nothing wrong, It is just an axis labelling thing.
The way your program has calculated the FFT, the answer is normaly interpreted
as 0Hz to 2*Nyquist.
Or if you prefer you could relabel it from left to right as:
(0-> nyquist), then (-Nyquist > 0)
If you want it to present as (-Nyq > Nyq) you need to rearrange the
answer. ("Roll" the right half of the array across to the left.)


Roll off rate vs filter order? The concept of "increasing roll off rate with increasing filter order" applies only to butter worth or also chebyshev and other filter types?

assuming it's a low-pass filter and assuming that the filter order equals the number of poles and there are no zeros in the s-domain transfer function, then the rolloff at the edge of the passband through the stopband is -20 dB/decade (or -6 dB/octave) multiplied by the order of the filter.
now Butterworth and Tchebychev Type 1 filters have no zeros.
Tchebychev Type 2 and Elliptical filters have zeros, so the ultimate rolloff is not there, the gain remains at the stopband level.


How to plot the discrete time signal $u[2k - 3]$? I am having a hard time plotting this graph. Do I express it as $u[2(k - 1.5)]$? But then still the step occurs for $k = 1.5$ which isn't possible for discrete time signals as they take integer values for $k$. So how to approach this?

$u[n]$ is the unit step function with the definition
$$ u[n] = \begin{cases}{ 1 ~~~, ~~~n \geq 0 \\ 0 ~~~,~~~ n < 0 } \end{cases} $$
Then you will have for $u[2k-3]$ as
$$ u[2k-3] = \begin{cases}{ 1 ~~~, ~~~ 2k-3 \geq 0 \\ 0 ~~~,~~~ 2k-3 < 0 } \end{cases} $$
$$ u[2k-3] = \begin{cases}{ 1 ~~~, ~~~ k \geq 2 \\ 0 ~~~,~~~ k< 2 } \end{cases} $$
I think this last line is easy to see. You can plot it wrt the variable $k$.


Characterizing real-time sampled signal I have a signal which is sampled in real-time. At every time interval I receive a value. I need to characterize this signal as active or inactive (to take priority in the control actions). By active it means that we observe big changes over a specific time-window.

One way to do it as I consider, is to keep a history of N samples and calculate the mean value and the standard deviation with an update algorithm. Based on the std I can say if it is inactive or not. The problem is that I have tens of thousands of signals and the sampling period is not the same on all of them. So, keeping a variable length sample history for all the signals will be difficult.

For the mean value I can use averaging with forgetting factor. The forgetting factor can be connected to the sampling period of each signal to have a similar observation window on all signals.

Now the question:

Is there any metric similar to averaging with forgetting factor that can give me the same information as std? It has to be low on memory usage and computationally light (since it's a real-time operation).

Any papers, books or references are welcome!

Thanks in advance.

Try something like this:

$\begin{eqnarray}
\mu(n) &=& (1 - \alpha_1) \mu(n) + \alpha_1 x(n) \\
\bar{x}(n) &=& x(n) - \mu(n) \\
s(n) &=& (1 - \alpha_2) s(n) + \alpha_2 \bar{x}(n)^2 \\
\sigma(n) &=& \sqrt{s(n)} \\
\end{eqnarray}$

This is equivalent to sending your input signal to a 1-pole DC blocking high-pass filter with a pole at $\alpha_1$, then sending the result to a RMS detector whose response time is set by $\alpha_2$ - to measure how much the signal is wiggling once its central trend has been removed.

Observe that $\sigma(n)$ is consistent with the definition of the standard deviation.


What is the Fourier Transform of an integral If we have a function that is an integral over the interval 0 to T

Is is correct to say that its Fourier Transform is :

When integrating functions with dirac delta you have to remember that
$$\int \delta(x - a) f(x) dx = f(a)$$ so, it works as  change of variables
and for the 2D dirac, integrating over an axis gives the 1D delta
$$ \int \delta(x - a, y) f(x) dx = \delta(y) f(a, x)$$
$$\begin{eqnarray}g(x,y)
  &=& \int_{0}^{T} f(x,y) \delta(x - \nu_x t, y) dt \\
  &=&\delta(y) f(\nu_x t, y)
\end{eqnarray}$$
This does not depend on $x$ or $u$, it is expressed in terms of $f(\nu_x t, 0)$, if you want to express in terms of $F(u, v)$ you must calculate the inverse fourier tranform for these for this point.
$$\begin{eqnarray}G(u, v)
&=& \int \left(\int \delta(y) f(\nu_x t, y) e^{-2i\pi i vy} dy \right) e^{-2i\pi i ux} dx \\
&=& \int \left( f(\nu_x t, 0) \right) e^{-2i\pi i ux} dx \\
&=& \delta(u) f(\nu_x t, 0) \\
&=& \delta(u) \int \left( \int F(u,v) e^{-2i\pi u v_x t} du \right) dv
\end{eqnarray}$$


Why do we quantize high frequencies coarsely, shouldn't it be the opposite? Also, in Huffman Coding, why do we use fewer bits for high frequency data. I don't get the intuition. Shouldn't it be the opposite or does this have to do something with the psychoacoustic model of the human hearing which says that our cochlea is sensitive to lower frequencies more than higher, hence we use up higher bits for lower frequencies?

In Huffman coding frequency refers to the number of times a particular sequence of discrete source symbols has appeared in the input. It is not the same concept as frequency in the continuous space such as sound frequency. As such, it has nothing to do with psychoacoustics or the cochlea.

We use fewer bits for the higher frequency source symbols because frequent source symbols have lower entropy and assigning fewer bits to low entropy symbols gives better compression.


What does 'canonical' mean? What does 'canonical' mean?
e.g. 'write IIR or FIR in canonical form'

I get the gist of it but maybe someone here has a precise definition. In a non-DSP context it would mean something like 'strictly by the book', and perhaps in a DSP context this would mean 'without considering implementation issues'. Is this true?

A canonical implementation of a digital filter has the minimum number of delay elements. The same filter cannot be implemented with less memory. Since memory is often not the only concern, a canonical implementation is not necessarily always the best implementation on a given platform.


How to understand the filter coefficients in FIR filter Here it is the difference equation we have:

y(n) = -0.25x(n) - 0.15x(n − 1) + 0.78x(n − 2) − 0.15x(n − 3) - 0.25x(n - 4)

I understand that is FIR. The coefficients are selected in the way that certain input signal frequencies are removed by such convolution system. Take Matlab audio file gong.mat for example, we read the file then we have input signal x and sampling rate Fs, then through such system, we have output signal y. How do we understand the filter coefficients that make some input signal frequencies are filtered? Different coefficients lead to different filter type? Can we estimate what frequencies are removed?

2-tap FIR filters can be understood with a simple intuition (I use Python for plotting, but freqz is the same there).

Let's start with two taps with the same value: b = [1, 1]. This lets through constant signals (e.g. [1, 1, 1, ...]) but will more or less block rapidly changing ones (e.g. [1, -1, 1, -1, ...]), i.e. it is a lowpass:

w, h = scipy.signal.freqz([1, 1], 1)

plt.plot(w, abs(h))

On the other hand, b = [-1, 1] will do the opposite: let changing signals pass but block constant ones, i.e., we have a highpass.

w, h = scipy.signal.freqz([-1, 1], 1)

plt.plot(w, abs(h))

I'm not sure it is straightforward to expand this intuition to more complex signals like yours. But looking at the shape, I would expect some band-pass behavior (because we have both high- and lowpass components) with a peak frequency around one third of the sampling frequency (because a signal like [1, 0, 0, -1, 0, 0, 1, ...] will be able to pass). This apparently isn't so far from the truth:


Does scaling property imply superposition? For a system to be linear,it follow the principles of scaling and superposition.Does scaling imply superposition?If so why are two different conditions given for linearity?If not can u specify an example for which only one of them is satisfied.

The answer is "no" in general. Let us see the details.

The notion of linearity goes beyond systems. Let us start from vector, or linear spaces. They deal with two objects: vectors and "scalars". In a vector space, one can add vectors, and multiply a vector by a scalar. More formally, a vector space over a field $F$ is a set $V$ together with two operations, vector addition and scalar multiplication. Those operations in $V$ satisfy a number of axioms. The field $F$ possesses two operations as well, satisfying different axioms, namely field axioms, ruling the addition or product of scalars. Those two operations are often assimilated to addition and multiplication. But in general the addition/multiplication in $F$ are not the same as the vector addition/vector scaling in $V$. Indeed, if $f \in F$, $v\in V$, then $f.v \in V$, but is not anymore an element of $F$. The scalar has been cast to a vector. This plays a role in the standard confusion in defining linearity.

Now take for instance the field of rationals $F=\mathbb{Q}$ with the two standard operations. Take $V=\mathbb{R}$ the vector space of reals over $F$ with the standard real addition and real product. To relate this to the introduction, for instance, $1.\sqrt{2}$ is not rational anymore, although the multiplication seems to be the same, because you mix a scalar ($1$) and a vector $(\sqrt{2})$.

Suppose you have a system $S$ that provides the following output: if $v_1\in V$ is rational, $S(v_1) = v_1$; if $v_2\in V$ is irrational, $S(v_2) = -v_2$. For any rational $q$, $q.v_1$ is rational, and  $q.v_2$ is irrational. Resultingly, $S(q.v_1) = q.v_1 = q.S(v_1)$, and  $S(q.v_2) = -q.v_2 = q.S(v_2)$. Hence, $s(q.v) = q.S(v)$ for any $v\in V$, so $S$ verifies the scaling property.

However, $v_1+v_2 $ is irrational. Hence, $S(v_1+v_2) = -v_1 -v_2$, which is different in general from $S(v_1)+S(v_2) = v_1 - v_2$, take $v_1 = 1$ and $v_2 = \sqrt{2}$ for instance.

So scaling does not imply superposition (in your sense) in general.

But there exists somehow converse statements. In other domains, one sometimes calls the scaling "homogeneity", and with additivity we get  the superposition principle for a system $S$: $S(q_1.v_1+q_2.v_2) = S(q_1.v_1)+S(q_2.v_2)$. Keeping with standard fields, it can be shown that additivity implies homogeneity with rational scalars, by playing on $S(v+v) = S(v)+S(v) = 2.S(v)$ etc. To go beyond rational scalar often requires additional assumptions like continuity, which could be troublesome depending on the sets you choose, see for instance Why doesn't superposition imply linearity? Why is homogeneity needed?


Polyphase decimation matlab I am new to DSP, and currently trying to implement a simple polyphase decimation program in matlab.

M is decimation factor, x is input (for now it is a square wave 8096 samples), y_polyDec is result from polyphase, output is result from regular filter then decimation.

code :

FilterBank_LPF = fir1(127, 1/M);
y_polyDec = zeros(8096/M, 1);

for i = 1:M

    p(:, i) = FilterBank_LPF(i:M:end);
    x(:, i) = input(M-i+1:M:end);
    y_polyDec = filter(p(:, i), 1, x(:, i))+ y_polyDec;

end

out_LPF = filter(FilterBank_LPF, 1, input);
output = downsample(out_LPF, M);

%===========================================

If I compare "output" and "y_polyDec" (either comparing the decimated sample values or frequency spectrum using fvtool) I notice there are differences in the values. I was wondering if my code is wrong, or am I still missing something.

Thank you

The problem is in the input signals to your polyphase filters. They should be delayed and subsampled versions of the input signal, like in the following code fragment:

M=8;    % downsampling factor
L=256;  % length of input signal, integer multiple of M
x=randn(L,1);  % input signal
h = fir1(127, 1/M);
yp = zeros(L/M, 1);

for i = 1:M,
    xtmp = [zeros(i-1,1); x];  % delayed input signal
    tmp = filter( h(i:M:end), 1, xtmp(1:M:L) );  % polyphase filtering
    yp = yp + tmp;  % accumulate outputs of polyphase filters
end

% compare with standard filtering
y = filter ( h, 1, x);
yM = downsample( y, M );

n=1:L/M;
plot(n,(yM-yp)/max(yM))


IDFT - influence of sample numbers I'm quite a newbie in dsp domain, and there's something I can't understand clearly.

I'm doing a IFDT of a spectrum regularly spaced with complex values : the spectrum is a complex conjugate symetric frequency domain sequence so the output of the IFDT should be real.
It is when using a sampling rate of 50MHz with 2048 samples.
But when I change the number of samples (for exemple 1024), it does not work anymore.

Can someone explain why it won't work ? I'm clearly missing something...

How are you changing the number of samples ? If you are just truncating the 2048 samples $X$ to 1024 samples $X^t$, then the samples won't be complex conjugate anymore.

First sample of the 2048 length sequence will be considered as DC and then the remaining should be conjugate symmetric around $N/2$. That means, $X[k] = X^*[2048-k], \forall k \in \{ 1, 1023 \}$ and $X[1024]$ must be real. Then IDFT will give real $x[n]$.

When you want to downsample it to 1024 samples, then, first of all there will be aliasing in time -domain but still you can get real time domain samples if the downsampled 1024 length $X^t$ maintains the above explained conjugate symmetry.


Unipolar Transmission Spectrum I am currently trying to understand mathematically the amplitude shift keying modulation technique. My question is, say I wanted to compute the Fourier series of a basic pulse train, it is quite simple since it is periodic. But in digital transmission (Unipolar) the signal is random, how does one compute the spectrum of this transmission. I hope my question makes sense.

Here is a snapshot of some notes I found online. The data signal is a random signal of 1's and 0's. How can anybody interpret the spectrum of it? Thanks.

The spectrum of a pulse train with a defined pulse with T that is indeed random in duration of zeros and ones would have the magnitude envelope of a Sinc function (PSD with envelope Sinc^2) with first null at 1/T.  Since the data is random the spectrum will be continuous. Similar to this is something like GPS where we transmit Pseudo-random patterns, but the pattern repeats (assuming no data modulation which again adds randomness). In this case, whenever there is a repeating pattern, you would see discrete tones at the repetition rate and every harmonic of the repetition rate, and the envelope would still be the Sinc function as described with nulls at 1/T. (Thus with doing forensics on such a spectrum you can discern the length and repetition rate of the sequence based on the spacing of the spectral lines, and the number of "chips" in the sequence based on the location of the null, etc.).

Finally multiplying by the carrier just moves this spectrum (which is best viewed in my opinion as a positive and negative spectrum even at baseband, as then you can translate it from there to any frequency with no changes to its features).  The fact that this is on-off keying means that it will have a DC component, so the carrier will come through very strong at the RF output. In comparsion if the signal was bipodal (which would then be BPSK modulation), with an equi-probable stream of 1's and 0's you would have no DC at baseband (avg =0) leading to supressed carrier at RF.  The carrier power in the former case is wasted (no information added) so is less efficient in terms of BER vs power than the latter case.

The reason for this is the base square pulse has a FT that is the Sinc function, and your random data is represented as a stream of impulses (with random positive and negative sign) that is convolved in time with your pulse (effectively a zero order hold). If the impulses are random in a pattern that never repeats, the spectrum of that pulse stream will be a continuous flat spectrum: white noise.  Therefore in frequency you multiply white noise with the Sinc to get the spectrum before being upconverted. The mixer translates this spectrum to the carrier.

Your picture just shows the main lobe, but if your pulse is indeed rectangular, you will have all the sidelobes of the Sinc function out to the effective BW (based on the rise/ fall time of your pulse).

I hope this clears it up for you!


Performing root-music, how do i retrieve frequency and phase data of a signal? I'm trying an alternative method to the classical fft to do signal analysis. I perform root-MUSIC on a certain known signal + noise in matlab with the function 'rootmusic'. This gives gives the signal in radians/sample. This can then be converted to a frequency of a signal and at first sight it looks promising.
My question is wether I also am able to retrieve the phase of the signal from this function or use the frequency output to in turn obtain the phase of the signal at that particular frequency? I'm not an expert at what perfectly goes on inside this function but having looked into it, it seems that this methed solely provides a frequency and no phase data. Anyone care to explain further details of what or what is not possible with root-MUSIC or if there is a better/higher resolution alternative method to fft? At first glance I tried root-MUSIC as it is suggested as a quite popular substitute.

Subspace-based spectrum estimation methods give you a high resolution pseudospectrum that is used for frequency detection only. The complex amplitude (that includes the phase) can be calculated with linear least-squares once your frequency estimates are available.
In MUSIC and root-MUSIC,

You calculate the autocorrelation matrix of the sensor signals.
Calculate its eigenvalue decomposition. By looking at its largest magnitude eigenvalues, you can estimate how many frequencies are there in the signal.
Separate the noise subspace (i.e, the noise eigenvectors associated with the low magnitude noise eigenvalues) from the signal subspace.
Project a linear search-vector of candidate frequencies onto the noise subspace. The frequencies that are actually present in your signal will have 0 inner product with the noise eigenvectors since they are uncorrelated. So, the inverse of this projection (known as the pseudospectrum) will have peaks at those frequencies.
root-MUSIC is a more computationally robust because it models the pseudospectrum as a polynomial evaluated on the unit circle. The roots of this polynomial are our desired frequencies.

Okay, after this high-level description of the two methods, I want to discuss how you can calculate the complex amplitudes once you have the frequencies. So you have $N$ samples of a complex signal, $x$, composed of $L$ frequencies plus added noise $w$,
$$x(n) = \sum_{k=1}^{L} A_ke^{j\theta_k}e^{j\omega_k n} + w(n) \\
 = \sum_{k=1}^{L} \bar{A_k} e^{j\omega_k n} + w(n)$$
In vector form, this can be written as
$$\begin{bmatrix} x(1) \\ x(2) \\ \vdots \\ x(N)\end{bmatrix} =\begin{bmatrix} e^{j\omega_1 1} & e^{j\omega_2 1} & \cdots & e^{j\omega_L 1} \\
\vdots & \vdots & \ddots & \vdots \\
 e^{j\omega_1 N} & e^{j\omega_2 N} & \cdots & e^{j\omega_L N}
\end{bmatrix} \begin{bmatrix} \bar{A_1} \\ \bar{A_2} \\ \vdots \\ \bar{A_L} \end{bmatrix} + \begin{bmatrix} w(1) \\ w(2) \\ \vdots \\ w(N)\end{bmatrix} \\
\mathbf{x} = \mathbf{\Omega \bar{A}} + \mathbf{w}
$$
Now, we know $\mathbf{x}$ and we have estimated the frequencies $\omega_k$ from root-MUSIC, so we know $\mathbf{\Omega}$. Now, the least squares estimate of $\mathbf{\bar{A}}$ is
$$\hat{\mathbf{\bar{A}}} = (\mathbf{\Omega^\top \Omega})^{-1} \mathbf{\Omega^\top x}
$$
Of course for this to work, $N > L$, which is also a requirement for root-MUSIC. Now you can estimate the phase angles by calculating the angle of $\hat{\mathbf{\bar{A}}}$.


Fast empty page detection As part of a duplex scanning pipeline I need to detect empty pages such that they can be excluded from the OCR and the resulting PDF document.

Currently, I'm using a 'blurry' trim which works pretty well. Example with ImageMagick:

convert image-0001.png -shave 300x0 -virtual-pixel White -blur 0x15 -fuzz \
   15% -trim info:

(if the resulting dimensions are below a threshold then the page is considered empty)

But this is slow.

What are faster alternative approaches for accurately detecting empty pages in scanned images?

Things that complicate the empty page detection:

the noise - i.e. scanned dust etc.
margins/punch-holes that may show up as black lines
almost empty pages (e.g. just 2 words in 12 pt on an otherwise blank pages) must not be detected as empty
pages without any text but with images must not be detected as empty

What are faster alternative approaches for accurately detecting empty pages in scanned images?

OK, convert is universally known to be rather hungry for RAM and potentially rather slow (due to its internal image model handling). My guess is that if you don't use it but, for example, simply do your whole detection within a scripting language with an image processing library (e.g. Python3 with PIL), you'll be fast enough for this to be a viable approach, as you, as a nice addition, won't be constantly saving and loading images from disk, but keep them in RAM until you're done, which is orders of magnitude faster.

I like the blurring approach; it has the potential of removing noise sufficiently well. But:

Why don't you simply

threshold to convert to black&white (might want to do that adaptively based on the average brightness)
erode by a few pixels, just enough to kill dust and noise,
count the remaining black pixels, and
declare clear page based on a threshold?

Since punch holes and margins will be larger in black area than two words, I'd previously try to just exclude them with a simple static mask, if your raw material permits (i.e. if you're scanning shelves of filed documents, then all will have the punch holes in the same place), or by trying to "smartly" detect them (for example, look for circles in potential punch-hole places, diagonal lines in potential staple places and so on).


Is there a closed form solution for a temporal calibration of two different signals? Given two signals $f_1(t),f_2(t+dt)$ where $t$ is time and $dt$ is the time delay between the two signals,

is there any closed-form solution with respect to $dt$?
what are the efficient global solutions to this problem?

In my problem, the function is $f\in SO(3)$ but I can reduce it to $f \in R^3$ or just look at only one dimension.

I would do the analyis on your three dimensions separately.  If the answers you get agree, then you have your confirmed answer.  If processing needs to be minimized, then you can find the delay on one, and confirm with the others.

Since you have sensor data, it probably contains some noise as well.  Heavily smoothing your data (I like exponential smoothing for its efficiency) will mitigate the effects of noise, reduce the number of local maxima, and allow you to do a coarser correlation measure.

There is some dispute on the proper definition of correlation.  The one you should use for best results is:

$$ c = \frac{\vec x \cdot \vec y}{\| \vec x \| \| \vec y \| } $$

Where $\vec x$ is a subset of the smoothed $f_1$ and $\vec y$ is a subset of the smoothed $f_2$ for the two intervals and sampling density your are comparing.  Each vector should be shifted so the elements have a mean value of zero.  Take the average of all the element values and subtract it from each value.

First do a coarse brute search.  Select every 10th (or whatever spacing you choose) value from your functions to build your vectors.  Then calculate the correlation coefficient for a range of delay values, stepping by a coarse amount.  This will give you a coarse delay value.

Repeat this process, centered at your coarse delay value, using every sample and stepping by a single sample size.  You should only need to slide plus or minus half your coarse spacing side.  This will give you a fine delay value.

If you need subsample resolution, use the minimum value you found in the fine search and its two neigbors, model them as a parabola and find the minimum point there.

Hope this helps.

Ced


why always structured Images are considered as 2D signals? I am dealing with 2D signals that are random. To be more specific my input is a matrix of 1's and 0's. Moreover, at any position (n1,n2) the probability of having a 1 or a 0 is same. One example is shown below

When I look at books dealing with 2D signal processing, everything is explained wrt to image processing. In my understanding, image is a very special kind of 2D signal where one can expect to see some relationship between the pixels like EDGES. Such features are useful when doing image restoration,etc.

These books start with defining 2D signals and system but explain everything wrt to image processing application which is a bit restrictive.
Why is the case?

A 1D signal that is only noise with no structured feature to be seen is of no interest.

A 2D signal that is only noise or a superposition of uncorrelated 1D signal rows, with no 2D structured feature to be seen, is of no interest.


Decimating by last interpolating factor? (Multi-stage) Let's say I have a multi-stage interpolator at the transmitter, which interpolates by 2, 4 and 8 (that is, 64). When decimating (at the receiver) why should I start decimating by 8, then 4, then 2, instead of 2, 4 and 8? I mean what's the reason behind? What makes it better than doing it the other way round?

The idea is that you'll be saving computational steps in the filtering part of the decimation. Think about it, if you have 800 samples then downsample by a factor of 8, now you only have to low pass filter 100 samples for the first stage. If you instead tried to downsample by a factor of 2 first, then you'd have 400 samples to filter after the first stage. The order matters because of the number of computations taken in the filtering step.


Model necessary for PID control As I am not from control theory, I do not understand - after reading the Wiki page about PID control (https://en.wikipedia.org/wiki/PID_controller) and watching some videos - whether a model of the process which is due to be controlled is necessary for a PID control or not. For example it is written on one occassion:

  Control action – The mathematical model and practical loop above both use a "direct" control action for all the terms, which means an increasing positive error results in an increasing positive control output for the summed terms to apply correction.

And on another occassion it is written:

  But the PID controller is broadly applicable, since it relies only on the response of the measured process variable, not on knowledge or a model of the underlying process.

If I want - for example - to control a heating system using a PID-control such that the temperature in the building in at a certain value, is it necessary to have a model of the heating system or the building that maps the fuel consumption (e.g. gas or oil) to the generated heat and the temperature?

Update: What does this transfer function (mentioned below by Marcus Müller) consists of? He wrote " I'd argue that you want your plants transfer function to be continuous, and that's already a model"

Update 2:
I think my question is not too dificult and people with experience in that field should be able to answer it. Does a PID-controller need a model of the process? Unfortunately this question has not been answered yet :-( –

This was my recent experience where I found it was absolutely necessary to use PID control versus the PI control I would typically use: I was developing a TEC (thermo-electric cooler) oven control circuit that had extremely low frequency poles (sub Hz). With a typical PI control the response time would be limited by the lowest frequency pole. Such that in my case the resulting system with the dominant lower frequency pole would take minutes to settle: this can be estimated using the relationship for a first order system of the 10% to 90% settling time in seconds given by $t_r = 0.35/BW$ where $BW$ is the frequency in Hz of the pole. If a single pole is dominant a first order approximation will be reasonably accurate.

With the use of a PID solution, the achievable rise time can be extended out to the next higher pole, and in my case this resulted in a TEC control of under 10 seconds where the plant (TEC) on its own would take minutes.

In practice what you observe if how the "D" enables a very large overshoot of current (the TEC is current controlled), so as long as the current source can supply the needed current, the temperature is very quickly changed much more during transitions than the current that would be commanded under PI control only. This makes sense since the D is reacting to the differential change of the input.

So in summary I am giving one example where it was justified in using a PID control and what the decision factor was: I was faced with a plant that had two lower frequency poles and wanted to have a response time and settling beyond what would be limited by the lowest frequency of the two. By using PID, the response time and settling was then limited by the higher of the two frequencies resulting in a significantly improved response time.

For your last question, no a PID controller does not need a model of the process as you can experimentally determine the gain coefficients in most practical cases (specifically I outline an approach below for systems dominated by lower frequency poles). The process to experimentally set the gain coefficients is as follows:

Set P, I and D all to be zero and then monitor the step response (step the input and measure the output) while slowly increase P only.

Increase P to the threshold of instability-- this does not mean ringing but a sustained oscillation where the amplitude neither further increases or decreases. Once determined, divide P in half and use that as the gain value for P.  There will be expected offset errors which is why we need I.

With P set as determined above, slowly increase I to the threshold of instability. Then divide this value in half. With a typical pole dominant system such as I have experience with, the rise time will be quite fast but there will still be a significant ring down and here we will see the benefit of adding D instead of simply reducing I to eliminate the ringing (which would also correspondingly increase the rise time).

With P and I set as above, slowly bring up D to achieve the desired step response balancing an overshoot with rise time, and final settling time. You will see that the rise time as set above will not be impacted but the ringing will be eliminated. The ideal transition that has the best compromise will typically have a 8% to 15% overshoot. At this point all values can be adjusted from these starting points to optimize given the actual response of the system.

This is a general approach that works with many real world systems dominated by low frequency poles. This will not achieve good results with a more complex system with many poles and zeros in close proximity or an unstable plant with poles in the right half plane, but it is quick to test and determine the practicality of proceeding with this experimental approach.

However, it is certainly very practical and useful to have a model of the process and be able to simulate the operation. Nyquist's approach of modeling control loops in the frequency domain is very practical for black box systems in that we needn't have a mathematical model of the plant in terms of polynomial transfer functions or pole zero locations to mathematically determine stability using the Nyquist criteria. What we do need is the magnitude and phase of the plant which we can do experimentally (when the plant itself is stable, meaning all poles in the left half plane) by sweeping the plant with a sine wave and measuring the magnitude and phase of the response compared to the input. This is the transfer function. Having that alone is sufficient to design the control loop and determine all stability factors (again when the plant is stable on its own, otherwise things get a lot more complicated which I didn't cover here). Further what is great about this approach is the transfer function needn't be limited to polynomials with pole and zero locations which means exponential factors such as true time delays can be easily incorporated and stabilized. The system may be too complex or cumbersome to try and obtain an accurate mathematical model and here we see that simply having the magnitude and phase of the black box system will suffice. This was the great achievement in Nyquist's approach as to its practical application.

For more on the Nyquist approach to determining stability see Nyquist plot interpretation when curve hits the origin


convolution and groupdelay I have been given a sinusoidal signal x in time domain,asked to investigate the frequncy responce of the signal using matlab and observing the result,and using grpdelay(x,1) command in matlab to observe the group delay of the signal .and taking y=conv(x,x);which define the convolution of x signal by itself and then asked to use the grpdelay(y,1).
My question here is what change should I exspect to see in group delay figure of conv(x,x)?
or better to say how convolving of the signal by itself would effect the group delay of the signal.
Thanks

There is some confusion here. Frequency Response and Group Delay are properties of a System not the properties of a signal. You can only calculate the spectrum of a signal (magnitude and phase). A "signal" doesn't have group delay unless it's the impulse response of an LTI system.


Bitmap border/stroke alogirthm? I am looking for an algorithm that adds borders to a bitmap image(very much like the photoshop stroke effect) but one that will also have antialiasing to the borders(or one that copies alpha from the original bitmap to make the border look antialiased).

If you can't answer with an algorithm a direction of creating one myself would be great.

*****Update***

The part coloured black is the original ellipse and the part coloured green is the stroke.

The stroke is basically adding a contour to the bitmap or edges.

The Euclidean Distance Transform can produce dilations and erosions with suitable parameters and filtering.

The algorithm that Photoshop uses, and the one that is best suited for stroking, is to calculate the Euclidean Distance Transform in integers and without taking the square root (i.e. calculate distance squared). This can be made extremely fast using the technique from Meijster.

To get anti-aliasing, just use the fractional part of the number after taking the square root.

For example:

if (distanceSquared < radius * radius)
{
  (fill solid pixel)
}
else
{
   distance = sqrt (distanceSquared);

   if (distance < radius + 1)
   {
     alpha = 1 - (distance - radius);

     (blend pixel at alpha)
   }
}

This techinque is used in my open source library:

https://github.com/vinniefalco/LayerEffects.git


How to find analytic description of filtered signal I am looking for an exact analytic description of a filtered signal. I have an electronic circuit whose input is a monoexponential decay. First (1) the signal gets filtered by a simple RC-Lowpass. Then (2) a second order Bessel Low Pass Filter is applied.

The respective transfer functions are:
\begin{align}
TF_{\rm RC}(s) &= \frac{R}{s+\frac{1}{R C}}\\
TF_{\rm Bessel}(s) &= \frac{a_1}{a_1+b_1s+s^2}
\end{align}

The input signal is:
$$I_{\rm in}(t) = \frac{V}{R_1 + R_2} + R_1\cdot V_{\rm cmd}\cdot e^{-t\frac{R_1 + R_2}{CR_1R_2}}\frac{1}{R_2(R_1 + R_2)}$$

My idea was to calculate the output in the $s$-domain:
$$ V_{\rm out}(s) = I_{\rm in}(s)\cdot TF_{RC}(s)\cdot TF_{Bessel}(s) $$
and then apply the inverse laplace transform to obtain $V_{\rm out}$(t).
The inverse transform is what troubles me. I tried to calculate it with MATLAB, but it doesn't return an analytic function. I guess this comes from the complex poles of the bessel filter. Could you give me a hint on how to calculate the output?

Edit:
$$ TF(s) =  \frac{6.092e12}{s^3 + 5.367e04*s^2 + 1.038e09*s + 6.074e12}$$

$$ I_{in}(t) = \frac{10e-3}{500e6 + 10e6} + 500e6*10e-3*e^{-t*\frac{500e6 + 10e6}{10e-12*500e6*10e6}}\frac{1}{10e6*(500e6 + 10e6)} $$

$$ TF(s)*I_{in}(s) = \frac{4.9e-21*(1.2e24*s + 2.5e26)}{s*(s + 1.0e4)*(s^3 + 5.4e4*s^2 + 1.0e9*s + 6.1e12)} $$

PFE yields:

$$ TF(s)*I_{in}(s) = \frac{1.1e-9 - 4.8e-11i}{s + 2.2e4 - 1.1e4i} + \frac{1.1e-9 + 4.8e-11i}{s + 2.2e4 + 1.1e4i} + \frac{6.6223e-7}{s + 10200.0} + \frac{-6.6441e-7}{s + 10235.0} + \frac{1.9666e-11}{s}  $$

Inverse Laplace transform of the upper term will yield a complex result. I know that I can simplify the first two terms by bringing them to a common denominator and thereby get rid of the complex numbers. The inverse transform is however still complex.

Edit: calculation with symbolic parameters yields:
$$ H(s) = \frac{a*s+b}{c*s^2+d*s+f} $$
$$ H(t) = a*dirac(1, t) + (2*b*sin((t*(- d^2 + 4*c*f)^{1/2})/(2*c))*exp(-\frac{d*t}{2*c}))\frac{1}{(4*c*f - d^2)^{1/2}} $$

I don't quite understand where the dirac comes from, but this is basically the same result as doing the calculation with complex numbers and the taking the real part.

Assuming your algebra is correct, there is nothing wrong with complex numbers in the Laplace domain. Complex poles reflect an oscillatory signal - there is no oscillation at all if all poles are strictly real, just exponential growth or decay. The variable $s$ is itself a complex number, and the transfer function is defined over all $s$, i.e. the entire complex plane.


FIR difference equation vs convolution summation I'm sorry if my question is too simple, but I am self-studying FIR because we have to report it in our class. I am using the book of Proakis and Manolakis as reference.

Under 7.2 Structures for FIR systems:

  In general, an FIR system is described by the difference equation
    (7.2.1)

  or evidently, by the system function

   (7.2.2)

Under 7.2.1 Direct-Form Structure:

  The direct-form realization follows immediately from the nonrecursive difference
  equation given by (7.2.1) or, equivalently, by the convolution summation

   (7.2.4)

I can see that they almost look the same apart from the bk in 7.2.1 and h(k) in 7.2.4. Question is, how are they different or related to each other? Sorry if you find my question too basic. Thank you in advance for the help!

I should have read further before asking this question! Anyway, I've found out that there are different ways to implement the FIR system as mathematically described in Equation 7.2.1 in the post above. Among these realizations/structures is the direct-form, which is the simplest to begin with. Equation 7.2.4 shows the transformation of the general equation in 7.2.1 in accordance with the direct-form.


Calculate symbol rate for QAM modulation Not sure if this is the right place to ask, feel free to move it if any other site suits this better.
Anyway I am referring to QAM modulation used in US Cable TV. I know it comes under Annex B which is related to FEC.
I know that 64 QAM has a symbol rate ~= 5057000 Symbols/s and 256 QAM 5360000Symbols/s
I would like to know how this is calculated, also what symbol rate will, say 128 QAM have for example?

There is no inherent symbol rate for any modulation type.  64 QAM symbols transmit 6 bits per symbol, and 256 QAM symbols transmit 8 bits per symbol.  If the bit rate were to stay constant, then the 256 QAM symbol rate would be 6/8 of the 64 QAM symbol rate, but it appears that that is not the case (assuming your numbers are correct).  That may be because they are using more forward error correction (FEC) to mitigate the higher power requirements of 256 QAM.

Given that the bit rates are not the same, we can't guess what the 128 QAM symbol rate would be.

The general equation for a given modulation type and its symbol rate is:
$symbols/s = \frac{bits/s}{bits/symbol}$

Determining the bits per symbol is usually pretty easy.  It's log2 of the number of possible symbol points.  Thus, the bits per symbol of 256QAM (which has 256 possible points) is 8.

The trick is finding the bits per second.  Wireless transmitters often have the following blocks:

Data source => Compression => Scrambling => FEC Encoding => Pulse Shaping => Modulation

To know the final bit rate you need to know the bit rate of the data source, how much the compression reduces the data rate (if there is compression), and how much the FEC encoding increases the data rate.


How to Select Point Spread Function Empirically for Image Deconvolution? When the captured image is blurring, one way of obtaining a clear image is via image deconvolution technique. In order to perform deconvolution successfully, usually we need to pay attention to the following issues:

Image content, which determines what kind of apriori information you will use (Total Variation, Laplacian and so on) during the image reconstruction procedure.
Image noise level, which determines the balance factor between apriori information and data consistence term.
Point Spread Function (PSF), which is used to measure how blurring the image is.

Here my question is related to PSF estimation. The simplest way of finding the right PSF is to use different PSFs to perform deconvolution, and then select the best one visually. Then my question is how we can tell which PSF is the best visually.

One way, though not directly visual, is to observe the image statistics.
We have a pretty good idea about the statistics of Natural Images, more specifically, their Gradient Distribution (See Statistics of Natural Images and Models by Jinggang Huang, D. Mumford, What Makes a Good Model of Natural Images by Yair Weiss, William T. Freeman, The Statistic Distribution of Image Gradient - Mathematics StackExchnage and UCLA Stat232A-CS266A Statistical Modeling and Learning in Vision and Cognition - Chapter 2 - Empirical Observations: Image Space and Natural Image Statistics).

So what you can do is see how well the statistics of the result matches the models.
Actually many modern (Prior to the Deep Learning era) Deblurring Methods work according to this approach with very nice results.

Look at Blind Motion Deblurring Using Image Statistics by Anat Levin.


how to deal with signal to make it suitable for kalman filter I want to build a navigation system with gyroscope/odemeter and gps. I have konwn there is a famous method called kalman filter for this application. And now, I have done some experiment in matlab with kalman filter to get hide state vector which included vehicle's current postion/velocity/heading etc. But out of my expection, the result is not very good even more bad than origin data which comes from gps. May be you could guess the reason why they are so so bad. Yes, gyroscope and odemeter have uncertained error also as gps, so that's why you use kalman filter to deal with them. In some case that's useful enough.
But in some scenario case that's not working. Because the character of original data can't match with requriment of kalman filter. As I probably know, it suppose that filter works under the condition of which your data noise is white gaussian or some other types. so many artical tell me you must have good knowledge to character of this noise, but on the other hand they don't tell me how to know it. that's just my problem and ask you in here. how to do it and how to deal with it for kalman filter?

Kalman filter have basic requirements that the random processes like input signal, measurement signal,process noise and measurement noises all should be white Gaussian in nature.what Kalman filter actually do is that it uses the system model and input signal and generates its own output which is actually a predicted value now using the observation from the system it calculate the correction to be applied to the predicted output and hence final results are better than observations obtained from the system. Hence very first step- is that you need a system model. i am totally disagree with @Deniz notion that in non-linear case you should use particle filter. particle filter is used when the random process as mentioned above are not white Gaussian. Even if your system is non linear you can use Extended kalman filter or unscented kalman filter. As ur case is concerned i think you should apply extended or unscented kalman filter as velocities for navigation will keep on changing.


Frequency components much less prominent at the receiver for multiple tones I am trying to transmit a multi-tone signal from one USRP (N210) to another and have written an out of tree block to generate the multi tone signals using an ifft. The two USRPs are connected directly via a cable and a $50\ \Omega, 20\textrm{ dB}$-attenuator.

I am finding that the signal is received clearly for 1 and 2 tones but not for 3+ tones. I'm not sure what I can do digitally to make the frequency components are lot more distinct. The received signal is amplified by $10\textrm{ dB}$ and passed through a LPF and I've attached some FFT plots showing the results.

I don't have enough reputation to show the transmitted signals but they are just direct spikes evenly distributed through the bandwidth and very distinct.

The received signals at the USRP are complex and so the FFT plots shown are only the real components.

Flow graph

Relevant graphs for a series of different tones

Single tone

Receive:

Transmit:

Two tones

Receive:

Transmit:

Three tones

Receive:

Transmit:

The number sweep block is an out of tree module that produces the frequency array that is passed into the FFT block.

The two USRPs are connected via a MIMO cable to ensure synchronisation of their clocks.

The problem was from the fact that I was using a rectangular window for the ifft at the transmitter.

Using any other window at the ifft ensures there are no harmonics outside the 200kHz (which are interrupting the signal).

Here is an example of the signal using the rectangular window

And using a hamming window

Four tones

Received

Transmitted


For filters with impulse responses on a specific $x[n]$, how do you calculate DFT of $y[n]$? Basically I'm talking about N-point moving average/difference filters and not one where $\delta[n]=k$. For example given say a two-point moving average filter:

$$h[n] = \dfrac{\delta[n]+\delta[n-1]}{2} $$

and a defined input signal like

$$x[n] = \sin[2 \pi n]$$

$x(n)$ is discrete and the window size = the period size.

How would you find $Y_k$ (the coefficients of the DFT of $y[n]$)?

How would you find Yk (the coefficients of the DFT of y(n))?

By applying the DFT to Y.

So: math.

Mathematically, you filter something with a FIR by convolving the signal with the taps vector:

$$y = x * h$$

Hence, the DFT of $y$ becomes is:

$$\mathcal F \{y\} = \mathcal F \left\{ x * h \right\}$$

and thanks to the properties of the DFT, convolution in time domain becomes multiplication in frequency domain:

$$\mathcal F \{y\} = \mathcal F \{ x\}  \mathcal F \{ h \}$$

So you just apply the DFT to both your signal and your taps vector and multiply.

Seriously, this is very basic DFT theory, and you will learn that if you continue your studies of the discrete Fourier transform.


Determining sampling rate for signal reconstruction I want to know how to determine the signal sampling rate required to reconstruct electrical signals whose frequency components are unknown, such as, for example, signals of charge/discharge voltage on a capacitor. I want to design a system capable of scanning this type of signal in real time.

Sampling a signal requires that it should be bandlimited. When it's not, you should apply a preprocessing analog bandlimiting (aka anti-aliasing) filter to enable its proper sampling. There are various concerns on how to choose the most practical anti-aliasing filter. The most important concern is the filter cutoff frequency which determines the minimum sampling rate to record the signal.

Your application is about capturing capacitor charge/dischage signals which are of spiky nature and hence posses very large bandwidths; quite large sampling rates might be necessary. In very rough terms, the dischaging signal's bandwidth is inversely proportional to the time-constant of the discharging network, assuming a first order simplified model.

If you know this time constant, then you should select the lowpass filter whose time constant is smaller than that. This will make sure that the anti-aliasing filter's bandwidth will be larger than the signal of interest's bandwidth.

So you should firt determine the time cosntat of your observed network, to predict a first order estimate of the minimum sampling rate required.


Criteria to distinguish FIR and IIR filters from difference equation? What criteria should be used to safely decide if it is an IIR filter or FIR filter from a difference equation? FIR filter are always stable, meaning all poles are inside the unit circle AND have pole-zero cancellations i.e have equal number of poles and zeros at the same point in the z-plane. IIR filters are conditionally stable, i.e always have poles on/outside the unit circle. Is this good enough to distinguish from the difference equation? what other criteria should be considered?

PS: Please don't say IIR filters have recursive difference equations, and FIR doesn't. It is not true.

If the Z-transform of the feedforward section is divisible by the Z-transform of the feedback section, the filter is FIR.

Consider your example: $y[n] = y[n-1] + x[n] - x[n-3]$.

The Z-transform is $\mathrm Y(z)- z^{-1}\mathrm Y(z) = \mathrm X(z) - z^{-3}\mathrm X(z)$, and the Z-transform of the response is $\mathrm H(z) = \mathrm Y(z)/\mathrm X(z) = (1 - z^{-3})/(1 - z^{-1}) = 1 + z^{-1} + z^{-2}$. This means that the filter can be realized as $\mathrm y[n] = \mathrm x[n] + \mathrm x[n-1] + \mathrm x[n-2]$, which is clearly FIR.


What does the visual representation of matched filtering of images look like? So I am using the Source Extractor to find point sources in my images and it offers the possibility to use matched filtering to enhance the detection results. In the documentation they also claim that the matched filter defaults to a convolution of the data, should the noise be equal across the field (please correct me if I got this wrong). While there is a nice and easy to understand visual reprentation for the convolution, where the convoluted pixel value is simply the sum of the pixel values "below" the kernel weighted by the factors in the kernel "above", I have yet to find a similar explanation for what is going on in matched filtering.
I also like how the matched filter is derived in terms of matrices, since I feel like this makes the most sense with image processing (unlike the derivation using itegrals I found in many places), but I haven't found a book yet where this is well explained. Can someone point me to a source where matched filtering is explained regarding image processing?

There's really nothing special about image processing here: the matched filter is still the hermetian, spatially-inverse filter (=convolution kernel) to the original filter (=convolution kernel).
Since images tend to be real-valued, most convolutional kernels are, too. In which case you just flip the kernel along both axes and are done.
If the original kernel is symmetrical along both axes, it's its own matched filter.
The explanations used for the one-dimensional matched filter are the same as for the 2D case:
Cauchy-Schwarz inequality says that, for any given filter energy, a filter that fulfills the above criteria leads to the highest SNR in the additive uncorrelated noise case.


Help with automated baseline removal This is more of a conceptual question, so apologies in advance if this is not the correct forum for such questions.

I am trying to find a repeatable and consistent method for removing the baseline of a large number of independently acquired spectra. After an extensive literature search, I am well acquainted with the various methods for baseline removal - including the favored ones for the particular spectra I am dealing with.

However, the method I have used so far (asymmetric least squares) requires me to set some hyperparameters which do not produce consistently good results across the entire dataset. I have tried variations of lambda and p, some of which produce better results on certain spectra, but no combinations as yet produce consistently good results throughout the dataset (>10,000 spectra). The application of the modelling project is such that manual tuning of these hyperparameters for individual spectra is not an option when putting these models into production.

Does the community have a view on a baseline removal approach that is more appropriate for my application?

Baseline removal for signals and background filtering in images  remain, as far as I am aware of, an ill-posed  problem related to source separation. Reasons are legion, mostly around the unresolved question: "how to measure the range of admissible baselines both qualitatively (for the signal specialist) and quantitatively (for the objective)?".

If from one measurement at a time, you want to separate two or three components, this is an under-determined problem. Using asymmetry already induces at least one asymmetry parameter. Using any prior or constraint (gaussianity, positivity, peak symmetry, sparsity) is likely to yield one more parameter each, generally linking two quantities of different dimensions. And even if my objective target is clear, which is seldom the case, many solutions could be valid. For instance, if I am interested in  peak area, the baseline below is fine, yet shocking for the data analyst:

At the present time I see two main options:

try to fit baselines globally, not individually, using optimization and algebraic tools (low-rank, NMFetc.) but as numerical data has limited precision and quantization (integer, float), some constraints have to be (machine-dependent) fixed,
try to learn baselines of "meaningful subsets", but you have to choose their size, how to weight them, etc. Even deep learning needs a little tuning,
try to guess parameters from signal-based heuristics, and a little grid searching.

I admit I still stick to the third option, not yet successful in most instances. Check  How to remove or filter the drift problem in measured Strain signal? for a related question and  answers.


What is Audio codec in simplified term? I am trying to understand what a codec actually is

Wikipedia says

  In hardware, audio codec refers to a single device that encodes analog
  audio as digital signals and decodes digital back into analog. In
  other words, it contains both an analog-to-digital converter (ADC) and
  digital-to-analog converter (DAC) running off the same clock. This is
  used in sound cards that support both audio in and out, for instance.

So this implies codec is something that comes with both ADC and DAC converters (so a codec should have an analog input with a corresponding digital output and a digital input with a corresponding analog output)

From my limited understanding, codec should have some kind of coder and a decoder, that is, for whatever digital representation I get from a particular analog input, the codec should should give me a different, encoded version of this digital representation, and it should also be able to reverse this encoding (since it's a codec, unless I got it wrong and the digital equivalent itself is considered a type of encoding?)

However, this doesn't seem to make sense to me in terms of the DSP that I am working with

The following figure is from the datasheet of the CSR8670 Bluetooth Audio chip

In a lot of a places, the datasheet, other documents and comments in the sample programs mention the audio input ADC as a codec, so I assumed codec here meant the wikipedia definition

Looking at this picture, we can see that it shows ADC and digital codec to be two separate entities. Moreover, both the inputs from the analog mics as well digital mics go through the digital codec

What exactly does this mean? I can use the debugger to check the DSP memory for the numerical representation of whatever audio that comes in through the microphone. Does this imply that the numerical representation that I see in the DSP memory is not the same as what was generated by the ADC?

What is codec here, in very simple terms?

As Marcus have pointed, the codec in that article refers to a narrower meaning of the term codec specifically applied to digital (computer) hardware systems.

The term codec is acronym for coder and decoder, hence a codec consists of two separate parts (at different ends) working together to perform an information processing task in between.

A coder is an algorithm (implemented in digital hardware, software or analog hardware) that is used to convert a source of information into a form that has some advantages for the purpose of representation, storage or transmission etc. The source after being processed by the coder is said to be as encoded. The most typical application for coders is compressing digial media such as music, speech, images, video so that they can be represented with less number of bits then their raw formats which is advantageous for storage space and transmission bandwidth for a fixed transmission time such as in real time communications or equivalently transmission time for a given fixed bandwidth for file downloading from internet. The algorithm of the coder is specified mathematically, and it can be implemented in software such as those dll files in a Windows operating system or inside specific program folders. It can also be implemented in hardware which is the case for such as the digital cameras where the hardware device is required to convert the captured raw video input stream into a chosen video format for storage purposes as fast and efficient as possible. A software codec running on a general purpose CPU could be weak for such embedded applications. Note that historically analog hardware base encoders were also used for transforming raw media analog inputs into feasible analog formats, the prime example being a TV encoding IC which converts RGB raw analog video stream into an NTSC/ PAL/ SECAM video format which enables a low bandwidth and black and white compatible transmission over RF broadcasting.

The decoder does the inverse; given an encoded data stream, it produces the original raw data stream as accurate as possible, exactly if the codec is a lossess type and approximately if lossy type. The decoder can also be implemented in hardware or software. The hardware implementation is always faster wheres the software is always more flexible.

As a general principle, the coder and decoder pairs are designed such that the encoders are more powerful and expensive, so that decoders can be less powerful and cheaper (form a commercial perspective) which is the case in multimedia broadcasting systems.

Coming to your figure, it's an accepted misusage of the term codec. Where the input audio from the digital microphone (which is another accepted misusage) is converted into a suitable form for representation and transmission in the following digital systems and vice versa. Hence conventionally an ADC chip which converts analog audio into digital PCM stream and a DAC chip shich converts digital PCM stream into analog audio output may be termed as a codec in the digital audio industry.


How to calculate biquad coefficients for the partial (1%, 10%) applying of pre-emphasis effect I'm using sox to emulate RIAA pre-emphasis and de-emphasis effects on audio files. For pre-emphasis I call it like this:

sox --plot octave in.wav out_biquad.wav biquad 1 -1.700724 0.7029382 0.2380952 -0.1718791 -0.0442981

The biquad coefficients are taken from http://jiiteepee.tripod.com/de-pre-emphasis.html. I wonder if there is any easy way to calculate the biquad coefficients if I want to apply pre-emphasis only partially - in the range from 0% (no pre-emphasis) to 100% (full pre-emphasis).

The case for 0% seems to be easy, coefficients 1 0 0 1 0 0 would work just fine. But I'm not sure how to approach calculation of biquad coefficients for the combination of zero_filter + float_fraction * full_filter.

So, sox defines its biquad operation parameters as follows:

biquad b0 b1 b2 a0 a1 a2

  Apply a biquad IIR filter with the given coefficients. Where b* and a* are the numerator and denominator coefficients respectively.

  See http://en.wikipedia.org/wiki/Digital_biquad_filter (where a0 = 1).

  This effect supports the −−plot global option.

A biquad filter is just a special case of an IIR filter, which is a linear system!

That means that you can just have two signal streams, one filtered and one un-filtered, and add them up, scaled by two coefficients with constant sum.

The only thing you'd have to care about is delay, which is a bit tricky to calculate for IIRs; now, luckily, the RIAA biquad guarantees phase errors of this filter to be below 30% at any given frequency, so you can probably get away with summing up with a single-sample delayed input signal.

A single sample delay is a FIR filter with coefficients

b0 = 0, b1 = 1

and any FIR is also understandable as IIR with all aN==0 but a0=1, since IIRs are linear, you can directly modify the coefficients of your biquad by mutliplying them with the aforementioned coefficient and adding 1-coefficient to a0 and b1


How to extract features of individual spikes in non-negative signal? I have data of non-negative (in the sense there's no signal below baseline) spiking waveforms, which are in the form of a 1D array of numbers:

Spikes that cross some threshold are considered real signal events and I want to be able to measure the properties of individual events such as their amplitude, rise-time, decay-time, duration, and half-width.

The signals are pretty clean so I imagine I don't need some fancy machine learning algorithm. Would a wavelet transform approach give me the results I need? If not, what's the best and simplest approach here?

Continuous wavelet transforms can provide information and parameters for sparse, piece-wise regular signals. One example related to yours is present in The Continuous Wavelet Transform in MRS, A. Suvichakorn, C. Lemke, A.  Schuck, J.-P. Antoine, as exemplified in the following picture.


Telling apart a right traveling non-periodic signal from a left-traveling non-periodic signal, with a background periodic signal I have a time-signal sampled over a one-dimensional spatial domain $x$. I know it is made of 3 components:

a right traveling non-periodic wave (characteristic time width $T_{R,NP}$, mostly 0 outside this width). This is prescribed at the left end of the domain.
a left-traveling non-periodic wave (characteristic time width $T_{L,NP}$, mostly 0 outside this width). This is prescribed at the right end of the domain
a periodic signal (right traveling in my case, of period $T_P$)

While traveling, all these signals slightly change shape due to non-linearities. Again changes are not drastic and fairly small. The periodic signal starts as a single harmonic but then more harmonics show up, which have small energy and the peak period is still $T_P$. Here is an example of the signal at a specific location (purple), together with my best guesses (I drew them by hand) of the two non-periodic signals. Two signals were prescribed on the left and traveled right, and one was prescribed on the right and travel left.

Which is the best way to tell apart the 3 components at different locations of my domain? In particular:

Is there a method that can take advantage of the known information about the direction of travel?
Sometimes (but not always) $T_{L,NP} > T_{R,NP} > T_P$. Would this help? I would prefer a general algorithm that does not rely on this information, but if a general algorithm does not exist, solving for these special cases would still be something.

If you can sample at any point in time, this shouldn't be that difficult.  The first step is to collect as series of frames in time.
In each frame you should be able to find a best fit sinusoid, or set of harmonic sinusoids, to your periodic signal.  Since this signal is assumed to have the same parameters across time, except the phase values which would need to be time adjusted, you can find an "average" value across all your frames for the periodic signal.
Subtract the best fit periodic signal from your frames.  This should leave the other two behind.  From your drawing, it appears that you have regions where one hump will b e zero, and others where the other hump is zero.  Those are the locations to read your humps at.
For more accuracy, once you have your hump signals estimated, go back and subtract them from your frames.  This should leave the periodic signal behind.  Once again, find the best fit sinusoid assembly.  This time your results should be more accurate.
Rinse and repeat.


What transformations have been applied to this image to yield four derivatives? In a few days I have an exam of image processing. I'm currently having a look at the example exam and I came across this question. "Which transformation of this image is being done here. Explain the components in the quadrants NA (second image)."

Na

Can anyone help me with this ?

This to me looks exactly like the two-dimensional Haar wavelet transform. All four quadrants are a part of the transform, even though different operations have been performed on them. Upper left quadrant is a result of low-pass filtering with a 2x2 two-dimensional rect function (Haar wavelet) and decimation be a factor of 2, upper right is the output of a horizontal high-pass filter, lower left is the output of a vertical high-pass filter, and lower right is the output of a diagonal high-pass filter.


Hessian Matrix. Second partial derivative test I have a question about the hessian matrix / determinant of the hessian matrix. For the Hessian matrix, one is able to determine the local curvature around a point. And with the Second partial derivative test one can test, using the determinant of the hessian, if that point is located on an extrema, or saddle point.

My question is
Wikipedia states:

  the Hessian approximates the function at a critical point with a
  second degree polynomial.

By 'hessian' does this mean that the hessian matrix approximates the function or the determinant? (i am confused as determinant of hessian and hessian matrix are sometimes used intertwined)

thank you in advance

The function is fully approximated if one uses all the derivatives (see Taylor expansion). With using the Hessian only, we can only make a second degree approximation (because it is second derivative matrix), which is geometrically the same as using the second order polynomial.


Needed: introduction to FFT I am looking for a great introduction to FFTs for my job. I studied math, but never heard of FFT until I got this job. There is not a lot of documentation at my place of work and the articles I have found online have been somewhat confusing, even the wikipedia entry.

So I am looking for a recomendation for a higher level introduction to FFT. I would perfer free resorces over paid ones, but I can't be too picky here.

Thanks in advance!

I found the following review paper very illuminating:

Fast Fourier Transforms: A Tutorial Review and a State of the Art
by P. Duhamel and M. Vetterli.

It covers a lot of things, from history to implementations.


FIR filter Gain of 0dB at Passband Frequency I have designed a FIR pass band filter. The impulse response looks good and on applying this on my input data, the filter works as desired. Now I have a requirement in my system to have exact $0$ dB gain at passband frequency and hence I normalize the FIR coefficients using the following method:

for (n=0; n < N; n++)
  Gain += FIRCoefficients[n] * cos(2 * PI_DOUBLE * PassbandFreq * n);

for (n=0; n < N; n++)
   FIRCoefficients[n] = FIRCoefficients[n] / Gain;

Even after this normalizing, the Gain at passband frequency is around $23$ dB. What am I doing wrong?

FIR coefficients are $h[n]$, the same as the impulse response.  there are $N$ non-zero taps.

$$ h[n] = 0 \qquad \text{for } n<0, n \ge N $$

frequency response is

$$\begin{align}
H(e^{j \omega}) &= \sum\limits_{n=-\infty}^{+\infty} h[n] e^{-j \omega n} \\
                &= \sum\limits_{n=0}^{N-1} h[n] e^{-j \omega n} \\
                &= \sum\limits_{n=0}^{N-1} h[n] \big( \cos(\omega n) - j \sin(\omega n) \big) \\
\end{align}$$

the magnitude-squared is

$$\begin{align}
\Big| H(e^{j \omega}) \Big|^2 &= \left| \sum\limits_{n=0}^{N-1} h[n] e^{-j \omega n} \right|^2 \\
                &= \left| \sum\limits_{n=0}^{N-1} h[n] \big( \cos(\omega n) - j \sin(\omega n) \big) \right|^2 \\
                &= \left| \sum\limits_{n=0}^{N-1} h[n] \cos(\omega n) - j \sum\limits_{n=0}^{N-1} h[n] \sin(\omega n) \right|^2 \\
                &= \left( \sum\limits_{n=0}^{N-1} h[n] \cos(\omega n)  \right)^2 + \left( \sum\limits_{n=0}^{N-1} h[n] \sin(\omega n) \right)^2 \\
\end{align}$$

you must square root that expression to get magnitude of gain

$$ \Big| H(e^{j \omega}) \Big| = \sqrt{\left( \sum\limits_{n=0}^{N-1} h[n] \cos(\omega n)  \right)^2 + \left( \sum\limits_{n=0}^{N-1} h[n] \sin(\omega n) \right)^2 }  $$

evaluate at at your known PassbandFreq $\omega_0$ (normalized angular frequency) and divide each $h[n]$ by that gain magnitude:

$$ h[n] \leftarrow \frac{h[n]}{\ \big| H(e^{j \omega_0}) \big| \ } $$


Fourier Series of Aperiodic convolution of periodic functions we were given the following classic exercise:

  Given two periodic signals $x(t), y(t)$ with fundamental period $T$
  with Fourier series coefficients $c_m^x, c_m^y$ respectively, find the
  Fourier coefficients of the signal $z(t) = x(t) * y(t)$ with relation
  to $T, c_m^x, c_m^y$.

Now, this can easily be solved when the aforementioned convolution is the circular convolution (integral over a period only).

However, in class our professor noted that it can be solved even when we have an aperiodic convolution (that is, convolution as an integral from $-\infty$ to $+\infty$). We argued that, in this case, that infinite integral doesn't converge, and he responded that, even though the convolution integral doesn't converge (i.e. might be infinite), the Fourier Series coefficients are still finite and can be calculated!!

Is this true? If yes, then is the relation the usual one: $c_m^z=Tc_m^xc_m^y$ or another and how do you prove that? If not, why? Austere mathematical proofs would be appreciated.

okay, now that the popcorn is eaten and my fingers are clean...

$$ x(t) = \sum\limits_{m=-\infty}^{+\infty} c_m^x e^{j 2 \pi (m/T) t} $$

$$ y(t) = \sum\limits_{m=-\infty}^{+\infty} c_m^y e^{j 2 \pi (m/T) t} $$

then

$$\begin{align}
 z(t) &= \sum\limits_{m=-\infty}^{+\infty} c_m^z e^{j 2 \pi (m/T) t} \\
\\
  &= x(t) \circledast y(t) \\
\\
  &= \int\limits_{-\infty}^{+\infty} x(\tau) y(t-\tau) \, d\tau \\
\\
  &= \int\limits_{-\infty}^{+\infty} \sum\limits_{m=-\infty}^{+\infty} c_m^x e^{j 2 \pi (m/T) \tau} \sum\limits_{n=-\infty}^{+\infty} c_n^y e^{j 2 \pi (n/T) (t-\tau)}  \, d\tau \\
\\
  &= \sum\limits_{n=-\infty}^{+\infty} \sum\limits_{m=-\infty}^{+\infty} c_m^x c_n^y  \int\limits_{-\infty}^{+\infty} e^{j 2 \pi (m/T) \tau} e^{j 2 \pi (n/T) (t-\tau)}  \, d\tau \\
\\
  &= \sum\limits_{n=-\infty}^{+\infty} \sum\limits_{m=-\infty}^{+\infty} c_m^x c_n^y \int\limits_{-\infty}^{+\infty} e^{j 2 \pi ((m-n)/T) \tau} \, d\tau \ e^{j 2 \pi (n/T) t} \\
\end{align} $$

this integral: $\int\limits_{-\infty}^{+\infty} e^{j 2 \pi ((m-n)/T) \tau} \, d\tau$ ain't converging when $n \ne m$ and certainly not when $n=m$.

you will not be getting finite values for $c_m^z$.


Are all isotropic bandpass filter kernels, wavelets? Are all isotropic bandpass filter kernels, wavelets?

Are there any properties that the radial frequency profile has to meet to be considered a wavelet?

Certainly not. ,,Isotropic'' means that the impulse response is independent of direction, that is, that the filter is rotationally symmetric. Wavelets are orthogonal basis functions (filters, maybe) with certain properties. Not every orthogonal basis function will be isotropic, and not every isotropic filter is orthogonal. It is even hard to see how a isotropic filter may be a wavelet.

Look up ,,wavelet conditions'' to learn more about how this orthogonality has to be understood to build an invertible transform. Banded Matrices with Banded Inverses might help the algebraically inclined.


Apply the Hamming window to a sinusoidal signal I am a novice in signal processing. I am writing a code to apply the hamming window to a sinusoidal signal. So first I have digitized my sinusoidal signal by choosing a sampling rate more than twice its frequency.

These are the values for the signal: F=500 Hz, T=2 ms, and sampling rate=2000 Hz.

Now, I guess I have to multiply the function of the window to the signal so:

for (i=0; i<windowSize; i++)
     data(i)=signal(i)*window(i)
     i=i+0.5;  //sampling time

My question is when we multiply the window to the signal, does the window function has to be digitized by the sampling rate of the signal or it can be digitized anyway we want? In the code I have digitized the window function according to the sampling rate. Thank you.

The window isn't "sampled" as it isn't a signal, really.  It's a vector that you construct, and it doesn't matter how you construct it.  The window vector just has to be the same length as your signal vector.


Logarithmic Scale Question I just have a small theoretical question that buggers me to be honest.

I have a sound file that I found its spectrum on linear and on logarithmic scale. Alright, I have done this, it's a plain standard procedure. The problem is that I cant really understand the plot on logarithmic scale. (it just looks like a typical logarithmic plot with dB, and it is correct by the way) .

Well, the visual differences are obvious, the linear plot frequencies are quite scattered and they peak in a non homogenous way, while I notice that on log scale the frequencies peaks are very close to each other and they don't seem to have so many extreme values. Could someone explain me this a bit ?

Cheers

The meaning of the logarithmic plot is that the signal power above 250 Hz is about ${10^{-10}}$ times less that the starting values (most of the power is concentrated at the lowest frequencies), so it's almost zero. The capability of a logarithmic plot is to show both, the high and low values at the same time, specially useful in acoustic field. The negative aspect is the visual distortion of the signal, you must always keep in mind that what are you looking at is a logarithmic plot.


How does projecting a few IR dots capture an image three dimensionally? The original Microsoft Kinect device IR projected a spacially encoded dot field and looked for distortion in the dot field with a single IR camera to create a depth map. But there were thousands of dots.

I've seen a few devices however, usually in motion capture, that only use a few dots. The one pictued below for example is used to determine the position, movement and rotation of a moving golf ball.

Do they use a similar process of looking for distortion of the projected pattern, or some other process?

The projected differaction pattern (formed out of Gaussians) has a special assymetric structure, which in turn would help the software resolve the stereo correspondence problem. After correspondence is obtained, the distance can be triangulated, as in a simple stereo setting.

For further info, please check the PrimeSense patent:
http://www.google.com/patents/US20090185274

Other stereo reconstruction techniques might use color coded patterns, blue light patterns, phase shifting projections and so on. Generally, tools projecting fewer dots enjoy muliple pattern projection and encode every pattern differently. This might for example obtained by rotating the projector or programming it to project different light stripes.

Depending on the object to be reconstructed, using multiple cameras might just be enough as the RGB image information might provide you enough data to disambiguate. In this case, you wouldn't even have to project a pattern.

The technique is generally called structured light and for further info you might check:
http://fofi.pagesperso-orange.fr/Downloads/Fofi_EI2004.pdf
http://eia.udg.es/~qsalvi/Tutorial_Coded_Light_Projection_Techniques_archivos/v3_document.html
http://en.wikipedia.org/wiki/Structured_light


Removing a non linear noise trend from data I have a dataset that looks as follows

Now, this is from an actual measurement (of the reflection coefficient of a microchip with four LC resonators connected in a periodic fashion) so there is obviously noise involved. For reference, this is what it should look like in theory (with different numerical parameters)

So there's a few things about the data. First of all there is gaussian(?) noise around the 'mean'. This is not really an issue; I'm only really interested in the positions/sharpness/depths of the resonances, and the noise doesn't complicate it too much. I suppose that if it would be very problematic, one could use Gaussian filter?

But then there's a more problematic issue, and that is also why I put accents around mean. The data is not in a straight line; there is some nonlinear trend over the range of the data. This is what my question is about; getting rid of this trend so that the data is on a straight line. Once that is done, getting it to have its mean at 1 is of course easy, and I'll have a dataset that can be fitted with the model.

My issue is that I don't really know where to start. I did some searching, but either my search terms are terrible or this hasn't really been discussed; it is quite specific after all. I'd like your advice on how to proceed, as personally I don't have any clue. Maybe I should fit the data with some general function? It does seem to be somewhat oscillatory, but not very nicely.

I apologize for not providing some code to generate similar looking data. I would, if I had an idea of how to go about that, but I don't as I'm not sure what kind of function would generate this type of trend. Instead I have uploaded the x-coordinates (Frequency in GHz) and the y-coordinates (complex reflection coefficient) here. As commented, this is a rather strange format for the data if you are not using Mathematica. Here is a different version which should work anywhere, and/but already has only the absolute value of the y-coordinates (as mathematica's way of storing complex numbers is a little nonstandard)

(Dropbox, should be safe and registration free to download it). I should definitely add that what is plotted above is the absolute value of the y-coordinates.

Note that it doesn't require registration. It'll pop up a registration dialog when you open the site, but you can close this, and then you can click the download button (although it seems greyed out, but it isn't). They write this here too dropbox.com/help/20

As for the source of the noise. This is one dataset of many, taken with a Vector Network Analyser. It's connected to the device with SMP cables, which are calibrated before doing the measurement. The calibration is supposed to take care of these trends, but in this case it didn't. So, does theory tell us something.. I'm not sure. I think the noise comes from the cables, and that it might have some oscillatory component to it, but beyond that I can't say to be honest, I have almost no experience with the system. The reason I don't just remeasure is that it is done at liquid helium temperatures, so I'd rather post process

The trend can be perceived as some low frequency signal, so a simple highpass filter should remove it. It will naturally slightly decrease the peaks, but by very little, because such narrow peaks are mainly composed of high frequencies.

Here is some concept code in Python (using the NumPy, SciPy and Matplotlib libraries)

import numpy as np
import scipy.signal as scisig
import matplotlib.pyplot as plt
# assumes it has just two columns with f, r values
f, data = np.loadtxt('CQDataAbs.dat', unpack = True)
# order4 IIR Butterworth highpass filter coefficients,
# the critical frequency is relative to the Nyquist frequency,
# determining the right value might require some experimentation
# note that using a higher order with such low cuttof will likely
# be numerically unstable, if you need higher orders,
# split it into several lower-order stages applied sequentially
b, a = scisig.butter(4, 0.0003, btype='highpass')
# forward-backwards filter application with edge padding (mitigates edge defects),
# results in zero phase lag and an effective filter order 8 as it is applied twice
filtered = scisig.filtfilt(b, a, data)
plt.plot(f, filtered, 'b.')
plt.show()

And here is the result

The shoot-up at the end may indicate that there are some higher frequencies in that area. Perhaps you can just disregard that region during fitting if possible.

If you were interested in the trend itself, you would obtain it by using a lowpass filter instead of a highpass filter. With zero phase lag application the result of highpass filtering should be almost identical (within numerical accuracy) to subtracting this lowpass trend from the original data.

Note that you may find that a different filter type might give you better results. A Butterworth IIR filter (with forward-backwards application to remove phase lag) is usually a good choice when the shape of the signal is important, because the frequency response within the pass-band is maximally flat. You could of course use some large FIR filter for better performance if you have enough processing power. With FIR filters you also don't have to worry about numerical instability that can occur in IIR filters of high orders.

Just in case you think about using a moving average (which is a FIR lowpass filter): That will likely perform terribly in this case, because in each window containing a peak the average will shoot up and subtracting that from the original data will greatly decrease the peaks. From a frequency perspective, the moving average filter does not mitigate higher frequencies very well, which are very significant within the region of the peaks.


Finding audio "tics" in stereo wav files using python I have thousands of audio files that I want to "batch" process searching for audio tics.  This is easily done using a DAW but that requires opening each file independently, searching for the tic (visually looking at the waveform or listening to the file play), and then cleaning the tic with signal processing tools within the DAW.  However, this approach is extremely time consuming.

I'd like to write a program using python to process all of the files and provide candidate files that contain tics in the audio.  I've been researching the web and found tons of references for plotting frequencies using fft after reading in the wav file but I don't want to plot the data I want to examine it in the program to determine if there is an audio tic (or potential tic) present.

Any help or ideas would be deeply appreciated.

I don’t know what DAW you use, but I’m pretty sure Waves and probably some other companies have plug ins for this kind of application.
After reviewing your images it seems that the tics are high frequency and high amplitude so there are some simple things you could try.
You could check each sample and compare the absolute value to some threshold, such that if it is higher, you have a tic.
Second, you could high pass filter the signal to isolate the tic, then do the previous step.
Third, you could calculate the RMS of the signal, then compare the absolute value (peak) of each sample to the RMS, such that a tic occurs whenever the ratio or difference exceeds a given threshold.
Lastly, you could do frequency magnitude analysis.  For this you would set up a window size N and step size M.  Every M samples you take N samples and apply a window function, then calculate the frequency transform magnitude of the N samples.  One idea would be to look at the rate of change of individual frequency components such that if it exceeds a threshold you have a tic.  This is probably the most CPU intensive and complex implementation.


Why is the raised cosine pulse for digital communications so popular? The raised cosine pulse is essentially a windowed sinc function, designed with some roll-off in mind. Why do we use this particular window, as opposed to other windows such as Hann, Hamming, Kaiser, etc? All of these have zero ISI when applied to a sinc function, and have tunable frequency characteristics just like the raised cosine. Is the raised cosine optimal for digital comms in any sense?

There is no good reason, except that it's familiar and "good enough" for many applications. Many other better pulses have been proposed; see for instance "Nyquist filters in non-ISI transmission", by Ping-Kuen Lam, E.W. McCune and M.A. Soderstrand, DOI:10.1109/MWSCAS.1997.666198.


Is it possible to apply inverse Fourier transform (I-FFT) to the following image? Suppose, I have this image in my hand and nothing else. Clearly, this is a Magnitude-plot of some unknown image.

Is it possible to apply an Inverse Fast Fourier Transform (I-FFT) operation to recover the original Grayscale image from this image?

N.B. For more information and actual C# source code, please follow this link.

Suppose, I have this image in my hand and nothing else. Clearly, this is a Magnitude-plot of some unknown image.

Sounds like you're losing the phase information in the signal; you only have the pixel values, correct?

  Is it possible to apply an Inverse Fast Fourier Transform (I-FFT) operation to recover the original Grayscale image from this image?

The important information in images tends to edges -- edges contain most detail in the image. This information is carried in the phase information -- when decomposed into sinusoidal signals, a square wave can have arbitrary location only if the waves contain phase information.

Additionally, real world images tend to look roughly the same in the frequency spectrum when only plotting magnitude. These images tend to have a large (in magnitude) low frequency components with (much) smaller high frequency components. It's reasonable to say that a significant portion of the entropy (or information) is carried in the phase information.

With a naïve inverse Fourier transform on the values obtained from the image, it is not possible (at least by experiment) to recover the original signal. However, the phase only inverse Fourier transform (where we divide by the magnitude to only preserve the phase).

The data behind the image was generated with

x = rgb2gray(astronaut())
X = fft2(x)
X_mag = np.abs(X)  # magnitude only
X_pha = X / (np.abs(X)) # phase only
y_mag = ifft2(X_mag)  # both y_mag and y_pha have max(y.imag) < 1e-15
y_pha = ifft2(X_pha)  # (some Fourier transform theory here + floating point precision)