Design 2nd order sections IIR filter for high orders (more than 500) I am currently building a Reverb that needs to fit on a DSP. I am using the model proposed in this link https://ccrma.stanford.edu/~jos/pasp/Freeverb.html.
It involved to create some Feedback Comb-Filters with high orders (more than 1000). For that I use a 2nd-order section IIR filter (IIRSOS) integrated in my DSP. I am using the Matlab function "tf2sos" that takes numerator and denominator of the Transfer Function (TF) of the desire filter and create a matrix of Second Order Section (SOS) coefficient. With this coefficients, I can make what I want.
Problem : for high filter orders (arround order 200) the "tf2sos" function appears to be totally unstable. However the Matlab object called dsp.IIRFilter that use also 2nd-order section IIR seems to handle it :
Impulse response for order = 200 :

Impulse response for order = 1500 :

On the 2nd picture you can see that the filter with the tf2sos coefficients is going totally wrong. Does someone knows how I can get a functional coefficient matrix for orders like 1500 ?
Thank you.

The Schroeder/Moore implementation uses sparse IIR filters. It's extremely inefficient to implement these as regular IIR filters or biqauds.

It's much easier and more efficient to implement these directly. See for example
https://ccrma.stanford.edu/~jos/pasp/Schroeder_Allpass_Sections.html

Both the comb filter and the schroeder allpass use the same structure as a regular IIR filter, but the delay(s) are not a single sample but much larger number. So the Schroeder allpass looks like any other 1rst order allpass except that the delay in there is big. If you want get technical about this

The algorithm was specifically designed to be computational efficient, hence it makes no sense using standard biquads. On the downside, it typically doesn't sound very good. Feedback delay networks tend to sound a lot better.


Why we have to multiply the images with $(-1)^{x+y}$ to center the transform in frequency domain filtering? To perform the filtering infrequency domain we perform multiplication of $(-1)^{x+y}$ why??

Below 1D argumentation also explains the 2D case.

First consider the DTFT property for the pair $x[n] \longleftrightarrow X(e^{j\omega})$

$$ e^{j\omega_0 n} \cdot x[n] \longleftrightarrow X(e^{j(\omega - \omega_0)}) $$

Then recognise that $(-1)^n = e^{j \pi n} $ which yields:

$$ e^{j\pi n} \cdot x[n] \longleftrightarrow X(e^{j(\omega - \pi)}) $$

The effect is such that the spectrum $X(e^{j \omega})$ is centralized (shifted) into the $\omega = \pi$ frequency.

In 2D, the DTFT of the image is shifted into the central zone.


What are some techniques to recover an original image from three different compressed versions? This question is theoretical. Let's say you have three different highly compressed images of a same source image, each compressed with a different algorithm and thus "lossy" in different ways.

What would be the possible techniques to recover the best possible image, using the three different compressed versions ?

First of all, there are many solutions. Something better under one measure may be worse under another measure. So first you need to think about what is your quality measure(s) to evaluate a restored image.

Commonly, people use the "mean square error (MSE)" (or its log version, known as the peak signal noise ratio). Assume the ground-truth image is $X$, and the recovered image is $\hat{X}$. Then MSE is defined as
$${\rm{MSE}}(X,\hat{X}) = \sum_{i,j}(X[i,j]-\hat{X}[i,j])^2/N$$
where $N$ is the number of pixels in $X$, and $X[i,j]$ denotes a pixels located at $i$th row and $j$th column in $X$.

Without loss of generality, say we have $n$ compressed counterparts are $C_1,C_2,\cdots,C_n$, and we reconstruct $\hat{X}$ as
$$\hat{X}=f(C_1,C_2,\cdots,C_n)$$
where $f(\cdot)$ is a method to recover the original image from its compressed counterparts. And your objective is to find an appropriate function $f(\cdot)$ and minimizes the MSE.

To simplify the problem, people often restrict the function $f(\cdot)$ to some fixed type but with parameters to be learned from training data. For example, one simplest function can be the linear combination i.e.
$$f(C_1,C_2,\cdots,C_n) = \sum_{k}w_kC_k$$
where $w_k$ is the weight of using the $k$th compressed image in the recovered image. This is a classic regression problem and can be solved easily (go and check the normal equation, which gives the closed form solution to find all $w_k$s). However, you should note the following things:

1). if you donot have any training data, the above method does not
work
2). if your image is very large or your have many training images,
computing the closed form solution is not a good idea, because of too
many computations. Instead, you shall seek for methods like the
stochastic gradient descent or similar ones.
3). if I were you, I will go and read more about linear regression
stuff. You will see interesting topics like how to prevent overfit by
adding regularization, and how to include the Bayesian idea in this
regression problem.
4). my guess says the linear regression method is not the best way,
but it is a very good start point. Later your may try other
$f(\cdot)$, like kernel regressions, or adding more terms, e.g. the
total variation term.


power spectral density plots Hi all,

How do we interpret a power spectral density plot. I have used modified covariance and burg methods and plotted using MATLAB. What does the peaks at some frequencies suggest? and what is the difference between positive and negative power in the graph?  how do we know that which frequencies variations are strong and at which frequencies variations are weak.. can anyone help me please..

Thanks

The plot is not of power.  The units of the $y$ axis are dB/Hz. The dB part means that a decibel scale is used, which means it is a scale relative to a reference value.

All that negative values on a decibel scale mean is that the actual power is less than the reference power (because the ratio is less than one).

The higher the decibel value, the more powerful the signal.

Your signal has a major peak at about 4Hz, and lesser peaks at 28Hz and 39Hz.  There might be something happening at about 12Hz, but it's hard to say without more information (or data).


A PID that can handle 2 inputs for my application I am currently working on a Drone that navigates to different points that it detects on the ground. To go from one point to another in a smooth manner, I am employing a PID looks at the error between the current and the desired location of the drone. Currently this is a two dimensional problem. Soon enough, it is going to be a three dimensional problem as I will try and have the drone navigate three dimensional spaces. Now I am unsure as to how to program a PID that can take in more than one input and give an two outputs (which are x and y in this case for the Cartesian geometry the drone is trying to navigate).

Can someone take me through the exact theory behind a PID that approaches the 2D problem and different solutions that are possible. I also want to know about any solutions that are present in C++ that can help me speed up creating this program (other languages are also okay, but I am mostly building this project in C++).

If you require two outputs, which belong to a movement in x and the other in y direction, you could start with the assumption, that both movements are independent an simply implement two PID controller, each caring for one direction.

It gets interesting at the point where you translate the output of these two controller to signals to the actuators. Here it may be wrong to assume that both movements are independent. If you give more insight on that, we can give better ideas for control algorithms

(more a comment then an answer, but I cant comment)


Implementing a Kalman filter for position, velocity, acceleration I've used Kalman filters for various things in the past, but I'm now interested in using one to track position, speed and acceleration in the context of tracking position for smartphone apps.  It strikes me that this should be a text book example of a simple linear Kalman filter, but I can't seem to find any online links which discuss this.  I can think of various ways of doing this, but rather than researching it from scratch, perhaps someone here can point me in the right direction:Does anyone know the best way of setting up this system?  For example, given the recent history of position observations, what's the best way of predicting the next point in the Kalman filter state space?  What are the advantages and disadvantages of including acceleration in the state space?  If all measurements are position, then if speed and acceleration are in the state space can the system become unstable?  Etc ...Alternatively, does anyone know of a good reference for this application of Kalman filters?Thanks

This is the best one that I know of

Full derivation with explanation

Kalman

This is a good resource for learning about the Kalman filter. If you are more concerned with getting the smartphone app working I would suggest looking for a pre-existing implementation of the Kalman filter. Why reinvent the wheel? For example if you are developing for android, openCV has an implementation of the Kalman filter. See Android OpenCV

Bradski and Kaehler is a good resource on image processing in general and includes a section on the Kalman filter including code examples.


Find similarity of sound signal So I was working on a program to comparing the similarity between the pronunciation of a word, First I got the audio file, and I use FFT on it to get the amplitude over frequency array, then I just square everything to get the Energy Spectral Density of it. The graph that can be drawn from those frequency of 4 sample of the same word I spoke myself look very similar to each other as shown in the picture I included here. All 4 sample have 2000 value, represent frequency from 0 to ~800Hz.

But since even when it's me saying the word 4 times, the frequency of the word cant be the same everytime I say that word, the graphs look visually similar but they are a few sample different from each other, like the 2nd highest point in the first graph is at the value number 1150, the 2nd graph have it at 1124, the 3rd have it at 1132 and the 4th have it at 1160, same thing happen with the first peak.

So my question is: how do I actually compare these 4 graph together to show how similar they are? Since what I was suggested to do the first time around was to using Mean Square Error but they didnt work well due to the slight different everywhere.

To compare the similarity between multiple pronunciations of a single word, the spectral power that you've computed is too general a feature and doesn't capture details of the shape of the vocal tract, where variations in pronunciation may come from.

I recommend you compute Mel Frequency Cepstral Coefficents (MFCCs), which is widely used in speech recognition systems because it captures the envelope of the short time power spectrum of the vocal tract. You can build on your FFT code, but there are also many good toolkits available.

Once you have MFCCs sequences for each pronunciation you can use Dynamic Time Warping to compute the corresponding "distance" or "similarity" between them.


Find the Differential Equation of an Analog Band Pass Filter I don't even know how to start. We're supposed to find the differential equations that links $V_{in}$, $V_{out}$, and $V_{io}$, as well as the values of $C_1$, $C_2$ and $R$ so the filter works in the audible spectrum. We're allowed to use Fourier transform, if that simplifies anything.

I thought of splitting the band pass filter as a high pass at first, and a low pass after.

That way I can relate $V_{io}$ with $V_{in}$, and $V_{out}$, with $V_{io}$, in two different equations.

Doing that... the results were:

$RC \frac{\mathrm{d} Vin(t)}{\mathrm{d} t} = Vio(t) + RC \frac{\mathrm{d} Vio(t)}{\mathrm{d} t}$

for the high pass filter

$Vio(t) = Vout(t) + RC \frac{\mathrm{d} Vout(t)}{\mathrm{d} t}$

for the low pass filter

Then, using..
$fc = \frac{1}{2 \pi R C}$

For the audible frequency, it'd be between $20$ hz- $20$ KHZ. So using that, I should be able to express both $C_1$ and $C_2$ as a function of $R$.

However, how do I actually determinate the value of $R$ required? or do I assume a set value?

If there's anything wrong, please let me know why. Thank you, have a nice day.

I'm afraid your first equation is wrong because it assumes a highpass filter without any load. However, there is a lowpass filter connected at its output and for this reason the current through the resistor leading to ground does not equal the current through $C_1$.

What you need to do is apply the complex voltage divider rule in the frequency domain, which will give you

$$\frac{V_{io}(\omega)}{V_{in}(\omega)}=\frac{Z(\omega)}{1/j\omega C_1+Z(\omega)}\tag{1}$$

where $Z(\omega)$ is the impedance of $R$ in parallel with the series connection of $R$ and $C_2$. I'm sure you can figure this out. The relation between $V_{out}$ and $V_{io}$ is simply

$$\frac{V_{out}(\omega)}{V_{io}(\omega)}=\frac{1/j\omega C_2}{1/j\omega C_2+R}=\frac{1}{1+j\omega RC_2}\tag{2}$$

The total response $V_{out}/V_{in}$ is simply obtained by multiplying Eq. (1) with Eq. (2). The time domain differential equations are obtained by multiplying Eqs (1) and (2) by the denominators of both sides and noting that a multiplication by $j\omega$ corresponds to differentiation.


Remove minor cycle from time series In the above image I have created a time series (in cyan) that is a simple addition of a major cycle (in green), a minor cycle (in cyan) and a linear trend (not shown separately). What I would like to do is filter this time series to remove the minor cycle so what is left would be the trend plus major cycle (in red). In practice no information would be available about the minor cycle, although the period of the major cycle would be known from applying code for the Hilbert transform.

One way in which I might do this is:

apply a bandpass filter set to the major cycle period, on the original data
subtract this bandpass from the original series
apply the Hilbert transform to the remainder to get period of minor cycle
apply another bandpass filter set to this minor cycle period, on the original data
subtract this second, minor cycle bandpass from the original data
and hopefully end up with the desired red output

First off, would this be a valid approach, and secondly is there another way in which I could filter out the minor trend?

Your approach seems awfully complicated. It looks like the frequencies of the "minor cycle" are much higher than the "major cycle" and the ramp so a simple lowpass filter should do the trick. If you don't know the frequencies, a simple FFT analysis can be done.

Subtracting filtered versions of signal is tricky since you need to carefully watch time alignment and phase distortion in the filter.


Nulling filter coefficients I'm wondering if I have calculated a nulling filter coefficients correctly:

$\ 1, -2\cos(0.44\pi), 1$

I have to make a nulling filter that filters out the frequencies $\theta = 0.44\pi$. I've been told that the formula for nulling filters is $y[n] = x[n] - 2\cos(\theta)*x[n-1]+x[n-2]$.

If I'm correct the two 1s come from the fact that $\ x[n]$ and $\ x[n-2]$ both have a coefficient of 1 and no other variable and the $\ -2\cos(0.44\pi)$ comes from the $\ - 2\cos(\theta)*x[n-1]$ and the fact that $\theta$ = 0.44$\pi$.

If the coefficients are correct, is all I have to do to find the output signal just multiply the input signal by each coefficient and add them together (convolution sum)?

You are looking at an FIR filter whose transfer function is

$$
\begin{align}
H(z)&=(1-e^{j\omega_0}z^{-1})(1-e^{-j\omega_0}z^{-1})\\[10pt]
&=1-(e^{-j\omega_0}+e^{j\omega_0})z^{-1}+z^{-2}\\[10pt]
&=1-2\cos(\omega_0)z^{-1}+z^{-2}
\end{align}$$

So as you can also see in the comment, there are two zeros at $e^{\pm j\omega_0}$ that null the frequencies $\pm\omega_0$.

Obviously, here we have $\omega_0=0.44\pi$. The filter's response is

The impulse response is $$h[n]=\delta[n]-2\cos(0.44\pi)\delta[n-1]+\delta[n-2]$$

Assume an arbitrary input $x[n]$. The output $y[n]$ is given by:

$$y[n]=x[n]*h[n]$$ Therefore $$y[n]=x[n]-2\cos(0.44\pi)x[n-1]+x[n-2]$$


Pink ($1/f$) pseudo-random noise generation What are some algorithms for generating a good pseudo-random approximation to $1/f$ (pink) noise, yet suitable for implementation with low computational cost on an integer DSP?

Linear Filtering

The first approach in Peter's answer (i.e. filtering white noise) is a very straightforward approach. In Spectral Audio Signal Processing, JOS gives a low-order filter that can be used to produce a decent approximation, along with an analysis of how well the resulting power spectral density matches the ideal.  Linear filtering will always yield an approximation, but that may not matter in practice.  To paraphrase JOS:

  There is no exact (rational, finite-order) filter which can produce
  pink noise from white noise. This is because the ideal amplitude
  response of the filter must be proportional to the irrational function
  $1/\sqrt{f}$, where $f$ denotes frequency in Hz. However, it is easy
  enough to generate pink noise to any desired degree of approximation,
  including perceptually exact.

The coefficients of the filter he gives are as follows:

B = [0.049922035, -0.095993537, 0.050612699, -0.004408786];
A = [1, -2.494956002, 2.017265875, -0.522189400];

They're formatted as parameters to the MATLAB filter function, so for the sake of clarity, they correspond to the following transfer function:

$$
H(z) = { .041 - .096z^{-1} + .051z^{-2} - .004z^{-3} \over{}  1 -2.495z^{-1} + 2.017z^{-2} - .522z^{-3}}
$$

Obviously, it's better to use the full precision of the coefficients in practice. Here's a link to what the pink noise generated using that filter sounds like:

For fixed-point implementation, since it's usually more convenient to work with coefficients in the range [-1,1), some reworking of the transfer function will be in order.  Generally, the recommendation is to break things up into second-order sections, but part of the reason for that (as opposed to using first-order sections) is for the convenience of working with real coefficients when the roots are complex.  For this particular filter, all the roots are real, and combining then into second-order sections would probably still yield some denominator coefficients >1, so three first-order sections is a reasonable choice, as follows:

$$
H(z) = { 1 - b_1z^{-1}\over{1 - a_1z^{-1}} } \space { 1-b_2z^{-1} \over{1-a_2z^{-1} } } \space
{1 - b_3z^{-1} \over{ 1 - a_3z^{-1} } }
$$

where

$$
b_1 = 0.98223157, \space b_2 = 0.83265661, \space b_3 = 0.10798089
$$
$$
a_1 = 0.99516897, \space a_2 = 0.94384177, \space a_3 = 0.55594526
$$

Some judicious choice of sequencing for those sections, combined with some choice of gain factors for each section will be required to prevent overflow.  I haven't tried any of the other filters given in the link in Peter's answer, but similar considerations would probably apply.

White Noise

Obviously, the filtering approach requires a source of uniform random numbers in the first place.  If a library routine isn't available for a given platform, one of the simplest approaches is to use a linear congruential generator.  One example of an efficient fixed-point implementation is given by TI in Random Number Generation on a TMS320C5x (pdf).  A detailed theoretical discussion of various other methods can be found in Random Number Generation and Monte Carlo Methods by James Gentle.

Resources

Several sources based on following links in Peter's answer are worth highlighting.

The first filter-based chunk of code references Introduction to Signal Processing by Orfanidis.  The full text is available at that link, and [in Appendix B] it has coverage of both pink and white noise generation.  As the comment mentions, Orfanidis mostly covers the Voss algorithm.
The Spectrum Produced by the Voss-McCartney Pink Noise Generator.  Way down near the bottom of the page, after extensive discussion of variants of the Voss algorithm, this link is referenced in giant pink letters.  It's much easier reading than some of the preceding ASCII diagrams.
A Bibliograph on 1/f Noise by Wentian Li.  This is referenced both in Peter's source and by JOS.  It has a dizzying number of references on 1/f noise in general, dating all the way back to 1918.


A very fundamental question about average normalized power of noise If I have a average normalized power spectrum of noise is like this
$W(f)=N_0$ while frequency is $-B$ to $B$.

Why the average normalized power is $1/4 N_0 (2B)$ insteat of $N_0 (2B)$?

If a filter with transfer function $H(f)$ has a noise process as input, then
the output noise power $P_{\scriptstyle{\text{out}}}$ is given by
$$P_{\scriptstyle{\text{out}}} =
\int_{-\infty}^\infty |H(f)|^2 S_x(f)\,\mathrm df
$$
where $S_x(f)$ is the power spectral density of the input noise process.
If the input process is modeled as a white noise process, then $S_x(f)$
is a constant, call it $K$ for now, and so we have that
$$P_{\scriptstyle{\text{out}}} = K
\int_{-\infty}^\infty |H(f)|^2\,\mathrm df
$$
For an ideal low-pass filter with cutoff frequency $B$ Hz,
$$ H_{\scriptstyle{\text{lowpass}}}(f)
= \begin{cases} 1, &-B \leq f \leq B,\\
0, & |f| > B,\end{cases}$$
we have $P_{\scriptstyle{\text{out}}} = 2KB$. More generally,
for an arbitrary linear time-invariant system with maximum
passband gain equal to $1$, the noise-equivalent bandwidth of the system
is defined as
$$B_{\scriptstyle{\text{ne}}} = \frac{1}{2}\int_{-\infty}^\infty |H(f)|^2\,\mathrm df$$
and so we get $P_{\scriptstyle{\text{out}}} = 2KB_{\scriptstyle{\text{ne}}}$.
Put another way, the noise-equivalent bandwidth of an ideal LPF with
cutoff $B$ Hz is just $B$ Hz.

It is conventional in the engineering literature to write $N_0/2$ instead
of $K$ because then the output noise power expression becomes
$N_0B_{\scriptstyle{\text{ne}}}$. The constant of proportionality
$N_0$, measured in watts/Hz as Jason R pointed out, is convenient
for use in ordinary language when conversing with other engineers.
If we have a filter with noise
equivalent bandwidth $B_{\scriptstyle{\text{ne}}}$, and
the thermal noise power at the input of the filter  is $N_0$ watts/Hz,
then the output noise power is $N_0B_{\scriptstyle{\text{ne}}}$.
That's the way people talk during discussions. Nobody will say that an ideal
lowpass filter with cutoff $B$ Hz has "bandwidth" $2B$ even though
the transfer function has value $1$ for $f \in [-B,B]$ and thus
the support of $H_{\scriptstyle{\text{lowpass}}}(f)$ has length
$2B$.

Turning to your question, it is not a good idea to write $W(f) = N_0$
for $f \in [-B,B]$, but if the paper you are reading does insist on this notation,
then the output power is $N_0(2B)$. If the paper says that it is
$\frac{1}{4}N_0(2B)$, then it is mistaken. Another possibility is
that there are other conditions that are mentioned in the paper that
make $\frac{1}{4}N_0(2B)$ the correct answer, but you
have not included these conditions in your question. For example,
if the passband gain of the ideal LPF is $\frac{1}{2}$, then
the expression $\frac{1}{4}N_0(2B)$ for the output power
is perfectly correct.


Extracting the carrier signal from AM I am new to DSP and leaning about amplitude modulation.

I am trying to demodulate a 5 Hz signal that is modulated by a 50 Hz carrier. I am trying to extract the carrier, not the the 5 Hz signal.

A simple band-pass filter at the desired frequency frequency did not work and I am not sure why. I understand that phase locked loops can also be used.

Are phase locked loops a better approach to extracting the carrier signal? What may be some other useful approaches? I am looking for an approach that makes no assumptions on the carrier and signal frequency.

Thanks

The image below depict the AM signal along with the failed attempt to extract the carrier signal. Below that is the frequency components of the signal (which indicated that the carrier is clearly not suppressed).

To extract only the carrier with a band-pass filter, the filter pass-band has to be narrower than the modulating frequencies, or else some of the envelope modulation will "leak" through the filter.  Also note that a very narrow filter has a long settling time, so you will need a long enough test signal.


Oversampling and decimation: What filter to use? I have a sensor producing (more or less) bandlimited data with a cut-off of about 45Hz, with a roll-off and AWGN. I have an ADC that samples said signal at 800Hz, with a single-pole anti-aliasing filter at about 200Hz. The problem is, I only have enough communication bandwidth to send samples at 100Hz and therefore some decimation is necessary.

Currently, I simply have an 8-sample moving average filter and send every 8th sample. This feels dirty and suboptimal. Surely there must be a better way.

Is there an accepted "best" thing to do in this instance?  Should I, for example, do a low-pass FIR filter to squeeze out as close to 50Hz of signal bandwidth as possible?  Or, is there some sort of optimal estimation scheme will do better?

The aim is to implement several channels (9 channels) on a smallish microcontroller (ARM Cortex M4, for example), so the computationally cheaper the better!

Should I, for example, do a low-pass FIR filter to squeeze out as close to 50Hz of signal bandwidth as possible?

Yes, that is exactly what you should do.  That is an extremely low data rate, so even with a wimpy processor I would think that should be able to do a pretty good filter.  Especially since you only need to calculate $\frac{1}{8}$ of the filter outputs.


Estimate High Frequency I'm smoothing a signal with an little algorithme and I can control the amount of smoothing dynamically. Now it could be interesting to drive that amount with a measure depending of the frequency content. in other words if my signal is containing quick frequency it is interesting to have less of smoothing especially is the main signal is stationary.

What could I use for that? I tried mobile variance but I'm not sure this is the best method. I'm trying to avoid FFT methods or other that could dephase the signal.

second question, does kurtosis can be something relevant for that ? I think no but I would like to be sure.

thank in advance : )

PS: this isn't in real time, just sayin'

You could simply run your signal through some sort of high pass filter or shelf and compare the energy before and after the filter. If the ratio is high, you have a lot of high frequencies in your signal and should probably do less smoothing. If the ratio is low, there are not a lot of high frequencies and you can do more aggressive smoothing.

The frequency selectivity can be dialed in by the shape of the filter

As a simple example, the R code below runs two filters: a high pass one with unity gain at $f_s/2$ and a low pass one with unity gain at $0$.  Two test signals are used: one with more low pass content, one with more high pass content.

The output is:

[1] "X: 344.6108  LPF1: 32.1535  HPF1: 15.1717  LPF2: 17.4986  HPF2: 27.7821"

The low pass filtered test signal x_test_lpf retains more energy from the low pass filter than the high pass one (32.1535 vs 15.1717).  The high pass filtered test signal x_test_hpf retains more energy in the high pass output than the low pass one (17.4986 vs 27.7821).

R Code

#26668

x <- runif(1000,-1,1)

Nfilt <- 20

x_test_lpf <- c(x, filter(x,rep(1/5,5), circular = TRUE))
x_test_hpf <- c(x, filter(x,(-1)^(0:4)/5, circular = TRUE))

x_hpf_1 <- filter(x_test_lpf,(-1)^(0:(Nfilt-1))/Nfilt, circular = TRUE) # Unity gain at fs/2
x_lpf_1 <- filter(x_test_lpf,rep(1/Nfilt,Nfilt), circular = TRUE) # Unity gain at 0

x_hpf_2 <- filter(x_test_hpf,(-1)^(0:(Nfilt-1))/Nfilt, circular = TRUE) # Unity gain at fs/2
x_lpf_2 <- filter(x_test_hpf,rep(1/Nfilt,Nfilt), circular = TRUE) # Unity gain at 0

print(paste("X:", sum(x^2), " LPF1:", sum(x_lpf_1^2), " HPF1:", sum(x_hpf_1^2), " LPF2:", sum(x_lpf_2^2), " HPF2:", sum(x_hpf_2^2)))


Why do we need DFT when we already have DTFT? We already have DTFT so why do we need DFT?
Please don't write DTFT is used for discrete, aperiodic signal which transfers in the frequency domain to periodic and continuous, but in DFT is discrete periodic transfer in discrete periodic, these things I already know.

The answer is the same to the question: "Why do we need computers to process data when we have paper and pencil?"
DTFT as well as the continuous-time Fourier Transform is a theoretical tool for infinitely long hypothetical signals.
the DFT is to observe the spectrum of actual data that is finite in size.


What does a "vector" in a hidden Markov model mean? I know that a Hidden Markov Model (HMM) is used in speech recognition and understand it to some degree. However, what I don't know is how input (speech) is "transformed" to a vector which in later used in HMM.

How do you get a vector from a sound input? Is this vector readable by a human?

The way speech recognition is carried out with HTK (or any other tool) is sort of similar to the way speech recognition is carried out in the brain. When you hear a word, you instantly break it down into its constituent phones and then compare the phones with an internal mental "model" of the phones. These "models" are constructed over years of listening to speech and gives you the ability to distinguish between similar-sounding sentences like "How to wreck a nice beach" and "How to recognize speech". Speech recognition with HTK or any other model-based scheme works in a similar manner. Here, in a few steps, is how you do it:

You take the input speech signal and convert it into a feature vector representation.
Take a large number of sentences and perform step 1 on each of them.
Use the feature vectors in step 2 to build a statistical model for each of the phones/words in the sentences (there are a limited number of phones/words as against an infinite number of ways of saying them - so you reduce the unknowns by modelling).
When a new word comes in, break it into phones and compare with each of the known models. The sequence of phones with the highest probability wins!

All the above steps are critical to the successful completion of any speech recognition task. By decomposing a sound into its feature vector, you are taking it into a model-space, giving it a representation that makes it more suited to making-a-model-out-of than other representations (say the time-amplitude representation). Most such representations lie in the frequency, or the time-frequency domain. One of the most popular such representations is the MFCC (Mel Frequency Cepstral Coefficient). In a way, this technique mimics the human hearing response with a set of filters. An input signal is decomposed with this set of filters that have a logarithmic spacing of their center frequencies. The MFCC coefficients of any one sentence (say) are then used to model each of the phones that the sentence is made of. As an example, consider,

Sentence: HI.   Phonetic description: hh aa ey
When you feed the MFCC coefficients into HTK, it will associate the MFCC coefficients of a portion of the sentence with hh, another with aa and so on. When this is repeated many times over, the models for the phones begin to form.

HTK uses the tool HCopy to convert an input sentence into its feature-vector representation. There are many "flavours" to MFCC's as well (E_D_A or E_D_A_Z representations). It would be a good idea to read up on the documentation for HCopy within the htkbook.

The MFCC coefficients are written to a file with extension .mfc by HTK. It is not possible to read that file using any one of the text editors because (I think) the coefficients are written in binary. You can try to read the files with C though.

HTH.


How to isolate overlapping FM signals? Analogue radio stations often have a range of frequencies that they can use. Different transmission towers will transmit on different carrier frequencies in this range, because the signal takes time to propagate through air and if they used the same carrier then the received signal would be... Well, it wouldn't be well-formed FM, and it would almost certainly be echoey and distorted.

But what if they did all use the same carrier? You'd need multiple aerials placed some distance apart (would 3 be enough?) and would have to do something clever with them... But would it be possible to extract a signal from that?

If it's not too different a question, would it be possible to extract the individual signals if the transmission towers were transmitting something different to each other on the same carrier frequency?

What if they did all use the same carrier?

What you're describing is a single-frequency network (SFN). These are in common use for things that are not stupid FM broadcast.

The whole truth is that reception from different senders in an SFN just look like heavy multipath propagation, where the sent signal just takes multiple paths of different length to the receiver, and thus the receiver sees the signal as what would come out of a tapped delay-line.

If you consider it like that, the effect is clear: the different delays with which the signal can reach the receiver lead to something that for all practical purposes is an FIR filter. If you look at it from the viewpoint of digital communications, you get Inter-Symbol Interference (ISI) as delayed versions the last sent symbol still reach the receiver while the shorter paths already brought on the next symbol.

If you look at the same thing from the perspective of spectrum analysis, it is a frequency-selective channel.

Either way you look at it, you counter that with an equalizer that requires knowledge of how the channel looks (i.e. the different propagation delays, or the frequency-dependent behaviour of your channel, which are just two descriptions of the same thing, linked by the Fourier transform) and can then do its best to reverse that effect.

Now, sadly, FM radio is not digital communications. You don't get symbols, and by the very nature of FM  being a modulated carrier, you don't even get any kind of "known, not very repetitive" input with which you can estimate what your channel looks like.

This needn't be a problem – if the resulting delays don't make the channel too harsh, you don't lose much SNR at the receiver. We can measure that by comparing the delay spread (difference between shortest and longest path) to the inverse of the channel bandwidth. If the channel bandwidth is "sufficiently" smaller than the inverse delay spread, we're probably in luck and can still consider the channel "sufficiently" flat. "Sufficiency" here depending on your needs, so let's just do a order-of-magnitudes calculation.

An FM channel is roundabout 150 kHz wide. Invert that, and you get $6\frac23$ µs.

Electromagnetic waves travel with the speed of light in vacuum, and for all practical purposes here, in air, too. So those microseconds correspond to a distance difference to two transmitters of $c_0\cdot \frac{20}3\cdot10^{-6}\,\text{s}=3\cdot 10^{8}\frac{\text{m}}{\text{s}}\cdot \frac{2}3\cdot10^{-5}\,\text{s}=2\,\text{km}$.

Since we can be pretty sure that different FM radio transmitters aren't only 2 km apart, we can directly see that the resulting radio channel will be wildly non-flat and frequency selective. In fact, since the distance between these will be orders of magnitudes larger, we can assume that even given good channel information, we might have frequency-dependent deep fades that we simply can't properly "equalize away" (you get nothing through at those frequencies, so you can't even amplify anything to make it "flat" again).

In other words, FM and SFN don't play well together. There's been attempts to still implement SFNs for FM broadcast (I think Netherlands tried that), but as far as I know, they failed – they can only work by "aligning" the interference for limited locations, and if your network can't be too large and there's no strong variation of propagation environments for different transmitters (which would explain why small and flat Netherlands was an attempted adopter).

That's one of the most important reasons the world is leaving analog broadcasting behind. DAB and DAB+ are SFN-based, as their heavily multi-carrier oriented OFDM system with channel bandwidths of 1 kHz allows for much larger delay spread, and their frame structure and channel coding allow for heavy error correction, so that even deeply faded subcarriers don't matter that much.

So, yes, you've found a (if not the most important) spectrum management reason why FM should, if you ask me, die as soon as possible: it enforces the non-reuse of spectrum, and that's a waste of a precious resource.

Now, you mention multiple antennas:

Well, yes, if you have multiple antennas, you can use phasing to make them have some directivity, i.e. you can "simulate" having a single antenna that only receives from one direction. Hm. You buy multiple antennas, build a complex and large phased array and try to operate in a SFN, and the best-case scenario is you get the same performance as in a MFN with a single antenna. Doesn't sound too tempting, and, you'd need to have an idea in which direction you want to "look". It's notoriously hard to estimate that out of, again, an analogue signal with no known non-periodic structure. In a situation where you're at a fixed position (moving obviously breaks knowledge of directions) and adjust antenna phases just to get broadcast radio: Have you considered satellite radio?


Additive White Gaussian Noise (AWGN) and Undecimated DWT One of the benefits of DWT is that it is an orthonormal transform.

There are statements that the energy of noise component mainly concentrates on the high-frequency (detail) part and distributes homogeneously. The energy of noise component is included in more wavelet coefficients with smaller amplitudes, while the energy of the useful signal concentrates on fewer wavelet coefficients with bigger amplitudes. See Wavelet Denoising of Well Logs and its Geological Performance.

As far as I understand, that is due to the orthogonality property of DWT.

Does this behavior of white Gaussian noise also hold true:

for undecimated DWT (like SWT or MODWT).
for less redundant forms of DWT (comparing to undecimated version) - DTCWT, for example.

As far as I understand both of them aren't orthogonal transforms.

One property of Orthogonal Transformations is that White Noise stays White (Uncorrelated) under Orthogonal Transformations (One could say it's a property of White Noise).

Usually this property assist with deriving properties of operators (Things are simpler).

Yet this property doesn't not hold for Transformations in general which might "Color" the noise in the new base.


Inverse system of a box function I learned about the inverse system.

Suppose I have an impulse response $h(t)$ which is a box function.
($A = 1, T = 2$)

If I take the Fourier transform and referring the transform pairs,
$$H(j\omega) = 2\operatorname{sinc}(\omega),  \\ H(j\omega)H_1(j\omega) = 1, \\ H_1(j\omega) = \frac{1}{H(j\omega)} = \frac{1}{2\operatorname{sinc}(\omega)}$$
where $\operatorname{sinc}(u) \triangleq \frac{\sin(u)}{u}$ , right?

Then, what is the reverse of this Fourier Transform $\frac{1}{2\operatorname{sinc}(\omega)}$?

Or instead of using the Fourier Transform, should I use a different approach to get the inverse system of the box function?

As pointed out in the comments, the inverse system is not realizable. However, it is possible to write down an expression for the impulse response of the inverse system. We want a system that results in a Dirac impulse when convolved with the given rectangular pulse $h(t)$. This can be done as follows. We delay the rectangular impulse by $T/2$, so it starts at $t=0$. Then we add an infinite number of delayed versions of it such that the result becomes a step function $u(t)$. Each delay is an integer multiple of $T$:

$$u(t)=\sum_{n=0}^{\infty}h(t-T/2-nT)\tag{1}$$

The system that achieves this has an impulse response given by

$$f(t)=\sum_{n=0}^{\infty}\delta(t-T/2-nT)\tag{2}$$

Now we have a step function, from which we can generate a Dirac impulse by taking the derivative. So the inverse system is given by the concatenation of the system with impulse response $f(t)$ given by $(2)$ and an ideal differentiator. The total impulse response of the inverse system is then

$$g(t)=(f\star\delta')(t)=\sum_{n=0}^{\infty}\delta'(t-T/2-nT)\tag{3}$$

where $\delta'(t)$ is the distributional derivative of the Dirac impulse.

Clearly, the system with impulse response $g(t)$ given by $(3)$ is not realizable, but as a theoretical concept it is interesting, just like other ideal systems (e.g., ideal frequency selective filters, differentiators, Hilbert transformers, etc.).


mathematical transfer for different camera attributes Suppose I have two camera shoots(Image1,Image2) for the exactly same scene (each pixel is alignment correctly) but from different cameras and different light conditions. I want to remove this light difference. I tried two way till now:

S,V component adjustment

Gamma correction
I am getting result around 7/255 difference after adjustment.

I am looking for a mathematical transform "F", something like Adjusted_Image=F(Image1)  which leads to |Image1-Image2| ~ 0
The "F" function will contains must-tuned parameters (a,b,c,....).
Any ideas of the formula of "F"?
many thanks,

This Process is Known as Color Correction.It is Explained in more detail here.

To Give a short Explanation:

Lets suppose

You have a Input Color Matrix represented as

    | R_O1 G_O1 B_O1 |
A = | R_O2 G_O2 B_O2 |
    | R_On G_On B_On |

and a Reference Target Color Matrix represented as

    | R_T1 G_T1 B_T1 |
B = | R_T2 G_T2 B_T2 |
    | R_Tn G_Tn B_Tn |

Color Correction Matrix X is computed as follows:

            **A X = B**

    | R_O1 G_O1 B_O1 |     | A11 A12 A13 |   =  | R_T1 G_T1 B_T1 |
    | R_O2 G_O2 B_O2 |   * | A21 A22 A23 |      | R_T2 G_T2 B_T2 |
    | R_On G_On B_On |     | A31 A32 A33 |      | R_Tn G_Tn B_Tn |

X = (B.A^T)(A.A^T)^-1

Some one please help me editing the formulas!


Is there just a carrier in each frequency band in FDMA? The figure is the FDMA,and i know that FDMA is that dividing the bandwidth into many blocks(frequency band),i want to ask Is there just a carrier in each frequency band in  FDMA?i mean one carrier for User1 ,one carrier for User2,and one carrier for User3,and so on.Or for the user N,there are lots of carrier in it,if i want to send a signal,then i have to use a carrier?

That's undefined.

FDMA just says that different users use different parts of the spectrum.

What these users do within "their" frequency ranges is up to the system to define.


How to demodulate PM modulation I'm using a sin wave $s_{sc}(t) = d(t)\cos(\omega_{sc}t + \varphi_{sc})$ to modulate a carrier using phase modulation. I'm working on Matlab so I use the pmmod function to do so. It is not too complicated to understand what this function is doing, probably computing something like $\cos(\omega_c t + h s_{sc}(t) + \varphi_{ini})$ (with h the modulation index).

In order to demodulate it, I'm using the pmdemod function of Matlab. But for this one, I have absolutely no idea of what is happening in background. I tried to analyze the results obtained after a pmmod-pmdemod sequence, and this is what I observed on Matlab:

Regarding the convergence time, it seems that Matlab use a PLL in background. but what kind of PLL is it ? I'm very interested about the digital methods used to demodulate PM modulation and more particularly the maths behind it. But I cannot find any article/document/book about it (Or I just don't know what to search on google, because when I type "phase modulation" I just find article about CPM modulation or digital PSK/FSK modulation methods...).

Any help would be helpful :)

Phase modulation is demodulated using phase discriminators otherwise called phase detectors- anything that can translate a change in phase to a change in magnitude.

The most common for real functions is a simple multiplier followed by a low pass filter as given by the product of two sinusoids in quadrature:

From the relationship

$$\sin(\alpha)\cos(\beta) = \frac{1}{2}\sin(\alpha + \beta) + \frac{1}{2}\sin(\alpha-\beta)$$

Then by multiplying the phase modulated carrier with a local oscillator at the carrier frequency we get:

$$\sin(2\pi f_c t + \phi(t))\cos(2\pi f_c t)$$

Which after the low pass filter which eliminates the sum term we get:

$$= \sin(\phi(t))$$

When the phase deviation is small we get a linear result out of the discriminator since for small angles $\sin(\phi(t)) \approx \phi(t)$.

Note that such a product detector is also sensitive to amplitude (AM), so for purposes of phase demodulation only it is common to hard limit the modulated signal prior to the phase detector to remove all AM components.

The challenge in an actual receiver is to determine the carrier frequency accurately and track it through changing conditions such as Doppler and clock drift, thus a phase lock loop can utilize the phase detector and adjust the local oscillator (LO) such that the longer term average phase difference (error) between the LO and the received carrier of the modulated signal is zero, while the signal itself which is changing faster than the loop can track (outside the loop bandwidth) would pass through.

I have more details on other phase detectors including what is commonly used for complex signals at these posts:

Phase measurement

phase difference detection

Phase synchronization in BPSK

How to correct the phase offset for QPSK I-Q data


Difference between Fourier Transform and DFT? - Example I have read many excellent answers to similar questions, but never one this specific.  Here is another way to ask it.

Why is the modulation transfer function (MTF) of $\textrm{rect}(x/5) = \textrm{sinc}(5x)$ not equal to the DFT of a vector of zeros with 5 ones in the center?

Relative MATLAB picture and code:

x = -512:511; x = x / 1024;
sq = zeros(1,1024); sq(511:515) = 1;
sqMTF = abs(sinc(x*5));               %'Modeled'
sqDFT = 1/5* abs(fftshift(fft(sq)));  %'Measured'

I can make the two curves overlay each other by multiplying the DFT result by $\textrm{sinc}(x)$, which I intuit should be related to windowing of the original signal and related spectral leakage.

But shouldn't this be a convolution in the frequency domain, as opposed to a multiplication?

Jack.

The difference is between $\mbox{sinc}$ and the periodic version obtained using the DFT.

See this answer for a comparison.

  It strikes me that asinc = ratio of two sincs.

The $\mbox{asinc}$ function is a ratio of two $\sin$ functions:
$$
\frac{\sin(N\omega/2)}{\sin(\omega/2)} = \frac{\sin(N\omega/2)}{\omega/2} \frac{\omega/2}{\sin(\omega/2)}
$$
so, yes, you could see it as the ratio of two $\mbox{sinc}$ functions.

  Underlying question: can I model the sampling of a continuous signal (by, say, a digital camera) in the frequency domain as the DFT of the signal convolved with a Dirac delta train times the DFT of the pixel footprint (a rect)?  If so, why is the DFT of the pixel footprint an asinc instead of a sinc?

So you seem to be trying to model taking a picture with a digital camera?   Camera modeling is a little tricky. Even Gonzales and Woods seem to gloss over it.

The fundamental issue is the DFT of a $\mbox{rect}$ ($\Pi$) is a $\mbox{asinc}$.

If you're doing a discrete-time Fourier transform (DTFT), then it's not, but usually when dealing with computed FTs, you want the DFT.


What's the best technique to determine if the contents of two images are the same? I need to run comparisons on the output of hundreds of frames of simple 3D-rendered characters. Perceptually, many of the frames are exact duplicates. Based on this, I assumed something like a perceptual hash and comparison using pHash would work great. And in some cases it does...

However, I'm noticing that it doesn't tolerate minor (read: practically impossible for the human eye to see) color or contrast variations. Because of the nature of the rendered frames I received, there are some frames that are apparently different enough for pHash to conclusively decide that they are not the exact same.

Here are two example images along with a difference image:

pHash gives me a 0.473464% deviation between these images.

All those variations are enough to throw off pHash. Actually, pHash is able to tell that they are very similar, but I need a consistent low threshold to use for programmatic comparison. Other frames which actually deviate a tiny bit are sometimes considered more similar than the example above, even though they have perceptual differences. Which means I can't rely on pHash alone.

Here are two frames that are actually different (notice that his hands have moved a tiny bit):

Yet pHash gives me a difference percentage of 0.0188496%, which is less than 0.473464% from above, meaning that the frames that are perceptually different are actually considered more similar according to pHash.

I tried scaling and blurring the images a good bit to reduce any color variations, but I ended up having to do a ton of blurring to get closer. I don't want to do this too much, because I do want significant color changes to be considered as a unique frame.

For example, if there is a flash of light that illuminates the character, even if his character model is rendered in the exact same position, those frames should be considered unique:

Next, I tried Canny edge detection. This worked okay. Two problematic frames that are, for all intents and purposes, precise duplicates exhibited less differences:

The only issue I have with this method is that it only considers edges, which normally would be sufficient. If the character model changes at all, different edges should be detected. But this pretty much ignores significant color changes as mentioned above.

So my question is what is the best way for me to compare images like this and effectively counter insignificant color/contrast/noise variations? Seemingly pHash is ideal because of its perceptual nature, but I'm uncertain what I need to do as far as preprocessing to guarantee consistency from pHash (or some other algorithm, if it's better).

Any suggestions?

A simpler approach (at least to start with) would be to evaluate the absolute difference between the two images and specifically, its mean value or maybe even just the sum of it (across all pixels), and then take its histogram across your sample of pairs of images.

The histogram will follow a bimodal (or multimodal) distribution with two main modes: The lower mode corresponding to the pairs of images with slight differences and the higher mode (higher values of difference, therefore more error) to pairs of images with drastic differences. You can then use this information obtained from the histogram to work out a threshold. For more information please see: http://en.wikipedia.org/wiki/Multimodal_distribution

A few details:

You don't have to work on the colour images for this, you can simply
convert them to greyscale. The way you do this will depend on the
tools you are using.
The difference image that you depict in your post seems to be a
bitmap, the absolute difference between two grayscale images would
also "look like" a grayscale image, please see:
http://homepages.inf.ed.ac.uk/rbf/HIPR2/pixsub.htm

Hope this helps, happy to amend the response if its not exactly you are looking for, I was a bit surprised to see you going for edge detection straight from the beginning so maybe I have not understood exactly the definition of "frame difference" you might be thinking about.


Time series forecasting on dsp embedded systems I have a time series analysis algorithm that uses arima modeling and fits a model to the time series and gets the residual values. It is more of a statistical signal processing algorithm than a digital one. I have a matlab code for it which uses functions from econometric toolbox.

I need to program it on an embedded system. My first instinct was to do this on dsp embedded system, but I couldn't find anything about parametric methods in the dsp embedded books. Most of the books only cover the FIR, IIR and adaptive filtering. Is it possible to do something like this on embedded systems?

and if yes, assuming that I use one of ARM cortex M series microControllers, is there any specific header file that I can use?

I hope my question is clear enough but I can explain it more.

Thanks,

Giselle

clarification:

I was trying to say that what I have is more into parametric estimation part of digital signal processing but I guess I said it the wrong way.

I have a machine learning algorithm (more specifically an outlier detection algorithm). my input is a time series data, lets say from an ADC and it is continuous but I am doing segmentation on the data. I think the segmentation will also help with the computational complexity and the limited RAM resources (If that was what you meant by estimating required RAM resources).

The limited resources of microcontrollers is my exact problem. I do have libraries in both C++ and python that will provide me with functions that I need for my code, but adding them all in the code will make it so huge. and that is why I thought maybe dsp systems are my best choice.

I'm not sure why you oppose statistical signal processing to digital signal processing, they are not mutually exclusive. If you perform statistical signal processing with non-continuous quantized samples, you perform digital signal processing.

Second of all, if you are using a microcontroller, you are likely to have limited RAM resources. Can you clarify what kind of algorithms you are using? You need to estimate the RAM resources required for your algorithms.


Sampling below Nyqust Frequency My understanding of Nyquist's theorem is that you need to sample at a rate that's twice the bandwidth of a signal to fully recover it.

What exactly does the bandwidth mean in this context and what happens if you sample below this frequency?

I am trying to measure a physical signal that is filtered at 100 kHz, but which only has content I care about below 1 Hz. Since I only care about the low frequency components, can my bandwidth be considered 1 Hz for this purpose? Is there any reason to sample the signal at 200 kHz?

One does need to (not a necessary condition)! But if you do, you can recover (at least theoretically) the continuous signal from the discrete samples (sufficient condition).

The bandwidth most usually  denotes the span between the zeroth and the maximum frequency (in the analog signal). Sometimes, for a band-pass signal,  it denotes the span between the minimum (positive) and the maximum frequency, and there are similar theorems: twice the bandwidth (sometimes a little more depending on the location of the minimum or the maximum) is enough.

If you sample below, aliasing might occur, but not every-time, this depends on the structure of the spectrum.

You care only about 1 Hz, it is dangerous to simply resample: it is advised to filter a little above 1 Hz (say 1.5) and then resample consistently with 1.5 and the precision of your filter. Could be 4 Hz (above $2*1.5$ Hz). Those figures should not be taken for granted, just an idea of the whole method.

Why what it set to 100 kHz initially? Several possibilities, for instance:

people did not know initially in which part of the spectrum the information was located.
the acquisition scheme  was designed in a "the more data you have/the more information" mood


Exam Question - end-to-end delay from video capture to video playout Hi I'm preparing for my exam tomorrow and I've come across a 17 Mark question which seems trivially easy, so I think I'm missing a key point giving the rest of this exam is fairly difficult.  The question is below and my attempt is below that.
Question

The frame sizes are below:

My attempt:
I looked up video capture on wikipedia and it defined it as the process of converting an analog signal to digital - so I took the point at which frames were captured as being after they were digitised.  We are told to ignore the time taken for compression and decompression, so the next part is just to send the data.  The total amount of data in the frame sequence is 54KB (I'm not sure why the frame sequence starts and ends with an I-frame to be honest - doesn't seem to make sense).  So to transfer 1s of video will take 0.27s.
To be play a second of frames smoothly the the playout buffer must contain 25 frames minimum.  The sequence of frames is IBBPBBI which has a total size of 54 KB.  In one second the playout is 54 x 25 = 0.432 MB/s.
I'm pretty sure that we would need to make some estimations to solve the rest of this.
As always I'm open to any suggestions!  Thanks for reading :D

You can't start playing any B frame until the next following I or P frame has been scanned, compressed, transmitted, received, and decompressed.


Python toolboxes for state-space estimation via subspace estimation Is there (open-source) toolboxes for state-space estimation via subspace estimation in Python?
I know this is used in Matlab's n4sid function, but I didn't found any Python's implementation (even in Scipy's signal module or in Python Control Systems Library).
It would really helps not to have to code it :)

You can use pyvib to do frequency based subspace identification. Beware that there is no estimation of the initial state. It is possible to do optimization of the identified model, if the data is not perfectly linear. See the implemenentation, maybe you can use it, in case you want to do your own implementation.

Somewhat incomplete example. Take a look at the example to get a working code:

git clone --depth 1 https://github.com/pawsen/pyvib.git
export PYTHONPATH=pyvib

from pyvib.subspace import Subspace
from pyvib.signal import Signal
# partion the data so the format is (npp,inputs/outputs,R,P)
sig = Signal(uest,yest,fs=fs)
sig.lines = lines  # which freq lines should be used
sig.bla()  # computer best linear approx.
nvec = [2,3]  # model size to scan over
maxr = 5  # max number of rows
model = Subspace(sig)
models, infodict = model.scan(nvec, maxr)
errvec = model.extract_model(yval, uval)  # extract best model on fresh data
model.estimate(n, r)  # or do direct estimation if you know sys size

I have written the code. You are more than welcome to criticize and suggest improvements.


Transformation of coordinate set I have a problem in which i have a video that was taken from the side of the subject, something like this

 and i need to transform the coordinates of the subject to be as if the photo was taken from above like this

I have the real set of coordinates that i can relate to from picture 2 to picture 1 , i only need to transform the coordinates of the moving subject, not the whole picture.

I tried using polynomial regression in order to approximate the right transformation, but it failed to work (i can see that it's not the true movement).

also i have a problem that the video is taken with a fish eye lens.

what is the right way to approach it? or if you can direct me to the appropriate papers? (i never done image processing before).

many thanks,

I can't help with the fisheye lens issue, but assuming that effect isn't too large, the mapping between the above images is a perspective transformation.

We can transform a point on the plane $z=1$ into 3D space using a $3\times 3$ matrix.

$$\left( \begin{array}{c} X\\ Y\\  Z\\ \end{array} \right) =
\left( \begin{array}{c} s_x\ \ h_x\ \ t_x \\
                         h_y\ \ s_y\ \ t_y \\
                         p_x\ \ p_y\ \   1 \\ \end{array} \right)
\left( \begin{array}{c} x\\ y\\  1\\ \end{array} \right)$$

$s$ is scaling
$h$ is shearing
$t$ is translation
$p$ is perspective tilting (positive number means positive side tilts away)

We can project this object in 3D space back onto the plane $z=1$ simply by dividing the coordinates by $Z$. This maps every point in space to the intersection of the plane $z=1$ and the line connecting the point and the origin. This is basically the process of taking a photograph.

$$X_{proj} = X/Z$$
$$Y_{proj} = Y/Z$$

What we want to do is find the reverse transformation that takes us from the projected coordinates back to the original coordinates where the points only existed in the plane $z=1$.

It turns out that we can solve for all of the transformation parameters using 4 points (the four corners of a rectangle, for example) and solving the following system of equations (see here for the derivation and some example MATLAB code).

Making certain assumptions like the object not changing size or shape ($s_x=s_y=1$, $h_x=h_y=0$) can simplify these equations.


Why Isn't the ML Estimator (MLE) in MIMO Spatial Multiplexing Obtained by the Least Squares Solution? In the simplest scenario of MIMO spatial multiplexing: $$\mathbf{y} = H\mathbf{s} + \mathbf{n}$$
where:     $\mathbf{s}=[s_0,s_1,...s_{M-1}] \\\mathbf{y}=[y_0,y_1,...y_{N-1}]$

$\mathbf{n}=[n_0,n_1,...n_{N-1}]$ is zero mean complex normal independent Random vector

$M$ transmit antennas, $N$ receive antennas where $N\geq M$

The maximum likelihood solution is:
$$\tilde{\mathbf{s}} = \arg\min_{\mathbf{s} \in QPSK^M} ||\mathbf{y}-H\mathbf{s}||^2$$

Now, my textbook says that due to the fact that in this scenario $H^{*}H$ is not a diagonal matrix, the least squares solution followed by proccesing is not optimal, but I can't see why this is true.
To my understanding, denoting the LS solution as  $\hat{s}=(H^{*}H)^{-1}H^{*}\mathbf{y}$
I can write:
$$  ||\mathbf{y}-H\mathbf{s}||^2 = ||\mathbf{y}-H\hat{s}+H\hat{s}-H\mathbf{s}||^2 = \\||y-H\hat{s}||^2 + ||H(\hat{s}-s)||^2 + 2Re[(y-H\hat{s})^{*}H(\hat{s}-s)]$$

Now, by substituing the LS solution the $Re$ part vanishes.
And the problem reduces to
$$\arg\min_{s\in QPSK^{M}}|\hat{s}-s|^2$$
So I don't get why the ML is not optimal is this case.
Can anyone please clarify?
Thanks

edit: screen shot from the book

You must attention to the written text.
It doesn't say the ML isn't optimal, what it says is that the problem isn't regular LS problem but Least squares problem with Constraints.

The constraints make analytic solution infeasible and hence in order to solve it one must go through any signal in the space of valid signals and mark the one which minimizes the LS problem.

What they mean that the regular LS solution result doesn't necessarily yield a solution in the space of valid space of signals and even taking the closes one to the solution isn't optimal.

So we have the problem the constrained problem:

$$\begin{align*}
\arg \min_{s} \quad & \frac{1}{2} {\left\| H s - y \right\|}_{2}^{2} \\
\text{subject to} \quad & s \in {QPSK}^{M}
\end{align*}$$

Which its solution is optimal.

What they say isn't optimal is solving:

$$\begin{align*}
\arg \min_{s} \quad & \frac{1}{2} {\left\| H s - y \right\|}_{2}^{2} \\
\end{align*}$$

And then find $ s \in {QPSK}^{M} $ as you'd do in MIMO.
Namely, the greedy method (Do LS while ignoring the constraint and then apply the constraint) isn't optimal.

Yet the problem is that there is no analytic or efficient way to solve the Constraint LS problem but doing brute force search which isn't feasible.

I agree the text isn't clear about making this observation.
Yet it well known in the optimization field that solving in this greedy 2 step method isn't guaranteed to be optimal.

Remark
I found a link to the PDF - Doron Ezri - MIMO OFDM Lecture Notes. What I talked about can be seen in part 5.4.


How to find out the noise level of a signal Is it possible to find the amount of noise or SNR of a signal from its time series output? For example; let x = randn(100,1)  and y = awgn(signal,-5,'measured'); be two signals. Here signal is an unknown source signal of length 10,000 data points.

In Matlab this is what I did: snr(x) = 0.5177 and  snr(y) = -34.3797 dB Shouldn't I get -5 dB as SNR?

I am unable to understand what these values convey as they are a ratio: How is 0.5177 coming? Since I am generating y with an SNR of -5 dB, so snr(y) should be -5 dB.  Please help in clarifying if these answers are right and how to make sense of it.

I just don't think you're thinking this through all the way, and using the toolbox in MATLAB before you understand the basics. SNR is signal-to-noise ratio and is defined as $\text{SNR}=\frac{\text{Signal Power}}{\text{Noise Power}}$. When you call snr(signal, noise) in MATLAB, all it is doing is calculating the signal power, mean(abs(signal).^2), and the noise power, mean(abs(noise).^2), taking the ratio of them and returning the result in decibels (dB).

In the comments you are asking about how to calculate the SNR of x (the noiseless signal)...well its dividing by zero because the noise power in a noiseless signal is zero. You'd be doing $\frac{\text{Signal Power}}{0}=\infty$, you can go ahead and try using the snr function by doing snr(signal, zeros(size(signal))) but I'm willing to bet it either returns Inf or gives an error.

Note: Read the MATLAB docs on snr(signal) and you'll realize it is probably not returning what you expect. Docs are here: https://uk.mathworks.com/help/signal/ref/snr.html


How to implement a digital oscillator? I have a floating-point digital signal processing system that operates at a fixed sample rate of $f_s = 32768$ samples per second implemented using an x86-64 processor.  Assuming that the DSP system is synchronously locked to whatever matters, what is the best way to implement a digital oscillator at some frequency $f$?

Specifically, I want to generate the signal:
$$y(t) = \sin(2\pi f t)$$
where  $t=n/f_s$ for sample number $n$.

One idea is to keep track of a vector $(x,y)$ which we rotate by an angle $\Delta\phi = 2\pi f/f_s$ on each clock cycle.

As a Matlab pseudocode implementation (the real implementation is in C):

%% Initialization code

f_s = 32768;             % sample rate [Hz]
f = 19.875;              % some constant frequency [Hz]

v = [1 0];               % initial condition
d_phi = 2*pi * f / f_s;  % change in angle per clock cycle

% initialize the rotation matrix (only once):
R = [cos(d_phi), -sin(d_phi) ; ...
     sin(d_phi),  cos(d_phi)]

Then, on each clock cycle, we rotate the vector around a little bit:

%% in-loop code

while (forever),
  v = R*v;        % rotate the vector by d_phi
  y = v(1);       % this is the sine wave we're generating
  output(y);
end

This allows the oscillator to be computed with only 4 multiplications per cycle.  However, I'd worry about phase error and amplitude stability.  (In simple tests I was surprised that the amplitude did not die or explode immediately--maybe the sincos instruction guarantees $\sin^2+\cos^2=1$?).

What is the right way to do this?

You're right that the strictly recursive approach is vulnerable to the accumulation of error as the number of iterations increases. One more robust way this is typically done is to use a numerically-controlled oscillator (NCO). Basically, you have an accumulator that keeps track of the instantaneous phase of the oscillator, updated as follows:

$$
\delta = \frac{2\pi f}{f_s}
$$

$$
\phi[n] = (\phi[n-1] + \delta) \mod 2\pi
$$

At each time instant, then, you're left with converting the accumulated phase in the NCO to the desired sinusoidal outputs. How you do this depends on your requirements for computational complexity, accuracy, etc. One obvious way is to just calculate the outputs as

$$
x_c[n] = \cos(\phi[n])
$$

$$
x_s[n] = \sin(\phi[n])
$$

using whatever implementation of sine/cosine you have available. In high-throughput and/or embedded systems, the mapping from phase to sine/cosine values is often done via a lookup table. The size of the lookup table (i.e. the amount of quantization you do on the phase argument to sine and cosine) can be used as a tradeoff between memory consumption and approximation error. The nice thing is that the amount of computations required is typically independent of the table size. In addition, you can limit your LUT size if needed by taking advantage of the symmetry inherent in the cosine and sine functions; you only really need to store one-fourth of a period of the sampled sinusoid.

If you need higher accuracy than a reasonably-sized LUT can give you, then you can always look at interpolation between table samples (linear or cubic interpolation, for example).

Another benefit of this approach is that it is trivial to incorporate frequency or phase modulation with this structure. The frequency of the output can be modulated by varying $\delta$ accordingly, and phase modulation can be implemented by simply adding to $\phi[n]$ directly.


Trying to identify the best method for IIR realizations I'm relatively new to implementing IIR filters in raw programming formats (i.e. C or C++) and I've run into the concept of IIR filter realizations.

So my question: is there a 'best format (works 99% of the time)' for implementing these filters? i.e. an implementation that will typically be 'best' regardless of the application? Or is the actual implementation highly dependent on the application?

Here's an array of smaller yet related questions:

Is there ever an application in IIR filters where implementing the direct form of the filter is beneficial? Or should I always attempt to realize the filter as a cascade of subsystems?
Regarding cascade and parallel realizations, is cascaded always better?
Assuming that the exact filter implementation is dependent on the application, is it worth the effort to spend time and analyze the filter implementation from the beginning, or should I always aim for implementing as a cascaded subsystem and only deviate from that if the results are poor?
What is the purpose of transposed systems? I understand mathematically the differences, however, are there actual physical benefits of implementing a transposed vs nontransposed system?

I don't think "regardless of application" works here. Every application has different requirements and hence the optimum choice would be different.

However in gross oversimplification

Direct Form I is the "safest choice". Direct Form II is slightly better since it has only half the amount of state variables. Stay away from Direct Form II and Transposed Form I since the transfer functions to the state variables are only determined by the poles and there can be massive gain scaling issues
Cascaded is the easiest choice. It's fairly simple to convert between poles and zeros and actual filter coefficients. Parallel requires calculating the roots of a large polynomial which is a numerically ill defined problem.
There are different ways to pair poles and zeros to create second order sections. Best one is this : http://scipy.github.io/devdocs/generated/scipy.signal.zpk2sos.html
Use floating point whenever possible. Fixed point IIR requires very careful analysis and application specific implementation.
Whether optimization is worth it really depends on the application. Requirements can be all over the place: precision, speed, memory, mips, ease of us, stability against outlier condition, etc.


why wavelet coefficients are zeros after decomposition I am learning wavelet theory for image processing. To understand the theory, I write one Matlab program to decompose one black-white image. The program is as follows

 Image = zeros(256, 256, 'uint8');
 Image(101:200, 101:200) = 255;
 figure; imshow(Image);
 [cA1, cH1, cV1, cD1] = dwt2(Image, 'db1');
 Image1 = [cA1, cH1; cV1, cD1];
 figure; imshow(Image1, []);
 [cA1, cH1, cV1, cD1] = dwt2(Image, 'db2');
 Image1 = [cA1, cH1; cV1, cD1];
 figure; imshow(Image1, []);

The first decomposition using the argument db1 produces zeros for all wavelet coefficients. The black-white image has the transition from 0 to 255 along horizontal and vertical directions and should have high-frequency component. Why are zero wavelet coefficients generated?
If I change the argument from db1 to db2, the result will show horizontal and vertical lines in the subbands.

Edit:
I decompose the image further and at level 3 edge-like structure can be seen on the subband. Wavelet is one kind of filter bank I read from one book. It means the frequency in the image can only be seen in the higher frequency band. From the below link, it told how to get frequency response of Daubechies filters.

http://www.mathworks.com/matlabcentral/answers/44879-daubechies-filters-frequency-response

I don't understand how is the frequency for images defined. From one book, it says that "Spatial frequencies are measures of the number of oscillations (cycles) from dark to light per unit of length of the image". How can I know the frequencies in the input image.

If you recall, db1 is the Haar Wavelet.  The Haar Wavelet takes either an average of pixels within local windows for the approximation coefficients (or the LL band), or a difference of pixels within local windows for the detail coefficients (or the LH, HL and HH bands).

Be advised that the input image that you specified only consists of two intensities: 0 and 255.  Also, you set a square grid within this image to be 255 and it is uniformly shaped.

For self-containment, this is what your test image looks like:

This uniformly shaped object within a square image is important as part of the reasoning why you are not getting any output for the detail images (HL, LH, and HH).

The best way to describe why you're seeing an output for db2 and not db1 can be shown visually.

This slide is from the University of Toronto's CS 320: Introduction to Visual Computing course, specifically, the Discrete Wavelet Transform lecture:

You are well aware that when you take the 2D DWT, you produce 4 sub-images that are half the resolution of the original image.  The first output of dwt2 are the approximation coefficients where each output pixel is an average of a 2 x 2 window.  The other outputs (second, third and fourth) are detail windows that take two pixels within the window and subtracts these with two other pixels in the window.

As such, the reason why you are not getting an output with db1 is because all of your calculations for the detail images will cancel out.  Specifically, you will get 2 x 2 windows of either completely 0, or completely 255, and when you calculate the detail images, you will get 0s for the output regardless.  You would take add two 0 values, or two 255 values and you would subtract these with two other 0 values, or two other 255 values respectively, thus causing the output to be 0 regardless.

The db2 wavelet is a more complicated transform which is a weighted sum of non-uniform coefficients and so you will certainly get outputs for the detail images rather than just a simple differencing of the 2 x 2 windows.

I would like to stress that if you have a more complicated shape that is non-uniform, db1 will certainly not give you a zero output.  Try it on any test image that comes with MATLAB, like with cameraman.tif.

Hope this helps!


Can the Hanning Window be represented in the time domain I know the Hanning window, $w_{Hann}(n) = 1 - \cos \left(2\pi \frac{n}{N} \right)$, is typically applied to discrete Fourier transforms. $N$ being the total number of points in a transient signal that is to be transformed with the DFT. However, I'm trying to have a look at how this window (which is what I use to measure) affects the spectral line shape of the resultant Fourier transform.

For now I simply say my transient signal is a damped oscillator of the form
$$f(t) = e^{-t/\tau} \cos(\omega_0 t).$$The Fourier transorm of which is well known to be a Lorentzian profile. However As I measure my data with a Hanning window I want to take the Fourier transform of
$$f(t) w_{Hann}(t)$$
This means I need the Hanning window in the time domain, which I assume to be simply
$$w_{Hann}(t) = 1 - \cos \left(2\pi \frac{t}{T} \right)$$ where $T$ is the averaging time for my Fourier transform.

Is this a good approach and a sounds assumption I can make with the Hanning window in a time domain rather than for a discrete index of points?

So what you mean is that you want the continuous-time Hann window instead of the discrete-time window?

$$w_{Hann}(t) = 1 - \cos \left(2\pi \frac{t}{T} \right)$$

is not correct, since it goes to 0 at 0.

Blackman & Tukey give the continuous time window as:

$$D_2(\tau)=\frac 1 2 \left (1+\cos \frac{\pi \tau}{T_m} \right )$$

for $|\tau| < T_m$ and 0 elsewhere.

Wolfram Alpha

https://en.wikipedia.org/wiki/Hann_function

https://archive.org/details/TheMeasurementOfPowerSpectra/page/n58


When should we use the full Gabor Function and when do we use only the real part? In this YouTube video, we see that the guy is only using the real part of the Gabor Filter to implement his project.

My questions are,

When should we use the complex Gabor Function?
When should we use the real part of the Gabor Function?
When should we use the imaginary part of the Gabor Function? Or, do we ever need it?

The real part of the complex Gabor function is symmetric. This means it will give a large response at the centre of symmetric features such as lines.

The imaginary part of the complex Gabor function is anti-symmetric. This means it will give a large response at the centre of anti-symmetric features such as edges.

The real and imaginary parts are sort of in quadrature. This means when you take the magnitude of the response to the complex Gabor filter it will be high at the location of both edge and line features, which can be useful.


Keyword spotting - constant phoneme length I wish to implement a keyword spotting algorithm (in speech), on the basis of what was published in this article (Apple's Machine Learning Journal).

The article describes a neural-network-based method for keyword spotting, in which a network is first trained to identify twenty phoneme-classes and is triggered if the the sequence of phonemes that were inferred by the network correspond to Apple's 'Hey-Siri'.

It states that:

  The microphone in an iPhone or Apple Watch turns your voice into a
  stream of instantaneous waveform samples, at a rate of 16000 per
  second. A spectrum analysis stage converts the waveform sample stream
  to a sequence of frames, each describing the sound spectrum of
  approximately 0.01 sec. About twenty of these frames at a time (0.2
  sec of audio) are fed to the acoustic model, a Deep Neural Network
  (DNN) which converts each of these acoustic patterns into a
  probability distribution over a set of speech sound classes: those
  used in the “Hey Siri” phrase, plus silence and other speech, for a
  total of about 20 sound classes

What is puzzling to me is the bolded line.
It seems that the network is fed a constant number of frames (MFCC features) without taking into account that phonemes are uttered differently (the length varies) even by the same speaker (let alone by different ones).

Do I understand it correctly?
What would be (or is), a good way of dealing with the variations in length of phonemes?

There is no relation between fixed frame length and varying phoneme length. Since frames are just 0.01 seconds and average phoneme length is about 0.1 seconds, there would be a frame located around a center of the phone. And in that point of time DNN detector that spans 20 frames around central frame will detect the phone of interest even if it does not have any time-varying architecture. The fully connected DNN they use could do that pretty reliably.

Of course they could use something more advance and adapted to varying length sequences, for example, multilayer bidirectional LSTM, that would result in better detection accuracy but at the same time it would be much more computation-consuming and thus battery-consuming, it would also affect reaction time, thats why they prefer simpler architecture.

Their approach would not work for some significant phoneme length variations. For example you can test that "heeeeeeey siiiiiiiiri" does not work that reliably. At the same time one could argue that "heeeeeeey siiiiiiiiri" is not the same as "hey siri" for a human either and should not activate Siri. For a balance between energy consumption and detection accuracy such simple solutions are a good choice.


Problem with Kalman filtering accelerometer data I'm having some trouble implementing a Kalman filter in MATLAB. I have an Android phone connected sending data from accelerometer for 10 seconds. After i have the data I take out the x-axis vector. What I want after that is to get better readings using the Kalman filter.

The output I am looking for is both the acceleration and the velocity of the phone. When I run the code I get almost no change in acceleration (is between -0.5 and 0.5 m/s^2) and way to much change in velocity (increasing to very high numbers). The tests I have been running are when the phone is laying still on a table. So the correct values in acceleration and velocity should be 0.

The problem is that since i'm only measuring acceleration my state vector will be the same vector as the measurement vector with the only difference that i also have velocity in the state vector. Both the control input is the acceleration and the measurement is the acceleration.

So when I update my state vector the difference that is going to be multiplied with the Kalman gain is always zero.

So my question is: Should the measurement noise / the process noise be added to the measured acceleration "z" or should the noise be added to the state prediction vector?

If there is any offset in the acceleration measurement, then a Kalman filter won't help without the addition of some velocity and/or position reference measurements.  Double integration of random-walk noise in the acceleration measurement will eventually lead to a reported position of your device outside the solar system.


Continuous wavelet transform subsampling at lower frequencies? I am not completely confident in my understanding of wavelet theory, but since I am currently in the process of creating a CWT scalogram tool, I have to ask:

Is there a theorem allowing us to subsample at low frequencies and obtain correct scale/frequency results?

The core of continuous wavelet transform

$$ C_s(\tau,a) = \int s(t) \psi^*_{\tau,a}(t) dt$$

discretization is that (this is only one of the numerous possibilities) one can find constants $a_0$ and $b_0$, whose product   $a_0 b_0 <C_\psi$ ($C_\psi$ depending on the wavelet), such that the following sampling pattern

$$c_{j,k} = C_s(k b_0 a_0^j,a_0^j), \quad (j,k)\in \mathbb{Z}$$

keeps all the information contained in $C_s(\tau,a)$. The subsampling at low frequencies is bigger, i.e. one need less samples since large wavelets smooth the signal.

The varying sampling space from low (top) to high (bottom) is illustrated below:


Adaptive Digital Filter Block Diagram Question I'm currently attempting to study up on adaptive digital filters. My book presents the diagram I've included below and I'm having trouble understanding conceptually what it's indicating. The problem deals with noise cancelation. The idea is that someone is driving and makes a phone call. The x(k) is their voice input. There's a reference mic at v(k) which picks up road noise. I know that the ultimate goal is to filter road noise from our transmitted voice signal.

The desired output is obviously:

$$
d(k)=x(k)+v(k)
$$

The error in this case is:

$$
e(k)=x(k)+v(k)-y(k)
$$

Taking a quote from my book, it states

  If the speech x(k) and the additive road noise v(k) are uncorrelated with one another, then the
  minimum possible value for $e^2(k)$ occurs when y(k) = v(k), which corresponds to the road
  noise being removed completely from the transmitted speech signal e(k).

I don't understand how e(k) is our output of the system though. It seems to me that if we minimize our error, then it approaches zero. This means that $d(k)-y(k) = e(k)=0$ Consequently if our output is the error and we've minimized it, it seems like we're outputting 0 not a transmitted signal e(k) with the road noise removed!? I guess I'm asking why our desired output d(k) isn't our output....why is the error the output?

Can somebody help me understand this conceptually? Thank you for your help! Please let me know if I need to clarify anything.

Judging from the figure, the situation is slightly different from your explanation in the question. The noise $v(k)$ is the actual noise in the signal, not the noise picked up by the reference microphone. So the noisy signal is $d(k)=x(k)+v(k)$. If you knew $v(k)$ you could simply subtract it from $x(k)$ without the need to use an adaptive filter. What you have is another noise signal $r(k)$, which is a filtered version of the noise $v(k)$. This unknown filter is depicted by the "black box" in the figure. It is filtered because the transfer function from the noise source (road, tires, etc.) to the reference microphone is different from the transfer function to the microphone recording the speech. What the adaptive filter is trying to do is estimate the noise in the speech signal from the reference noise, i.e. it tries to invert the unknown filter in the black box. This can be achieved by minimizing the power of the error signal $e(k)$. The reason is that if you assume that speech and noise are uncorrelated, the output of the adaptive filter can only reduce the noise component in $d(k)$, not the speech component. So the power of $e(k)$ becomes a minimum if the output of the adaptive filter $y(k)$ equals $v(k)$. You don't need to worry that the speech signal is removed because $y(k)$ cannot model the speech signal at all, because noise and speech are assumed to be uncorrelated. So ideally the error signal $e(k)$ contains only clean speech. Note that due to noise and speech being uncorrelated, the power of $e(k)$ can never become zero. Neither can the error signal itself become zero (for all $k$).


Cancel Drift after numerical integration I am trying to integrate angular acceleration obtained from a set of accelerometers positioned specifically at opposite corners of a cube, based on the paper
EcoIMU: A Dual Triaxial-Accelerometer Inertial Measurement Unit for Wearable Applications

I am getting the angular acceleration on each axis.  This signal is quite noisy.

Then after integrating angular acceleration to angular velocity using the trapezoidal rule, I get signals which drifts heavily and randomly.

I understand that noise and also numerical integration are causing the effect. Other than low pass filtering the data, are there any other methods to reduce noise?

And the major factor for the drift is numerical integration, how can this be handled?

Numerical integration may reduce noise under some conditions: the noise is zero-average, ergodic, its properties do not vary too much over time. Since numerical integration seems linear, pre-filtering the acceleration with a linear filter, or post-filtering after integration are likely to yield very similar results (up to border effects), since linear operations commute.

The separation of a useful signal from a complicated drift is a complicated topic, without further information of the signal, the noise and the drift. These are a few venues you can use, alone or combined, to improve your signal:

filter the acceleration with more efficient (nonlinear) filters. For signals with time dynamics, Kalman filters may embed integration drift (Kalman filter: Modeling integration drift). In econometrics, people use a family of Holt-Winters filters to smooth or extrapolate time series with noise, drift or seasonality.
use more advanced numerical integration methods, like Simpson that can filter more noise or taking into account the dynamics, see for instance in Direct numerical integration methods, the Wilson Theta method:

  the Wilson Theta method assumes that the acceleration of the system
  varies linearly between two instants of time

post-correct the drift (called trend, wander, baseline in other domains), with specific techniques based on sparsity of derivatives (BEADS: baseline placement) or using piecewise linear models like  LOWESS:
.


Moving average frequency response over an image I'm studying image denoising by linear filtering with cross-correlation, in particular with a moving-average kernel (K x K kernel of all equal elements which sum is 1).
For clarity, I'd like to refer to 1D moving average filters (and the 2D extension is just done by using 2D Fourier transforms):
I read the dsp SE answer at this link that finds the magnitude of the frequency response of a moving average filter, but that post uses a discrete time Fourier transform, whereas I've seen in a book the use of the discrete Fourier transform: the only difference is that the discrete Fourier transform version is obtained by sampling the discrete-time Fourier transform $H(\omega)$ in the points $\omega_k = \frac{2\pi \cdot k}{N}$.
The book I'm reading also states that the discrete Fourier transform has "meaning" only when $k$ is such that $k \in [-\frac{N}{2}, \frac{N}{2}]$ (and of course $k$ is an integer).
I don't understand how to interpret the sequence given out by the DFT with respect to the frequency response of that filter and why the DFT has "meaning" only if $k \in [-\frac{N}{2}, \frac{N}{2}]$ and not $k \in [0, N-1]$. I also suppose that N is the size of the 1D moving average filter, but I'm not sure.
EDIT: moreover, how does the size of the moving average filter relate to the frequency response?
Thanks in advance.

OK.  There's a lot of things going on in your question, and they aren't really pertinent to the title.  Here goes:

... but that post uses a discrete time Fourier transform, whereas I've seen in a book the use of the discrete Fourier transform ...

OK.  The book was taking a short cut.  Hopefully it explained what it was doing and why.  For an infinitely-large image, the discrete-time Fourier transform (DTFT) is exactly correct; for any image the discrete Fourier transform (DFT) will have undesired edge effects.

the only difference is that the discrete Fourier transform version is obtained by sampling the discrete-time Fourier transform $H(\omega)$ in the points $\omega_k = \frac{2 \pi k}{N}$.

This is incorrect.  It is a not-bad approximation if some conditions are met:  For a time-domain signal $x(n)$ that is only nonzero for some interval $n \in [n_{min}, n_{max}]$ and $n < N$, and which has a DTFT whose spectrum is very small for $\omega \notin [-\pi, \pi]$, then the approximation is good to the extent that the spectrum outside of $[-\pi, \pi]$ is small.

I don't understand ... why the DFT has "meaning" only if $k \in [-\frac{N}{2}, \frac{N}{2}]$ and not $k \in [0, N-1)$.

The DFT has meaning for any contiguous interval of $N$ points (and even sets of $N$ points that aren't contiguous, but meet certain conditions -- but that's getting really weird).  But that meaning gets further and further away from easy physical truths the further you get from $k \in [-\frac{N}{2}, \frac{N}{2}]$
For any DFT, $H(k) = H(k \mod N)$.  This drops out of the defining math.  It is why the DFT only has "meaning" for some set of $N$ points -- that's just all there is.  This is a reflection of the fact that the DFT is a lossless transform, and you start with $N$ points.  Once you look at $N$ output points, there cannot be more information in additional points.
For any DFT of a real-valued (not complex) signal, $H(k) = H(-k)^*$, where $X^*$ is the complex conjugate of $X$.  So for that case, while the DFT output has "meaning" outside the interval $k \in [-\frac{N}{2}, \frac{N}{2}]$, it's not a terribly useful meaning.  In fact, because $H$ is complex, and it contains $N$ non-zero real or imaginary components in the interval $k \in [0, \frac{N}{2}]$, that's often the only part you look at if you're starting with real-valued data (and there's FFT routines that are designed to take real-valued data and be just a bit faster than their fully-complex counterparts because of the inherent symmetry).

moreover, how does the size of the moving average filter relate to the frequency response?

Generally speaking, if you have two filters of the same overall shape, the broader one has a narrower frequency response, and visa-versa.  There's an exact relationship when you're talking about a continuous-time signal and taking it's Fourier transforms.  Everything becomes hand-wavy when you get into discrete time, but the principle more or less applies.
So -- a wider moving average filter has a lower cutoff frequency, pretty much proportional to the number of points.


Resize High Resolution Camera Capture to Low Resolution Camera Capture? Say, you captured a high-res image of a scene and a low-res image of the same scene, using a lower-res camera, ensuring that the lower-res camera is using the same technologies, except for the image sensor's resolution.

If I wanted to drop my higher-res image to the same lower-res image via digitally processing the image, what would be the right approach?

I reckon interpolation or naive downsampling are not the right way-to-go. They seem kinda illogical. I was thinking of using the RAW images and, perhaps, deploying these downsampling therein, and, then, using the pre-processing algorithms for converting the image into the final output. Is there some standard approach available, though?

If you are using a “single-shot” sensor you could simply just not debayer the raw channels.  That would cut your sampling without adding the extra overhead that is already baked into the image acquisition.  A benefit of this approach is that you could probably boost your frame rate at the same time.


Is ideal band pass filter (brick wall filter) linear phase? I'm very new in digital signal processing.
I have multiple sensors and the way I filters the signals in post processing is:

take FFT of the signals.
put zero on out range of interesting frequency (like brick wall filter).
take IFFT of the spectrum.

I repeat this for each sensor, and compare the phase relation between them in time traces. Thus, the phase must not shift during the filtering process. I wonder the above method distorts the phase or not. Is the ideal band pass filter linear phase?

Assuming you are refering to LTI (linear time-invariant) systems to implement the filter.

The impulse response $h[n]$ of the ideal brickwall bandpass filter :
$$H(\omega) = \begin{cases} 1 ~~~,~~~ |\omega-\omega_c|<W \\ 0 ~~~,~~~\text{o.w.}\\ \end{cases}$$
is
$$
h[n] = 2 \cos(\omega_c n) \frac{ \sin(W n) }{ \pi n }
$$

which is real and even symetric about origin; i.e., $h[n] = h[-n]$,  hence it's a zero-phase (in the passband), but noncausal filter.

When this impulse response is truncated and shifted right enough to make it causal: $h_c[n] = h[n-d]$ then the resulting bandpass filter will be linear phase in the passband:
$$ H_c(\omega)= e^{-j\omega d} H(\omega)$$ with a linear phase term of $$\phi(\omega) = -\omega d$$

The phase is not defined for those frequencies where the frequency response is zero.

Note that the truncated filter is no more the ideal filter. Also note that you cannot shift the ideal filter to make it causal (requires infinite shift). And finally note that the zero-phase ideal filter is also a linear phase filter.


How to visualize and understand Gabor filter using MATLAB I am reading a paper, in which the author defines gabor kernels as :

$$
\psi_{\mu,\nu}(x,y)=\frac{\lVert \mathbf{k}_{\mu,\nu} \rVert}{\sigma^{2}}\exp \left(-\frac{\lVert \mathbf{k}_{\mu,\nu}\rVert^{2}\lVert \mathbf{z}\rVert^{2}}{2\sigma^{2}}\right)\times\left[\exp\left(i\mathbf{k}_{\mu,\nu}^{T}\mathbf{z}\right)-\exp\left(-\frac{\sigma^{2}}{2}\right)\right]
$$
where $\mu$ and $\nu$ define the orientation and scale of tghe Gabor kernels, respectively, $\mathbf{z}=(x,y)^{T}$, and the wave vector is defined as

$$
\mathbf{k}_{\mu,\nu} = \left(k_{\nu}\cos{\phi_{\mu}},k_{\nu}\sin{\phi_{\mu}}\right)^{T}
\quad\text{with}\quad k_{\nu} = \frac{k_{\rm max}}{f^{\nu}}, k_{\rm max}=\frac{\pi}2,f=\sqrt{2},\quad\text{and}\quad\phi_{\mu}=\frac{\pi\mu}{8}$$

Now, I want to plot the graph of this function using matlab to visually see and understand the Gabor filter. I have very basic knowledge in using matlab.

Specifically, I am interested in viewing these Gabor Kernels at single scale $(\nu=0)$ and four orientations $(\mu\in\{0,2,4,6\})$ with $\sigma=1$

Right now I am trying to understand Gabor filters. Any help will be greatly appreciated.

Thanks in advance.

The following is my implementation, of the above problem :

In my_gabor_filter.m, I have the following code :

function psi = my_gabor_filter(x,y,mu,nu,sigma)
    phi = pi*mu/8;
    f = sqrt(2);
    k_max = pi/2;
    k_nu = k_max/(f^nu);

    % This is the wave vector
    k_vec = [k_nu*cos(phi),k_nu*sin(phi)]';
    z_vec = [x,y]';

    k_vec_norm = norm(k_vec);
    z_vec_norm = norm(z_vec);

    exp1 = (k_vec_norm/(sigma^2));
    exp2 = (exp(-((k_vec_norm^2)*(z_vec_norm^2)/(2*sigma^2))));
    exp3 = (exp(complex(0,k_vec'*z_vec))-exp(-(sigma^2/2)));

    psi = exp1*exp2*exp3;
end

and, in plot_gabor_filter.m I did the following :

function plot_gabor_filter()
    mu = 0;
    nu = 0;
    sigma = 1;

    [X,Y] = meshgrid(-2:0.1:2,-2:0.1:2);

    Z = zeros(41,41);

    for i = 1:1:41
        for j = 1:1:41
            Z(i,j) = imag(my_gabor_filter(X(i,j),Y(i,j),mu,nu,sigma));
        end
    end

    surf(X,Y,Z,'EdgeColor','black')

    xlabel('x axis');
    ylabel('y axis');
    zlabel('g(x,y)');
end

and when I run plot_gabor_filter, I get the following output, which helped me in visualizing the gabor function :