IIR coefficients and difference equation implementation in C language I can't find anything on this topic, so either I'm in the wrong direction or else there is just nothing about it on the internet.

So let's say I have 3 $b_i$ coefficients ($b_0,b_1,$ and $b_2$) and 2 $a$ coefficients ($a_1$ and $a_2$) in a difference equation. How do I code the difference equation in C, so that I get the desired effect from the filter coefficients?

/**********************************************************************
*                                                                     *
*                                                                     *
*                                                                     *
*                        Copyright (c) 2005                           *
*                                                                     *
*                      Robert Bristow-Johnson                         *
*                                                                     *
*                     rbj@audioimagination.com                        *
*                                                                     *
*                                                                     *
*                                 EQ.c                                *
*                                                                     *
*                                                                     *
*   Utilities for computing EQ filter coefficients, plotting          *
*   frequency response, and EQ filtering sampled audio data.          *
*                                                                     *
*                                                                     *
**********************************************************************/

#include <stdio.h>
#include "__functions.h"

enum EQ_type
{
LPF,
HPF,
BPF,
BPF0,
notch,
APF,
peakingEQ,
lowShelf,
highShelf,
lowShelf1,
highShelf1
};

typedef struct
{
float Fs;
enum EQ_type type;
float f0;
float dBgain;
float Q;
float filterCoefficients[5];    /* b0, b1, b2, a1, a2   (a0 is normallized to 1) */
float graphCoefficients[6];     /* B2, B1, B0, A2, A1, A0 */
} sectionData;

int main(void);
int computeSectionCofficients(float Fs, enum EQ_type type, float f0, float dBgain, float Q, sectionData *data);
int clearArray(long N, float *array);
int plotSectionFrequencyResponse(int Npoints, float Fmin, float Fmax, float *frequencyValues, float *dBvalues, sectionData *data);
int filterBuffer(long Nsamples, int Nsections, float *input, float *output, float *filterStates, sectionData **sections);

#define BUFFER_SIZE 1024
#define GRAPH_SIZE 589
#define NUM_SECTIONS 4
#define FS 44100.0
#define FMIN 19.99955848414893
#define FMAX 19999.55848414893

#define F1 50.0
#define G1 0.0
#define Q1 1.0

#define F2 200.0
#define G2 0.0
#define Q2 2.0

#define F3 1000.0
#define G3 10.0
#define Q3 2.0

#define F4 5000.0
#define G4 0.0
#define Q4 1.0

sectionData section1, section2, section3, section4;         /* this has to be global */

void plotTheFrequencyResponse(void);
void runTheFilter(void);

int main(void)
{
computeSectionCofficients(FS, lowShelf1,  F1, G1, Q1, &section1);
computeSectionCofficients(FS, peakingEQ,  F2, G2, Q2, &section2);
computeSectionCofficients(FS, peakingEQ,  F3, G3, Q3, &section3);
computeSectionCofficients(FS, highShelf1, F4, G4, Q4, &section4);

plotTheFrequencyResponse();
runTheFilter();
}

void plotTheFrequencyResponse(void)
{
float freq[GRAPH_SIZE];
float graph[GRAPH_SIZE];

clearArray(GRAPH_SIZE, graph);
plotSectionFrequencyResponse(GRAPH_SIZE, FMIN, FMAX, freq, graph, &section1);
plotSectionFrequencyResponse(GRAPH_SIZE, FMIN, FMAX, freq, graph, &section2);
plotSectionFrequencyResponse(GRAPH_SIZE, FMIN, FMAX, freq, graph, &section3);
plotSectionFrequencyResponse(GRAPH_SIZE, FMIN, FMAX, freq, graph, &section4);
}

void runTheFilter(void)
{
float aFilterBuffer[BUFFER_SIZE];
float aFilterStates[(2*NUM_SECTIONS+2)];
sectionData *sectionPtrArray[NUM_SECTIONS];

sectionPtrArray[0] = &section1;
sectionPtrArray[1] = &section2;
sectionPtrArray[2] = &section3;
sectionPtrArray[3] = &section4;

clearArray(2*NUM_SECTIONS+2, aFilterStates);
clearArray(BUFFER_SIZE, aFilterBuffer);
aFilterBuffer[0] = 1.0;                         /* input an impulse to get impulse response */
filterBuffer(BUFFER_SIZE, NUM_SECTIONS, aFilterBuffer, aFilterBuffer, aFilterStates, sectionPtrArray);
}

int computeSectionCofficients(float Fs, enum EQ_type type, float f0, float dBgain, float Q, sectionData *data)
{
register float temp, *coef_ptr;
float A = __exp2(0.08304820237218*dBgain);
float A2 = A*A;
float Ap1, Am1;
float cos_w0 = __cos_pi(2.0*f0/Fs);
float sin_w0 = __sin_pi(2.0*f0/Fs);
float alpha = 0.5*sin_w0/Q;
float a0, a1, a2, b0, b1, b2;

switch (type)
    {
    case LPF:
        b0 =  0.5*(1.0 - cos_w0)*A2;
        b1 =  2.0*b0;
        b2 =      b0;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case HPF:
        b0 =  0.5*(1.0 + cos_w0)*A2;
        b1 = -2.0*b0;
        b2 =      b0;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case BPF:
        b0 =   0.5*sin_w0*A2;
        b1 =   0.0;
        b2 =  -b0;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case BPF0:
        b0 =   alpha*A2;
        b1 =   0.0;
        b2 =  -b0;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case notch:
        b0 =   A2;
        b1 =  -2.0*cos_w0*A2;
        b2 =   A2;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case APF:
        b0 =  (1.0 - alpha)*A2;
        b1 =  -2.0*cos_w0*A2;
        b2 =  (1.0 + alpha)*A2;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case peakingEQ:
        if (A > 1.0)
            {
            alpha *= A;
            }
         else
            {
            alpha /= A;
            }
        b0 =   1.0 + alpha*A;
        b1 =  -2.0*cos_w0;
        b2 =   1.0 - alpha*A;
        a0 =   1.0 + alpha/A;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha/A;
        break;

    case lowShelf:
        temp = 2.0*__sqrt(A)*alpha;
        Ap1  = A + 1.0;
        Am1  = A - 1.0;
        b0 =  A*(Ap1 - Am1*cos_w0 + temp);
        b1 =  2.0*A*(Am1 - Ap1*cos_w0);
        b2 =  A*(Ap1 - Am1*cos_w0 - temp);
        a0 =  Ap1 + Am1*cos_w0 + temp;
        a1 =  -2.0*( Am1 + Ap1*cos_w0);
        a2 =  Ap1 + Am1*cos_w0 - temp;
        break;

    case highShelf:
        temp = 2.0*__sqrt(A)*alpha;
        Ap1  = A + 1.0;
        Am1  = A - 1.0;
        b0 =  A*(Ap1 + Am1*cos_w0 + temp);
        b1 = -2.0*A*(Am1 + Ap1*cos_w0);
        b2 =  A*(Ap1 + Am1*cos_w0 - temp);
        a0 =  Ap1 - Am1*cos_w0 + temp;
        a1 =  2.0*( Am1 - Ap1*cos_w0);
        a2 =  Ap1 - Am1*cos_w0 - temp;
        break;

    case lowShelf1:
        temp = cos_w0 + 1.0;
        b0 =  A2*sin_w0 +    temp;
        b1 =  A2*sin_w0 -    temp;
        b2 =  0.0;
        a0 =     sin_w0 +    temp;
        a1 =     sin_w0 -    temp;
        a2 =  0.0;
        break;

    case highShelf1:
        temp = cos_w0 + 1.0;
        b0 =     sin_w0 + A2*temp;
        b1 =     sin_w0 - A2*temp;
        b2 =  0.0;
        a0 =     sin_w0 +    temp;
        a1 =     sin_w0 -    temp;
        a2 =  0.0;
        break;
    }

coef_ptr = &(data->filterCoefficients[0]);
temp = 1/a0;
*coef_ptr++ = temp*b0;                      /* b0/a0 */
*coef_ptr++ = temp*b1;                      /* b1/a0 */
*coef_ptr++ = temp*b2;                      /* b2/a0 */
*coef_ptr++ = temp*a1;                      /* a1/a0 */
*coef_ptr   = temp*a2;                      /* a2/a0 */

coef_ptr = &(data->graphCoefficients[0]);
temp = b0*b2;
*coef_ptr++ = temp;                         /* B2 */
*coef_ptr++ = -0.25*b1*(b0 + b2) - temp;    /* B1 */
temp = b0 + b1 + b2;
temp = 0.0625*temp*temp;
*coef_ptr++ = temp;                         /* B0 */
temp = a0*a2;
*coef_ptr++ = temp;                         /* A2 */
*coef_ptr++ = -0.25*a1*(a0 + a2) - temp;    /* A1 */
temp = a0 + a1 + a2;
temp = 0.0625*temp*temp;
*coef_ptr   = temp;                         /* A0 */

data->Fs = Fs;
data->type = type;
data->f0 = f0;
data->dBgain = dBgain;
data->Q = Q;

return 0;
}

int clearArray(long N, float *array)
{
register long n;

for (n=0; n<N; n++)
    {
    array[n] = 0.0;
    }

return 0;
}

int plotSectionFrequencyResponse(int Npoints, float Fmin, float Fmax, float *frequencyValues, float *dBvalues, sectionData *data)
{
register int n;
float logFreqStep = __log2(Fmax/Fmin)/(float)(Npoints-1);
register float phi, tau = 1.0/(data->Fs);
float *coef_ptr = &(data->graphCoefficients[0]);
register float B2, B1, B0, A2, A1, A0;

B2 = *coef_ptr++;
B1 = *coef_ptr++;
B0 = *coef_ptr++;
A2 = *coef_ptr++;
A1 = *coef_ptr++;
A0 = *coef_ptr;

for (n=0; n<Npoints; n++)
    {
    phi = Fmin*__exp2((float)n*logFreqStep);
    frequencyValues[n] = phi;
    phi = __sin_pi(tau*phi);
    phi = phi*phi;                                                                              /* sin(w/2)^2 */
    dBvalues[n] += DB_LOG2_ENERGY*(__log2((B2*phi + B1)*phi + B0) - __log2((A2*phi + A1)*phi + A0));
    }

return 0;
}

int filterBuffer(long Nsamples, int Nsections, float *input, float *output, float *filterStates, sectionData **sections)
{
register long n;
register int i;
register sectionData **section_ptr;
register float sampleValue, state1, state2, *state_ptr, *coef_ptr;

for (n=0; n<Nsamples; n++)
    {
    section_ptr = sections;                         /* reset section pointer */
    state_ptr = filterStates;                       /* reset state pointer */

    sampleValue = input[n];                         /* get input sample */

        state2 = *state_ptr++;                      /*   x[n-2]            */
        state1 = *state_ptr--;                      /*   x[n-1]            */
        *state_ptr++ = state1;                      /*   x[n-1] -> x[n-2]  */
        *state_ptr++ = sampleValue;                 /*   x[n]   -> x[n-1]  */

    i = Nsections;
    while (--i >= 0)
        {
        coef_ptr = &((*section_ptr++)->filterCoefficients[0]);  /* point to section filter coefficients */

        sampleValue =  *coef_ptr++ * sampleValue;   /*   b0*x[n]           */
        sampleValue += *coef_ptr++ * state1;        /*   b1*x[n-1]         */
        sampleValue += *coef_ptr++ * state2;        /*   b2*x[n-2]         */

        state2 = *state_ptr++;                      /*   y[n-2]            */
        state1 = *state_ptr--;                      /*   y[n-1]            */
        *state_ptr++ = state1;                      /*   y[n-1] -> y[n-2]  */

        sampleValue -= *coef_ptr++ * state1;        /*  -a1*y[n-1]         */
        sampleValue -= *coef_ptr++ * state2;        /*  -a2*y[n-2]         */

        *state_ptr++ = sampleValue;                 /*   y[n]   -> y[n-1]  */
/*
*      at this point sampleValue is the output of the section just completed
*     and will be the input to the section about to be computed unless this
*     is the last section, in which case sampleValue is the input.
*     state1 and state2 are the output states of the section just completed
*     and will be the input states of of the section about to be computed.
*/
        }

    output[n] = sampleValue;                        /* store output sample  */
    }

return 0;
}


How to measure a detector's performance under extremely complex real-world environment? I am trying to formalize my understanding of how best to measure the performance of a binary detector. (See previous question here). I had looked into ROC curves for doing this however based on previous research and feedback here I think ROC curves are too simplistic. Here is the problem I face:

Context:

I have a microphone that is setup in a park. The mic takes in inputs and an algorithm runs in the background that is supposed to detect the rustle of leaves. (Animal approaching, something in the bushes, etc). It does this by insert time-frequency-detection-scheme here, but that is not the issue. The end result of course is a 'Yes' or 'No', (So a binary detector), indicating whether or not leaves are indeed rustling.

Main Question(s):

How can/should I characterize this detector's performance, vis-a-vis true positive and false positive rates? (If that is the wrong way to go about it, then what is the right way?)
I understand that to generate a false-positive analysis, requires a false-positive model. In this case, the false-positive model is not something simple like white noise, but spans almost any signal possibility. Therefore, how can we model all those possibilities?

Thoughts:

My current solution/thought is to measure a true-positive rate and false-positive rate to characterize detector performance. (For given settings/thresholds).

Thus, for given and fixed detector settings, I test its true-positive rate by rustling leaves (within a specified range), and see if it detects the incidents. This will give me a true-positive rate. (A consistent 'Yes').

However I am stumped as to how to measure the false-positive rate.

I do not understand how to do this, because the mic is outside and subject to anything - cars, people, animals, thunder storms, just life. In other words, just about any combination or permutation of signals imaginable is possible. How then does one test its false positive rate?

on an ROC curve your plotting FP(x-axis) vs TP(y-axis) They are calculated by:

True-positives = # of correctly detected positives/# of actual positives

False-Positives = # of positives that are not true-positives/ number of positives

Now on your ROC curve, because you did not mention incuding this into your question above, you most likely did not include the "coin-flip" line. That is the line y = x (from x = 0-1)

This line represents how well your algorithm performs by just guessing randomly (Yes, No)

If your algorithm is a decent detector, the ROC curve will it will be above the coin-flip line, if it is below this line it is better at being an anti-detector, i.e., when there is no event and would be a poor detector for your application.

I've worked extensively on this for my master's thesis on onset detection. Below is an example of ROC curves plotted against the coin-flip line.


Applying discrete Wigner-Ville distribution to Lena picture? It seems that 2D Wigner-Ville has been used to test the compression of data.
The interference is however a problem.

Can someone send me a link to a code where this has been done?
I want to study the code and change it little study the compression and Wigner distribution.

Here, seems to be basic steps to make 512x512 matrix out of the picture.

So it is now only the application of the distribution to the image.
Then, JPEG2000 compression of the resulted image to test.
And make a plot PSNR(dB) as a function of Coding Rate (bits per pixel).

Where and how can you apply Wigner-Ville to the image?

I need to apply the function and Fourier transformation to every single pixel of the picture, here.

I have the Lena picture in 512x512 in Matlab, done by

lena512 = imread('lena512.bmp');
lena512 = double(lena512);

I am reading the package here and thinking where Wigner has been possibly been applied.
Not sure if used at all.
However, it speaks there something about autocorrelation so it can be possibly applied in code:

  The long-range correlations in the original series are eliminated by
  the wavelet transform, and that the wavelet coefficients form series
  with virtually no autocorrelation or cross-correlation.

Where Wigner Distribution has been applied to the Lena picture?

The developers of Lucid, which is based on time-frequency analysis, report interesting results. Technical details are provided in this PowerPoint file.


Modulation scheme for a noisy channel with a complicated transfer function I want to build an audio modem with a speaker in a room and a smartphone. The channel is quite bad:

Far field, so there is potentially significant reverb
Typical residential background noise
phone microphone quality, phone case, and room transfer function can be all over the place
microphone handling noise and hand partially blocking the mic

While the transfer function is unknown and can be quite variable (+/- 20 dB for narrow band features), it's safe to assume, that the channel is fairly constant during transmission.

On the positive side, it can be slow, I only need to get a handful of bits across. So I'm looking for something that's better than morse code, but doesn't require a locally flat transfer function.

Far field, so there is potentially significant reverb

That's basically a multipath channel with a large delay spread.

In communications engineering, the measures taken to revert that effect would be called equalization.

So, that's less of a modulation problem, but a channel estimation problem. Estimating channels with large delay spread / high frequency-selectivity is relatively hard, and so is equalizing them once you know the channel. You can deal with that by sending a known sequence ("preamble") against which the receiver correlates.

Since the result of that correlation is the autocorrelation function of the preamble, convolved with the channel impulse response, you can get pretty far with that. Especially considering you're free to pick a preamble that is spectrally white, i.e. has a sharp (dirac) impulse as ACF.

Based on that, you could build a zero-forcing equalizer, which is essentially simply taking that impulse response, inverting its Fourier transform, and transforming back to time domain.

Problem with zero-forcing equalizers is that you get strong noise amplification at the frequencies the channel is essentially "deaf". In audio engineering, that's what happens when the guy at the FOH doesn't pay attention, and suddenly you get a really unpleasant tone at some frequency when someone starts using their microphone.

Other options include MMSE equalizers, LMS adaptive filters, etc.; there's a large corpus of literature on both equalization for communications, and equalization for making audio work.

I'd actually, for simplicity reasons, recommend a different approach:

Instead of trying to equalize this very long, and hence very frequency-selective channel, you'd simply take it, and divide it into many narrower channels.

If you choose to pick these equidistantly and use a synchronous boxcar-shaped pulse shaping on each of these, you end up with OFDM (Orthogonal Frequency Division Multiplexing).

The idea is that instead of this complex equalizer that you need to use the full channel with a high bandwidth data signal, you get a lot of narrower channels that you divide your data to, and because they are so narrow, an equalizer is just a phase and amplitude correction – much, much easier!

And to boost the benefits, OFDM is mathematically just an IDFT - DFT pair, and can hence be implemented really efficiently using FFTs.

  Typical residential background noise

yeah, so noise. Some frequencies might be especially disturbed, but most of this noise will be relatively broadband.

So, channel coding with a relatively long code might be a good idea. Probably with an outer and an inner code, to convert burst noise to something benign and then throw a high-rate code on the rest.

  phone microphone quality, phone case, and room transfer function can be all over the place

Yeah, so that's covered by equalization.

  microphone handling noise and hand partially blocking the mic

burst Noise and fading.

  There may be mild time-variance, but probably not a lot

well, that's why you'll need to re-estimate your channel regularly.

If you want something to play with: GNU Radio is a software framework for live signal processing and can be directly attached to microphone sources or speaker sinks.

It comes with OFDM blocks, but there's also projects like gr-adapt that offer adaptive filters to compensate channels.

I'd start with the OFDM examples, and just replace the radio hardware block with a block that interpolates the complex baseband to match your audio sampling rate, and shift it up to a carrier frequency before discarding the imaginary part.

Do the reverse (add zero imaginary part, shift down to baseband, decimate to baseband bandwidth) at the receiver  – Tadah, OFDM audio testbed.

Pick an OFDM FFT length that is large enough to divide the channel into subchannels of  a couple dozen Hz bandwidth, at first (that system will then have very long symbols, and thus be sensitive to channel changes, but it should be robust against complex channel impulse responses).

Start with a really "easy" modulation – BPSK. If that works, and a constellation diagram says your SNR permits it, you can go for higher modulations like QPSK.

My guess is that staying within the confines of PSK is probably a good idea. OFDM isn't really easy on the analog chain – it requires a large dynamic range. In fact, I think I remember that if you want to avoid more than one in thousand signals clipping, you need a 10 dB headroom over average power.


How to prove: high-slope/more change is high freq part of signal Generally, everyone knows that "the part of signal with high-slope is the high frequency part."

What is the definition or theory behind it?
How to prove this well-known thing?
perhaps, I miss something apparent.

The Fourier transform of the derivative of a signal $x(t)$ is given by

$$\dot{x}(t)\Longleftrightarrow j\omega X(\omega)\tag{1}$$

where $\dot{x}(t)$ denotes the derivative of $x(t)$, and $X(\omega)$ is the Fourier transform of $x(t)$. From (1) we have

$$\dot{x}(t)=\frac{j}{2\pi}\int_{-\infty}^{\infty}\omega X(\omega)e^{j\omega t}\,d\omega\tag{2}$$

If we now assume that $x(t)$ is band-limited with cut-off frequency $\omega_c$ (in radians), then (2) becomes

$$\dot{x}(t)=\frac{j}{2\pi}\int_{-\omega_c}^{\omega_c}\omega X(\omega)e^{j\omega t}\,d\omega\tag{3}$$

From (3) we can determine an upper bound for $|\dot{x}(t)|$:

$$|\dot{x}(t)|\le \frac{1}{2\pi}\int_{-\omega_c}^{\omega_c}|\omega| |X(\omega)|\,d\omega
\le\frac{\omega_c}{2\pi}\int_{-\omega_c}^{\omega_c}|X(\omega)|\,d\omega\tag{4}$$

The upper bound (4) on the magnitude of the derivative of $x(t)$ increases with increasing $\omega_c$, so you can expect faster changes of $x(t)$ with increasing values of the cut-off frequency $\omega_c$.


How Does Gaussian Blur Affect Image Variance I am trying to find a way to analytically estimate how the noise (or variance in pixel values) will be affected for a zero-frequency image of random noise after applying a Gaussian blur.

In more detail, say I have an image of pixels with an average $ μ $ and random noise of variance $ σ_{0}^{2} $. The image has no features in it, but is simply that normal noise. I now apply a Gaussian blur (or Gaussian filter) with a kernel of standard deviation $ σ_{G} $. I would like a way to analytically estimate (and I know it will have to be an approximation) what the variance of the pixel values in the filtered image will be ($ σ_{f}^{2} $).

Wikipedia has an equation that clearly states what I need. However, they do not have any references and I will like to have a stronger reference for it. I have been unable to find this anywhere else. Any ideas where I'll be able to find this? Does anyone know if that is the correct equation? And maybe how to derive it?

Thank you in advance for your help! Every little contribution helps!

In Gaussian blur the value of each output pixel is calculated as a weighted sum of all input pixels:
$$\text{out}(x, y) = \sum^\infty_{j = -\infty} \sum^\infty_{i = -\infty} \frac{1}{{2\pi \sigma_G^2}} e^{-\frac{i^2 + j^2}{2 \sigma_G^2}}\text{in}(x+i,y+j).$$
We want to calculate the variance $\text{Var}[\text{out}(x, y)] = \sigma^2_f$ based on the variances $\text{Var}[\text{in}(x + i, y + i)] = \sigma^2_0$. We require that $\text{in}(x + i, y + i)$ are independent, which you should mention in your question. Based on the identity (variance of linear combination):
$$Var\left[\sum_i c_i X_i\right] = \sum_i c^2_i Var[X_i],$$
where $c_i$ are constants and $X_i$ are independent random variables,
$$\sigma^2_f = \sigma^2_0\sum^\infty_{j = -\infty} \sum^\infty_{i = -\infty} \left(\frac{1}{{2\pi \sigma_G^2}} e^{-\frac{i^2 + j^2}{2 \sigma_G^2}}\right)^2.$$
For a large $\sigma_G$, the squared Gaussian is rather smooth and the sum can be approximated by an integral as:
$$\sigma^2_f \approx \sigma^2_0\int^\infty_{-\infty} \int^\infty_{-\infty} \left(\frac{1}{{2\pi \sigma_G^2}} e^{-\frac{i^2 + j^2}{2 \sigma_G^2}}\right)^2\,\text{d}i\,\text{d}j
= \frac{\sigma^2_0}{4\pi \sigma_G^2}\\
\implies \sigma_f \approx \frac{\sigma_0}{2\sigma_G\sqrt{\pi}},$$
which is the same result as in Wikipedia. You can calculate the integral using Wolfram Alpha.

Note that the noise needs not be Gaussian for the above to apply.


Analyzing Spectrum of a signal? my_signal = 8192*cos(2*pi*6*(0:7)/8) + 1i* 8192*sin(2*pi*6*(0:7)/8);

my_signal_fft = fft(my_signal);

plot(abs(my_signal_fft));

This is matlab code. What do you notice about the spectrum? Is it symmetrical? Obviously not, because I think it has an imaginary part. Would you enlighten me for the analysis of this signal's spectrum and why having an imaginary part doesn't make it symmetrical? How is the other half gone?

Your signal has the form

$$x(n)=A[\cos(2\pi fn)+i\sin(2\pi fn)]=Ae^{i2\pi fn}$$

This is a complex exponential with normalized frequency $f$ (normalized by the sampling frequency). The bins of the signal's FFT of length $N$ correspond to normalized frequencies

$$f_i=\frac{i-1}{N},\quad i=1,2,\ldots, N$$

In your case $N=8$ and $f=6/8=f_7$ (using the Matlab indexation, which starts with index $1$). So the frequency of your signal $x(n)$ lies exactly on the FFT grid. This is why there is one peak at index $7$ and all other values are $0$ (apart from numerical inaccuracies).

If $x(n)$ were real-valued, its FFT would satisfy the following symmetry condition:

$$X_k=X^*_{N-k}$$

However, for complex-valued signals, such as yours, there is in general no symmetry.


Can I get some help plotting an IIR filter frequency response in Python? The problem I'm trying to solve is described here.

In my script, I used Euler's identity to calculate z^-n, such that z^-n = cos(omegan) - jsin(omegan). I find the magnitude of these values from a range of 0 to pi, and then plot them. However, the resulting graph aperiodically oscillates between 0 and ~ -8 dB, and I'm honestly not sure why.

I'd appreciate any help! The code is as follows:

import math
import matplotlib.pyplot as plt

# define constants
b0 = .125 a1 = .875 x = [] total = 0

while total <= math.pi: x.append(total) total += math.pi/256

return the magnitude of z-n

def find_z(n, omega):

# z = re^jw # e^jw = cos(w) + jsin(w) # z^-n = cos(wn) - jsin(wn)
ret = abs(math.cos(omega * n) - (1j * math.sin(omega * n))) return ret

h = []
for i in x: h.append(b0 / (1 - (a1 * find_z(-1, i))))

for i in range(len(h)): h[i] = 20 * math.log10(h[i])  plt.plot(x, h, color='blue', marker='o', linestyle='-',  linewidth=2, markersize=12)
plt.xscale("log") fig = plt.figure(figsize=(12, 6)) plt.show()

print(h)

We can't replace the learning experience of doing the actual homework, and coding questions should be posted on Stack Overflow, so below is not a complete solution but some hints in case this part wasn't clear:

The frequency response is the magnitude and phase of the complex result for $H(z)$ when you use $z = e^{j\Omega}$, where $\Omega$ represents normalized radian frequency (normalized as it is the radian frequency divided by the sampling rate). This is the unit circle on the z-plane.

That said, given the transfer function

$$H(z) = \frac{0.125}{1-0.875z^{-1}}$$

To get the frequency response simply $z$ with $e^{j\Omega}$ and solve for the magnitude and phase as you sweep $\Omega$ over the frequencies of interest (such as from $0$ to $\pi$ as given in the linked problem, which is the equivalent of sweeping from DC to half the sampling rate.

These links may help if the above is confusing:

$e^{j\omega}$ on unit circle

How/why are the $\mathcal Z$-transform and unit delays related?


Does Wigner-Ville distribution provide the highest resolution? For linear time-frequency representation, the resolution is limited by Gabor limit. However, if interference is allowed, the resolution may be higher. What is the possible resolution in such a case? In "Electroencephalography Feature Extraction Using High Time-Frequency Resolution Analysis ", they stated

  We  all  know  that  the  Wigner-Ville  distribution  achieves  to  the  lower  bound  of  the
  uncertainty  principle  and  have  the  highest  resolution  in  all  time-frequency  distributions

Why is this statement true?

The statement is true if you restrict it to all TF distributions in the Cohen Class. That is rather easy to see, because the Wigner distribution generates all other distribution in that class by convolution with a non-negative time frequency kernel and therefore attains the minimal uncertainty.

Beyond the Cohen class things get much more complicated. First of all, the very simple concept of uncertainty in the Cohen class is not necessarily useful for non-bilinear distributions. So it's not obvious how to even generalise the question. So statements that attempt to be more general than the Cohen class first need some clarification about mathematical terminology and possibly even concepts. Outside the Cohen class you will also find distributions without time-frequency translation invariance. If they're bilinear they can have a TF uncertainty that is just as good as the Wigner-Ville distribution. Dropping bilinearity too, you can even beat WV with reassignment distributions of all sorts.

So summing up, the statement is true with the (possibly implcit) restriction to the Cohen class, which is the most common theoretical framework for TF distributions. Without that restriction it is not accurate.


Upsampling and zero order hold Maybe to someone this question may seem quite easy but I have a problem grasping it. Namely, assume that we have two systems A and B with sampling rates $f_A$ and $f_B$ where $f_A>f_B$ and $R=\frac{f_A}{f_B}$. If we upsample system B by $R$ (i.e., $f_B=f_A$) using the zero order hold interpolation what is the relation between the time samples before ($n_{old}$) and after ($n_{new}$) upsampling? More precisely, we defined $\mathcal{N}_{old}=\{1,2,...,n_{old},..N_{old}\}$ and $\mathcal{N}_{new}=\{1,2,...,n_{new},..,R*N_{old}\}$ as sets of time samples before and after upsampling. Based on the definition of ZOH  following relation should hold:
$$
n_{old}=\frac{n_{new}}{R}
$$
i.e., every $R^{th}$ sample from $\mathcal{N}_{new}$ corresponds to the one time sample from $\mathcal{N}_{old}$
but if we for example take some $n_{new}$that is not multiple of $R$ we will get that $n_{old}$ is not an integer what doesn't make sense. What am I doing wrong?

Assuming a digital (discrete-time) sample rate conversion (upsampling) operation carried on old samples of a signal x[n] to upsample it by an integer factor R  by the utilization of a zero-order hold interpolation filter, you would essentially replicate (repeat) every new and empty (R-1) samples with the single old sample that they originate from. I guess the problem you have stems from the fact that you misinterpret the relation:
$$
n_{old}=\frac{n_{new}}{R}
$$

which is correct for continuous variables but not if you consider the fact that both $n_{old}$,$n_{new}$,and $R$ are integers. Therefore it cannot hold true for those samples $n_{new}$ which are not integer multiples of $R$

To remedy your confusion with the ZOH output index calculation, you shall better use the following interpretation as a relation between new and old sample indices:
$$ n_{old} = \lfloor \frac{n_{new}}{R} \rfloor$$

Where floor, with brackets $\lfloor$, $\rfloor$ is the mathematical operation that returns the largest integer that is less than or equal to its argument.

Or you could, equally, define a ruling to map old and new samples as:
index $n_{old}$ is equivalent to the integer $K$, for all those $n_{new}$ that satisfy $K R \leq n_{new} < (K+1)R$


improving pre-processing (noise removal) orange seven segment display for OCR reading About this project: Older images have the date they where taken placed on the image, in the bottom left hand corner. The goal is to extract this date and put it in the EXIF data of the scanned photo.

I am new to the challenge of working with images. I have a code playground and git repo I am working with and posting my results, and have made some progress. But I have hit a wall with dealing with the noise near the extracted text. I am reaching out for any assistance and direction to solve this noise problem.

What I have done thus far:
https://gitlab.com/51m0nt3r4aar/OCR_Img_Date_From_Scanned_Jpeg/blob/master/ImageManipPLayground.pdf

What I need help with: the Morph-Close did not remove the largest bits of noise, as seen in the last image, yet doing more morphs has not been successful. I am looking for any guidance to help solve this problem.

Thanks!

The shapes you are interested in are very regular while the noise you are trying to deal with will be less regular. Ideally, you will only have two shapes at the end of your processing: wide rectangles for the horizontal parts of the display and tall rectangles for the vertical parts. You can leverage this fact in how you shape your structuring element when you use opening to remove noise.

This code is for Matlab (or Octave, its free cousin), but should easy to implement in Python. Note that I screen grabbed the image from your GitLab file and resized it - so these dimensions likely won't match yours exactly. I also inverted the image so the foreground is white. Also, the images below were resized to 400% before posting.

pkg load image
img = imread("download.png");
img = img * 255;
imshow(img)

First, use a wide structuring element for opening so only wide rectangles are kept.

SE_h = strel('rectangle',[1,8]);
h_lines = imerode(img, SE_h);
h_lines_reconstructed = imdilate(h_lines, SE_h);
imshow(h_lines_reconstructed)

The repeat for tall rectangles:

SE_v = strel('rectangle',[8,1]);
v_lines = imerode(img, SE_v);
v_lines_reconstructed = imdilate(v_lines, SE_v);
imshow(v_lines_reconstructed)

Recombine and display:

combined = or(v_lines_reconstructed, h_lines_reconstructed);
imshow(combined);


Output of system given input and frequency response Given the signal: $$x(t) = −4\cos(2\pi 4t) + 3\cos \left(\pi t − \frac{\pi}{3}\right)$$

By sampling theorem $Fs = 2 Maxf = 8 Hz$

The question ask what is the output of system y[n] when the signal x[n] is the input and the frequency response is: $H(e^{jω}) = 2 + 2e^{j2ω}$

I try to do this, but i don't undestand how to separate magnitude and phase from this result:

Euler's formula is you friend :) $$e^{j\theta} = \cos(\theta) + j\sin(\theta)$$
After applying it, you will have $$H(e^{j\omega}) = \Re\{H(e^{j\omega})\} + j\Im\{H(e^{j\omega})\}$$ and then it's easy to find magnitude $$|H(e^{j\omega})|=\sqrt{\Re^2\{H(e^{j\omega})\} + \Im^2\{H(e^{j\omega})\}}$$ and phase $$\angle H(e^{j\omega}) = \tan^{-1}\frac{\Im\{H(e^{j\omega})\}}{\Re\{H(e^{j\omega})\}} $$


Object Recognition using bag of features using matlab I want to write a MATLAB program for simple object recognition using bag of features. In short, I want to first extract the features from an image, create a visual library using those features, then cluster up the features belonging to one part together, hence creating different parts. now use these parts for matching. i have the basic idea but i don't know much math behind this, also i don't know how i can go about implementing it. can anyone help me with the matlab code for this?

The, for now free available computer vision book, by Richard Szeliski contains an excellent introduction to this algorithm. Page 612.

EDIT:

I answered a similar question here on DSP with some example code.


Why do cutoff frequencies of a high-pass digital filter differ from the ones i used to transform it from a low-pass analog filter in MATLAB? I need to transform an elliptic analog low-pass filter to a digital high-pass filter but i'm having a problem because the resulting filter has lower cutoff frequencies than the ones i set. Here is the code sample from MATLAB:

Fs = 44100; %Sampling rate
Rp = 0.05;
Rs = 60;

Wp = 2*pi*10000; %Cutoff frequencies of my high-pass filter
Ws = 2*pi*9700;

WsE= Wp/Ws; %LP->HP transformation for the analog stopband frequency
WpE=1

k=sqrt(1-(WpE/WsE)^2); %Elliptic filter calculation of N
D=(10^(0.1*Rs)-1)/(10^(0.1*Rp)-1);
q0=(1/2)*((1-sqrt(k))/(1+sqrt(k)));
q=q0+2*q0^5+15*q0^9+150*q0^13;
N=ceil(log10(16*D)/(log10(1/q)));

[z,p,k]=ellipap(N,Rp,Rs); %Making the analog prototype
bE=k*poly(z);aE=poly(p);

[bE,aE]=lp2hp(bE,aE,Wp); %Transforming the low-pass into a high-pass
[bd,ad]=bilinear(bE,aE,Fs);

Ndigital = 20000;
[h,w]=freqz(bd,ad,Ndigital);
H=abs(h);

f = w/(2*pi)*Fs; % Ploting the filter
figure(1);
plot(f,20*log10(H), 'LineWidth', 2);
xlabel('Hz');
ylabel('|H(z)|');
pause

and i am expecting a filter with my desired frequencies but get this:

I can't pin-point the problem in the code. Why am i getting these results and how do i fix them?

The bilinear transform warps the frequency axis, that's why cut-off frequencies of the analog filter are not necessarily mapped to the correct frequencies in the digital domain. You can use the prewarping option by providing an optional parameter to bilinear.m. Check the corresponding mathworks documentation.


Finding local peaks in-between samples I have $n$ discrete samples of a seismic signal $y[n]$:

I want to find local maxima in the signal.

A naive test for if $y[n]$ is a maximum would be:
$$y[n]: maxima \textbf{ if } y[n] > y[n-1] \textbf{ and } y[n] > y[n+1]$$

However the maxima are probably located in between samples, e.g. there may be a maximum at $i=4.25$.

In order to find maxima in between samples I believe that I need to interpolate $y[n]$.

How do I find maxima using interpolation ?
What form of interpolation should I use ?

As you can see my signal is not very noisy, however it would be good if the method also did a bit of filtering so that the maxima exceed a treshold and have a certain width (no spikes).

My biggest problem however is just to find peaks in between samples.
Any suggestions for a good way to do this ?

Thanks in advance for any answers!

Getting a sub-sample resolution

A very cheap (in terms of code size) solution is just to upsample your signal. In matlab, this can be done with interp(y ,ratio). A slightly more complicated solution consists in naively detecting peaks ; and for each peak, fitting a parabola through y[peak - 1], y[peak], y[peak + 1] ; then using the point at which this parabola is maximal as the true peak position.

Regarding peak detection

A bunch of techniques which help:

As suggested by Hilmar, convolving the signal by a Gaussian or Hann window, the width of which is roughly equal to half the minimum interval you want to see between detected peaks. Since temporal accuracy seems essential to your application, make sure that you take into account the time delay introduced by the filtering, though!
Subtract to your signal a median filtered version of itself (with a fairly large observation window) ; and divide the result by a standard-deviation filtered version of itself. This gets rid of trends and allows the thresholds to be expressed in units of standard deviations.
For peak-picking, I formulate that using a "top-hat" filter. Define the top-hat filtered version of your signal as yt[n] = max(y[n - W], y[n - W + 1], ..., y[n + W - 1], y[n + W]) ; and use as peaks the points where y[n] == yt[n] and y[n] > threshold.

All this can be very efficiently implemented in Matlab with a few passes of nlfilter.


Options to enhance a poor image I have a C# desktop application. I am using EMGU as a wrapper for OpenCV.

I am receiving images in JPEG format that were originally taken by an analogue camera.

The quality of the image is poor (compared to digital cameras).

I want to see if I can use filters/methods derived from OpenCV/EMGU that could enhance the image.

What I mean by enhancement is to make the image a bit sharper.

I have played around with median, bilateral, contrast stretching to see if the image can be improved.

Nothing has worked.

Even if it can be confirmed to me that nothing I can do will enhance the image that will be an answer in itself. I can stop trying and persuade the client to use a digital camera instead...

This is the image:

This is the code that I apply to the original image to get the image you can see here:

        int factor = Shared.factor;
        int xVal = 288;
        int yVal = 360;
        using (Image<Bgr, Byte> color = new Image<Bgr, byte>(grid))
        {
            byte[] data = null;

              for (int x = 0; x < xVal; x++)
              {
                  for (int y = 0; y < yVal; y++)
                  {
                      grid[x, y, 0] = (byte)(grid[x, y, 0] / factor);
                      grid[x, y, 1] = (byte)(grid[x, y, 1] / factor);
                      grid[x, y, 2] = (byte)(grid[x, y, 2] / factor);
                  }
                }
                color.Data = grid;
                data = ImageToByteArray(color.ToBitmap());
        }

    private byte[] ImageToByteArray(Image imageIn)
    {
        byte[] data = null;
        try
        {
            using (MemoryStream ms = new MemoryStream())
            {
                imageIn.Save(ms, System.Drawing.Imaging.ImageFormat.Jpeg);
                imageIn.Dispose();
                data = ms.ToArray();
            }
        }
        catch (Exception ex)
        {
         //do something
        }
        finally
        {
            if (imageIn != null)
            {
                imageIn.Dispose();
            }
        }
        return data;
    }

As requested this is an original image but from a different view:

There is no way to regain any lost data at all. However, there are many phenomena that might distort the image in a recoverable way, not actually losing any data. You seem to try to undo a "point spread function", that is eg. some consistent blur. This would be not the only way an image can be distorted, however it is by far the most common.

In that case you would assume that the original image is convolved or filtered by an (initially unknown) point spread function. There may be the inverse operation to that. You could call this "equalization" of the frequency space of the image, then.

Not every inverse operation would be represented by an "edge sharpening" operation that is known a priori. You could maybe try to look at the image a virtually point-shaped light source gives, preferrably in front of a dark background. A laser pointer or a bright LED might serve as an example.


Position of centroids after convolution One of the properties of convolution is that the centroids of the functions being convolved should add. See equation (18) of Wolfram Alpha's section on convolution Convolution Properties.

For example, if we have a Gaussian centered at t1=300, and we convolute it another Gaussian centered at t2=200, the resulting Gaussian peak after convolution is centered at 500. MATLAB works fine for that.

I was trying to convolute Gaussian functions, one centered at 0 (orange), and other one at 500 (blue), but the resulting function is not centered at 500 in MATLAB but appears much later (yellow). How can one circumvent this problem in MATLAB? I want to keep one Gaussian at 0. Thanks.

'''W=[0:1:2500]'; % bin numbers

g1= normpdf(W, 500, 100); % First Gaussian with unit area (x, mu, sigma)

g2= normpdf(W, 0, 100); % Second Gaussian with unit area (x, mu, sigma)

Convg1g2=conv(g1,g2)

plot ( g1)
hold on
plot( g2)
hold on
plot(Convg1g2)'''

The reason is g1 is not a complete Gaussian pdf since you limited the range to not extend past 0.

The following will work for you:

W=[-2500:1:2500];
g1= normpdf(W, 500, 100);
g2= normpdf(W, 0, 100);
ConvG1G2 = conv(g1,g2, "same");
plot(g1)
hold on
plot(g2)
plot(ConvG1G2)

Results in:


Cepstrum: IFFT or FFT to enter the cepstral domain? I'm trying to implement source-filter separation of some pre-recorded speech in Max/MSP. I'm unsure of whether one uses an FFT to enter the cepstral domain, or an IFFT. Some of the textbooks and dissertations I've found that cover this topic say that I should perform an FFT to get into the cepstral domain and then an IFFT to exit back into the frequency domain; others say that I should perform an IFFT to enter the cepstral domain and then an FFT to get back into the frequency domain! The first way makes more sense to me, but something about my patch isn't quite working quite right, and I'm not sure what I'm doing wrong.

Here is how I understand the procedure (Please tell me if you notice anything incorrect):

1. Perform an FFT on the time-domain signal to enter the frequency domain.
2. Do Mel scaling (optional) and then logify.
3. Do an FFT (or an IFFT, depending on who you ask) to enter the cepstral domain.
4. Separate the low bins (the filter) from the high bins (the source / excitation signal).
5. Perform an IFFT (or an FFT, depending on who you ask) on both source and filter to exit back into the frequency domain.
6. Exponentiate both source and filter.
7. Modify amplitude of filter bins.
8. Multiply the modified filter with the source signal.
9. Un-Mel-scale (optional).
10.Perform an IFFT to exit back into the time domain.

In Step 3, should I use an FFT or an IFFT? (And then Step 5 will be the opposite of whatever Step 3 is.)

P.S.: I am not trying to modify phase information, only change some amplitude information of the filter portion of the signal (the lower FFT bins in the cepstral domain) before I add/multiply them back together.

For Cepstrum I always have used to this steps:

Apply hamming windows in the signal
Apply FFT
Get magnitude
Convert to log scale
Apply IFFT

The equation for cepstrum:

IFFT(log(abs(FFT(s))))

But you can use FFT or IFFT, take a look:

IFFT(log(abs(FFT(s)))) == real(FFT(log(abs(FFT(s)))))

The difference is the scale representation, if do you end using FFT you need extract just the real information, for both above equations you will get the same shape:

For IFFT(log(abs(FFT(s)))):

For real(FFT(log(abs(FFT(s))))):

This is a cepstrum shape example from 4096 points sine in 440hz sampled at 44100hz


question with chirped signal I have some difficulty understanding the following question.

I have written a code to plot continuous and discrete version of the chirped signal.
import numpy as np
from scipy.signal import chirp, spectrogram
import matplotlib.pyplot as plt
%matplotlib qt
N = 400
t = np.linspace(0, 2, N)
f0  = 0
f1 = 10
t1= 2
Fs = 8000
n = np.arange(N)
alpha = (f1 - f0)/t1
w = np.cos(2*np.pi*f0*t + np.pi*alpha*t**2)
plt.plot(t, w)
w_d = np.cos(2*np.pi*f0*n/Fs + np.pi*alpha*(n**2)/Fs)
plt.plot(n / (N/2),w_d, 'ro')

which gives,

Although N is completely taken randomly. It gave me a smooth curve for a continuous version. I don't know if my plot is correct or not. However, I am stuck now. what does the question mean? How to find the counts crossing the absicca? For zero-crossing, as far as I understood, the count would vary for N? I'd appreciate it if anyone could help me understand the question.

I can't interpret numpy code, but I don't think you should define an 'N' variable. For a 2-second signal you should have a total of 16000 samples. Here's how to do this in MATLAB (notice how I defined my "number of samples" variable ('Num_Samples')):
Fs = 8000;
Ts = 1/Fs; % Time between samples
T1 = 2; % Signal duration = two seconds
F0 = 0; % Start frequency
F1 = 10; % Stop frequency (at T1 seconds)
Num_Samples = 2/Ts; % Two seconds divided by 'Ts'
n = 0:Num_Samples-1; % Here's your time-domain index
Alpha = (F1-F0)/T1;
X = cos(2*pi*F0*n*Ts + Alpha*pi*(n*Ts).^2);

figure(1)
plot(n*Ts, X, '-bs'), grid on

Finding the number of zero-crossings in the 2 second duration of your chirp signal is nothing more than examining each generated 'x(nTs)' sample (starting with the 2nd sample) and seeing if its polarity is different from the previous 'x((n-1)Ts)' sample. Each time a polarity difference occurs you have a single zero-crossing.


Beat Detection with Arduino I'm trying to make a robot dance on several songs using Arduino.
I need to understand how to get the rhythm from the songs, something like detecting the beats, calculating the bpm..
I'm not expert but I read that probably I'll have to use some filters, is that right? or I can use a simple mic and opamp? do I have to analayze some frequencies?
I hope someone has some ideas and can explaing me how to do such a thing, thanks!! :)

You could use Sparkfun's spectrum analyzer shield. It will give seven bands measurements of the sound at the input jack. For BPM approximation you could select the band that suits you best for the music style you will be analyzing.


Necessity of Cyclic Prefix in OFDM I have a question on CP-OFDM.
In recent years people are talking about filter bank approach where they are trying to filter out the out of band side lobes using filters. Besides they are presuming that this approach may provide some robustness and spectral efficiency since they can get rid of CP. Is this true and what aspects we can remove CP and make spectral efficiency equal one.

Ref.
F. Schaich, T. Wild, Y. Chen , „Waveform contenders for 5G – suitability for short packet and low latency transmissions“, accepted for IEEE VTCs’14, Seoul, Korea, April 2014

In CP-OFDM, the transmitted signal is the output of the IDFT appended with a cyclic prefix. No pulse-shaping is used.

In so-called "new waveforms" proposed as candidate for 5G, pulse shaping is heavily used.

In FBMC, each subcarrier (i.e each complex exponential) is convoluted with a pulse with a certain shape (usually a bell shape with beginning and end at zero amplitude). As a consequence of convolution, the amplitude of each filtered subcarrier at the beginning and at the end of the multicarrier symbol goes to zero.
Consequently, in case of multi-path, this zero-valued beginnings and ends avoid ISI, since the part of the previous symbol overlapping most with the current one is its end, which as we said is zero-valued (or almost) and the current symbol has little energy at its beginning (zero or almost).

Hope this helps.


Pattern recognition for temporal data I'm trying to detect and classify non-speech sounds. Currently, I'm using a series of moving overlapped power spectrums from training sounds as the features I am looking for.

When I do analysis, I'm just computing the same amount of overlapped spectrums so that the number of features are the same. Right now the performance is not very good, it can only detect silence vs non-silence.

What techniques are there for this type of signal detection? One of my concerns is that for sounds of different lengths in the time domain would result in different lengths of feature vectors which so I can't use the same classifier, I'm stuck on this.

Are you trying to detect speech vs non-speech, or are there classes of non-speech sounds you are trying to discriminate from? I'm not clear from your question.

I think a decent first approach would be to block your signal into frames and compute the Mel-Frequency Cepstral Coefficients (MFCCs), as well as the delta-MFCCs (differences between adjacent frames' MFCCs) and delta-delta MFCCs (differences between MFCCs in frames that are two frames apart). This isn't the only way to do it, but without more specific knowledge of the problem domain, this is probably a good place to start.

Just googling should give you some good reference on how to compute the MFCCs if you are not already familiar with them. Basically you take the DFT, take the magnitudes, compute the energies inside triangular windows corresponding to human hearing, take the DCT of these coefficients, essentially as a compression step, and then discard the high order coefficients, usually only taking about the first twelve coefficients. I have an explanation of the meaning of the DCT step in this post: How do I interpret the DCT step in the MFCC extraction process?

You could then, say, use these coefficients as features for an SVM.


Adaptive Wiener Filter Coefficients Calculation I want to extrapolate a signal X of length 11, using the weiner filter coefficients W of length 7. The procedure I am using is as follows:

Compute the autocorrelation matrix upto lag 8 .
Using Levinson recursion, invert the autocorrelation matrix A.
Multiply the matrix of cross-correlations with the inv(A) to get W
Obtain Extrapolated Sample=Sum(W.*Signal)

My questions is that, if now I have a new signal Y of length 15 as:
Y=[X (4 new samples)]

So the Signal Y has first 11 samples of X and 4 new samples.

Now I want to compute weiner coefficients of this signal, given that I have already computed those for the signal X, can you recommend any possible approach to 'update' the filter co-efficients W given the 4 new samples, without having to do the whole weiner filter computation again thus saving some computational time?

Edit 1:

I have written the code for LMS Algorithm to iteratively 'adapt', when a new sample arrives, the filter co-efficients which were initially derived by Weiner Prediction.

The result is not same as the Weiner coefficients derived on the whole input.

The code is:

    `FilterOrder=7;

Window=7;

Filter=zeros(1,FilterOrder);

x = -20*%pi:0.1:20*%pi;

m=length(x);

StepSize=2.3e-4/1000000;

data=[cos(2*x(1:floor(m/2))) cos(3.3*x(floor(m/2)+1:m))];

Filter=WeinerCtrace(data(1:Window+1),FilterOrder);   // Initial Co-efficients

for i=1:50

Input=data(1+i:Window+i);               // Sliding Window of Filter Input

DesiredOutput=data(Window+1+i);         // New sample value

ActualOutput=Input(Window:-1:1)*Filter';

Error=DesiredOutput-ActualOutput;

while(1)

    NewFilter=Filter+StepSize*Error*Input;     // Gradient Descent

    if(sum((NewFilter-Filter).^2)<1e-20)

        break;

    end

    Filter=NewFilter;

end

Filter=NewFilter;

end
WCoefficients=WeinerCtrace(data(1:Window+i+1),FilterOrder)

Now, according to my opinion the values of WCoefficients and Filter should match as the WCoefficients are computed by applying Weiner Prediction on all the arrived data, while Filter is computed iteratively on introduction of every new sample.

Edit 2:
I have also implemented the RLS algorithm code given below, but the output of RLS still does not match the output of Weiner Prediction Coefficients for the same window:

FilterOrder=11;
Window=11;
t=0:1:1000;
f=linspace(.5,5,length(t));
fs=100;
data=cos((%pi/fs)*(f.*t));
[Filter,P]=WeinerCtrace(data(1:Window+1),FilterOrder);
w=1;
plot(data,'-r');
for i=1:500
    Input=data(Window+i:-1:1+i)';           // Sliding Window of Filter Input
    DesiredOutput=data(Window+1+i);         // New sample value
    ActualOutput=Input'*Filter;
    Error=DesiredOutput-ActualOutput;
    K=(P*Input)/(w+(Input'*P*Input));
    P=(P-(K*Input'*P))/w;
    Filter=Filter+(K*Error);
end
inp=data(i+1:Window+i+1);
WC=WeinerCtrace(inp,FilterOrder)
inp=data(i+1:Window+i);
for j=1:100
   NextValue=inp(11+j-1:-1:j)*WC;
   inp=[inp NextValue];
   plot(Window+1+i+j,NextValue,'r+')
end
inp=data(Window+i:-1:1+i);
for j=1:100
   NextValue=inp(1:Window)*Filter;
   inp=[NextValue inp];
   plot(Window+1+i+j,NextValue,'b*')
end

Can someone point out what am I doing wrong in the code or in anything else ?
Any suggestions/idea is welcome. Thanks!

The Wiener filter is, by definition

a) not adaptive and

b) not FIR / AR.

You seem to want an adaptive FIR filter.

There are many variants of this: LMS, NLMS, RLS (as you say), or the Kalman filter. Choose your poison! :-)

ALL of them are (generally) less computationally intensive that the Wiener filter.


Got stack in calculating state-space representation I got stack in the process of deriving a state-space representation of the following system:

There is an electrical oven described as follows:

control of the power supply $u$,

heating efficiency constant  $v$,

oven-product convection constant $p$,

oven-cover convection constant $c$,

cover-air convection constant $a$,

heater heating power $H$,

heating entering the product $P$,

heating entering the oven’s cover $C$,

heat loss to the surroundings $S$,

temperature of the product $y$,

temperature of the oven $T_o$,

temperature of the oven’s cover $T_c$,

temperature of the surrounding air $T_a$,

$$P = p(T_o−y),$$

$$C = \frac{c(T_o−T_c)}{10}, $$

$$S = a(T_c −T_a),$$

$$H = vu,$$

$$3\dot y  = P,$$

$$\dot T_o = H−P−C,$$

$$3\dot T_c = C−S$$

State variables are $y, T_o, T_c$,

Input is $u$.

On the left side I placed derivatives of state variables:

$$\dot y = \frac 13 pT_o - \frac 13 py$$

$$\dot T_o = uv - pT_o + py - \frac{c}{10}T_o + \frac{c}{10}T_c$$

$$\dot T_c = \frac{c}{30}T_o - \frac{c}{30}T_c - \frac a3T_c + \frac a3T_a$$

Let $x_1 = y, x_2=T_o, x_3=T_c$

$$\dot x_1 = - \frac 13 px_1 + \frac 13 px_2 $$

$$\dot x_2 = px_1 - (p+\frac{c}{10})x_2 + \frac{c}{10}x_3 + uv$$

$$\dot x_3 = \frac{c}{30}x_2 - (\frac{c}{30}+\frac a3)x_3 + \frac a3T_a$$

$$ \begin{bmatrix} \dot x_1 \\ \dot x_2 \\ \dot x_3 \end{bmatrix} =  \begin{bmatrix}
 -\frac 13p & \frac 13p & 0 \\
p & -(p+ \frac{c}{10}) & \frac{c}{10} \\
0 & \frac{c}{30} & -(\frac{c}{30}+\frac a3)
\end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} + \begin{bmatrix} 0 \\ v \\ ??? \end{bmatrix} u $$

The question is: what should be placed in ??? ?

Should it be $\frac a3T_a$? It is not connected with $u$.

And what to do with $\frac a3T_a$? Is it another input or should I erase it?

Model your deterministic input as a $2 \times 1$ matrix :

$$ \begin{bmatrix}
 u  \\
T_a \\
\end{bmatrix}
$$

Then your state space equations will be:

$$ \begin{bmatrix} \dot x_1 \\ \dot x_2 \\ \dot x_3 \end{bmatrix} =  \begin{bmatrix}
 -\frac 13p & \frac 13p & 0 \\
p & -(p+ \frac{c}{10}) & \frac{c}{10} \\
0 & \frac{c}{30} & -(\frac{c}{30}+\frac a3)
\end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} + \begin{bmatrix} 0 & 0 \\ v & 0 \\ 0 & a/3 \end{bmatrix} \begin{bmatrix} u \\ T_a\end{bmatrix} $$


Deriviation of the "Twiddle Sum" property I can't seem to understand how to derive the "twiddle sum" property:

$$\sum_{n=0}^{N-1}W_{N}^{kn}=N \ \delta[k\bmod N] $$
where $$ W_{N} \triangleq e^{\frac{j 2 \pi }{N}} $$
and $$ \delta[n] \triangleq  \begin{cases}
     1 & \text{if } n=0 \\
     0 & \text{otherwise}
     \end{cases}$$

I've tried doing the following:

$$\begin{align}
\sum_{n=0}^{N-1}W_{N}^{kn} & = \sum_{n=0}^{N-1}e^{\frac{j 2 \pi k n}{N}} \\
 & = \sum_{n=0}^{N-1} \left( e^{\frac{j 2 \pi k}{N}} \right)^n \\
 & = \frac{1 - \left( e^{\frac{j 2 \pi k}{N}} \right)^N}{1 - e^{\frac{j 2 \pi k}{N}}} \\
 & \\
 & = \frac{1 - e^{j 2 \pi k}}{1-e^{\frac{j 2 \pi k}{N}}} \\
 & \\
 & = \begin{cases}
    ?? & \text{if } k=mN \quad\quad m\in \mathbb{Z} \\
     0 & \text{otherwise} \quad k \in \mathbb{Z}
     \end{cases} \\
\end{align} $$

It seems I'm missing an $N$ before the delta, have I done a mistake ?

The most straightforward way to see this is to note that for $k=mN$

$$W_N^{kn}=e^{j2\pi mnN/N}=e^{j2\pi mn}=1$$

So the sum for the case $k=mN$ is simply

$$\sum_{n=0}^{N-1}1=1+1+\ldots+1=N$$

Note that the solution using L'Hopital's rule is a bit dubious because for $k=mN$ the formula for the geometric series is not valid because the terms are all equal to one (note the $r\neq 1$ condition in the link above). So you arrive at an expression that is only valid for $k\neq mN$, and then you somehow massage it into a solution that also works for $k=mN$, but this is bad math, even for engineers.


What is the difference between cubic interpolation and cubic "Spline" interpolation?. How to use it for upsampling purpose? After considering a couple of advices and suggestions for upsampling techniques here, I finally converged to use the cubic interpolation technique to estimate the voltage values corresponding to intermediate samples present between the original or previous samples. I know that spline interpolation is basically used for getting smoother curves, but what makes it different from the normal cubic interpolation technique as both of them use a 3rd degree polynomial to estimate intermediate values?.

Another implementation issue is that, for example, If I have some voltages corresponding to some samples say

    V = 3.674, 6.791, 8.888, 9.667.....
    Sample = 2, 3, 4, 5.....

Now, If we have to find the voltage information corresponding to the intermediate sample 3.5, then using cubic polynomial method, we arrive at 4 equations and 4 unknowns obtained by using the information provided by the neighboring samples closest to sample 3.5.

    V(2) = a + 2b + 4c + 8d
    V(3) = a + 3b + 9c + 16d
    V(4) = a + 4b + 16c + 64d
    V(5) = a + 5b + 125c + 125d

So, solving these equations I arrived at

   a = -4.428
   b = 4.3756
   c = -0.06299
   d = -0.04966

Using these values , we can calculate the voltage value at the sample 3.5 as

  V(3.5) = a + 3.5b + 12.25c + 42.875d
  V(3.5) = 7.186 Volts

Now my question - Is this method of interpolation suitable for large sampling rates?. How can I use this technique for upsampling a signal say from 10 sps to 100 sps for N = 1024 sample points?. I know that I have to develop a function that performs this cubic interpolation task between original samples (in C++), but I am just wondering how to implement it for upsampling for a continuous series of samples . Any suggestions, ideas or advices regarding the topic and its implementation would be appreciated. Thanks!

I understand that cubic interpolation can operate on 4 data points and the
more sophisticated technique I can think of is cubic spline. In case I am
using the normal cubic interpolation, how about I loop through the "N"
sample points i.e. 1024, for a condition below the "input sampling rate"
i.e. 10 sps considering 4 data points each and then performing the
interpolation function  based on the up sampling factor between each of
those 4 consecutive data points (Meaning - Interpolating/Estimating 10
values between each of those 4 data points) and then the function considers
the next 4 data points to perform the same operation and it goes on until a
100 samples has bee acquired i.e. the output sampling rate!

The difference between cubic interpolation as described in your question and cubic spline interpolation is that in cubic interpolation you use 4 data points to compute the polynomial. There are no constraints on the derivatives. Cubic spline interpolation computes a third order polynomial only from two data points with the additional constraint that the first and second derivative at the interpolation points are continuous. So if you have 4 points, then you compute 3 different polynomials (between points 1-2, 2-3, and 3-4), and these polynomials are smoothly connected in the sense that their first and second derivatives are equal at the given data points. This is not the case with standard polynomial interpolation.

But if you want an efficient implementation I wouldn't necessarily use polynomial or spline interpolation but I would go for a standard polyphase implementation of a low-pass interpolation filter. If $L$ is the upsampling factor you get $L$ parallel polyphase filters operating at the lower sampling rate, followed by $L$ upsamplers, the results of which are added. The low-pass filter can be designed using the window method (e.g. a Kaiser window). There is a lot of material available on the polyphase implementation of interpolators. This Master's thesis is very accessible. Just have a look at single-stage polyphase implementations using linear-phase FIR filters and ignore the rest.

EDIT:

If you want or need to use cubic interpolation, here's what you need to do. First of all, some coefficients in your equations are wrong, that's why also the result is not correct. These would be the correct equations:

    V(2) = a + 2b + 4c + 8d
    V(3) = a + 3b + 9c + 27d
    V(4) = a + 4b + 16c + 64d
    V(5) = a + 5b + 25c + 125d

However, you don't want to do it this way because if you always use the sample index as the $x$-coordinate you'll end up with a very large range for the coefficients of the system of linear equations (because you need to take the third power), and your coefficient matrix will sooner or later become ill-conditioned (asking for numerical problems). Second, and even more importantly, you would get a different system of equations for every new segment of the input signal, which would be very inefficient. You can do the interpolation by solving a $4\times 4$ system of linear equations with the same coefficient matrix each time. Just the left-hand side vector (containing the values of the given data points, V in your example) will be different for each segment.

You compute a cubic polynomial for each segment (i.e. for each range between two given data points) by considering the data points defining the segment and the two adjacent data points, just as in your example. If you define the cubic polynomial as

$$P(x)=a_0+a_1x+a_2x^2+a_3x^3$$

and if you define the $x$-axis values of the 4 current input data points as $-1$, $0$, $1$, and $2$, you get the following system of equations:

$$\begin{align}y_0&=a_0-a_1+a_2-a_3\\
y_1&=a_0\\
y_2&=a_0+a_1+a_2+a_3\\
y_3&=a_0+2a_1+4a_2+8a_3\end{align}$$

where $y_i$, $k=0,1,\ldots, 3$ are the current data points, and the current segment is between the points $y_1$ and $y_2$. The system above is readily solved as

$$\begin{align}
a_0&=y_1\\
a_1&=-\frac13 y_0 - \frac12 y_1 + y_2 - \frac16 y_3\\
a_2&=\frac12  (y_0 + y_2) - y_1\\
a_3&=\frac12  \left[(y_1 - y_2) + \frac13 (y_3 - y_0) \right]\\
\end{align}$$

So you just need to evaluate the above 4 expressions for the polynomial coefficients for each set of 4 data values. As soon as you have the coefficients, you can compute the interpolated output values in the current segment by evaluating the polynomial at $x_i=i/L$, $i=1,2,\ldots,L-1$, where $L$ is the upsampling factor ($L=10$ in your case). Then you forget the oldest input value, read the new input value, and you have a new set of 4 values from which you compute a new polynomial for the current segment, etc.

Here is a little Matlab/Octave script implementing this method:

x = randn(20,1);        % some input data
L = 10;                 % upsampling factor

% append points at both ends for first and last segment
N = length(x);
x2 = [x(1);x(:);x(N)];
Nseg = N-1;         % # segments

% matrix to compute polynomial coefficients from the current 4 data points: a = A*y
A = [0 1 0 0;-1/3 -.5 1 -1/6;.5 -1 .5 0;-1/6 .5 -.5 1/6];

% x-axis values of interpolated values in the current segment
mu = (0:L-1)/L;

% allocate output vector
Ny = Nseg*L + 1;
y = zeros(Ny,1);

% loop over all segments
for i = 1:Nseg,
    xt = x2(i:i+3);
    a = A*xt;
    yt = a(1) + a(2)*mu + a(3)*mu.^2 + a(4)*mu.^3;
    Iy = (i-1)*L + (1:L);
    y(Iy) = yt;
end

y(Ny) = x(N);

% plot result
tx = 0:N-1;
ty = (0:Ny-1)*(N-1)/(Ny-1);
plot(ty,y,'r-'), hold on
stem(ty,y,'r'), hold on
stem(tx,x), hold off


What is the definition of the bandwidth of a signal? I was reading over the famous Nyquist–Shannon sampling theorem and I was trying to understand what the term bandwidth meant rigorously (exactly) with no ambiguity. In this site the define it as follow:

  2) Bandwidth is the range of frequencies -- the difference between the highest-frequency signal component and the lowest-frequency signal component -- an electronic signal uses on a given transmission medium. Like the frequency of a signal, bandwidth is measured in hertz (cycles per second). This is the original meaning of bandwidth, although it is now used primarily in discussions about cellular networks and the spectrum of frequencies that operators license from various governments for use in mobile services.

which to me it means that given some periodic function:

$$ f(x) = \frac{1}{2}a_0 + \sum^{\infty}_{n=1} a_n cos(nx) + \sum^{\infty}_{n=1} b_n sin(nx) $$

we just need to extract the frequencies that are the largest and subtract them. For the sake of an example consider:

$$ f(x) = cos(2 \pi x) + cos(4 \pi x ) + sin( 3 \pi x) + sin( 7 \pi x) $$

the frequency of a sin/cos is $f = \frac{n}{2 \pi}$, thus, components of $f$ with the largest and smallest $n$ give the frequencies we want. So the band width $B$ of $f$ is:

$$ B = f_{max} - f_{min} $$

$$ B = \frac{n_{max} }{2 \pi} - \frac{n_{min}}{2 \pi} = \frac{7 }{2 \pi} - \frac{2}{2 \pi} $$

is this algorithm correct? Do I have the correct understanding of bandwitdh in the context of signals and systems?

I think less like an engineer and more like a mathematician so it was important for me to know exactly what bandwidth means in the context of Fourier series since its a totally unambiguous example to me (since all period functions can be expressed like that).

Also, as I was about to tag the question the definition of band width according to this site came up as:

  the difference between the upper and lower frequencies in a contiguous set of frequencies.

this reminded me that a lot of the time I see signals represented in what this community calls "frequency domain", as sometimes a continuous line (i.e. $f(c) = lim_{x \rightarrow c}f(x)$). This seems to be really unintuitive to me because the Fourier series is a countable summation but a continuous frequency domain plot would imply its a uncountable summation of the periodic function. In the light of this, how is "contiguous set of frequencies" a good definition? Or what does "contiguous" mean?

As you can see I just want to understand what bandwidth means in the end in a mathematical sense.

It's hard to give a definition of bandwidth without mentioning Fourier transforms.

I hope the following definition makes the greatest sense for a mathematician. Consider the space of Fourier transform functions $X(\omega)$ whose domain of support is a compact interval in $\omega$ (a single piece or two symmetric pieces rather than being scattered). Then we shall consider two cases:

case one is for real valued functions $x(t)$ whose Fourier transforms $X(\omega)$ are double-sided and conjuage symmetric about $\omega=0$ For such a function $x(t)$ the bandwidth is the length of the domain of support that resides in the positive $\omega$.

case two is for complex valued (analytic) functions $x(t)$ whose Fouerier transforms are single-sided and has zero value for $\omega<0$. For such a function $x(t)$ the bandwidth is the length of the domain of support that resides in the positive $\omega$.

Eventhough these two bandwidths are the same, the necessary Nyquist sampling rates are different; specifically it's twice larger for the real signals than the complex valued case.


Intermittent Harmonics, Is there a physical explanation? I am doing some frequency analysis of a vibration signal from a spinning rotor,

I am expecting to see peaks corresponding to the rotational frequency of the rotor. So the fundamental frequency should be ~0.7Hz, then I am assuming I should then see multiples of this corresponding to the harmonics?

I can clearly see peaks at ~ 1.4Hz, 2.8Hz, 3.5Hz, and 4.2Hz but nothing at 2.1Hz which I cant understand

I am not sure whether the explanation for this can be explained physically or if its related to my sensor setup / analysis methods? I dont think it is related to the frequency resolution since this is around 0.03 Hz which seems suitably low.

It seems likely to me that there is some kind of weighting to specific harmonics and the 2.1Hz harmonics is somehow damped but I can't think why this would be the case.

You have guessed it right. But before an explanation, you should make sure that your measuring setup is not the cause of this observation.

The fact that your fundamental frequency is 0.7 Hz. does not mean that you should see all harmonics of that frequency at 1.4 Hz., 2.1. Hz, 2.8 Hz. etc...

And furthermore, even if you would see those harmonics, their weights would most typically be different as suggested by the continuous-time Fourier series analysis of periodic waveforms.

The physical explanation is that your mechanical rotating system as a whole can be (approx) modeled as an LTI system with a self-frequency response composed of resonating or damping frequency modes...

And it seems that at 2.1 Hz., the system is not able to vibrate enough, some damping must be happening.

Nothing further can be stated at this point with the given information.


Filter signal for beamforming and phase response I was wondering this for quite some time but couldnt find a satisfied answer. Seems like not a lot worry about phase response in general. But maybe I'm thinking about it in a wrong way...

I want to apply a high pass filter to my microphone pressure signals to cut away the low frequency garbage before beamforming. So I was wondering about using filters (low, high, band) on your time signals before beamforming.

How do you guys correct for the phase response? And is it actually needed? If not, is only applying (e.g.) a butterworth filter adequate enough to perform correct beamforming afterwards (outside the frequency for which it is filtered)? Or is some correction in phase needed?

An ideal delay and sum beamformer is Linear, which implies that filtering before  beam forming is equivalent to filtering after beam-forming.   If you filter each sensor, the filters have to be matched very closely, so for a pure analog system, do you want to use one or N nearly identical filters.  Given temperature sensitivity and parameter drift, filtering after is much cheaper. No implemented beam former is perfectly Linear. Filtering prior has a number of advantages, such as reducing dynamic range requirements. Also for a digital beamformer, you need an antialiasing filter prior to the A/D converter so you filter before anyway.  30 years ago, switched capacitor filters were very useful. Today, for audio, over sampled delta sigma converters are not expensive and have huge dynamic ranges but you still have to do some signal conditioning between the mechanical sensor and the converter. The trade off for RADAR is more complicated. If your converters have enough dynamic range.

So the answer is that filtering prior to beamforming needs nearly identical characteristics. Some deviation is ok but you need to make that determination for your self . Ceramic and moving coils are not ideally Linear phase devices either and sensor location uncertainty has an effect.  You also need to consider that you might have nonacoustic paths between the source and your array,  Element coupling is another problem that can emerge.


EEG signal processing with sampling frequency What is the relation between sampling frequency  and the frequency of the signal after sampling? For example sampling frequency is 100 Hz,
that is the signal contain maximum frequency of 50 Hz. Then what will be the maximum frequency after sampling? Is it 100 or 50?
Or it is above it? How can we see it in MATLAB?

Intuitively... Assume to have a sampled signal, i.e. a sequence $\{s(n)\}$. The "fastest" oscillatory behaviour you can see is a component which changes "polarity" at any successive sample, for example consider
$$
s(n) = (-1)^n
$$

it is positive for $n=0$, negative for $n=1$, positive again for $n=2$ and so on. Trivially you cannot have something oscillating faster. The period of this signal is 2 (samples), infact you have $s(n) = s(n+2)$ for any $n$. Therefore its frequency is half the sampling frequency.


Developing a classifier for a periodic music I am trying to develop a classifier for Indian Taals, for which I would be focusing only on the Tabla.
Now once a extract a period of the Tala(a repetitive pattern), I want to develop some feature vectors for a Tala, which can be used to differentiate two Talas.

I am not sure what they should be. As of now, what I am planning to do is to find the beats of the Tala by beat tracking, or the onsets, and then take the median of the spectrogram values at the beat times to be my feature vector for a particular tala.  Also, there is the concept of Bayan/Dayan in Tabla which have different pitches. I am not very sure how to use this information into my feature vector.

I would also want to include the time(temporal characteristics) into my feature.

Any help would be appreciated.

Very first thing to do: there is literature describing systems doing exactly what you want to do. Look for it!

Your suggested approach - representing the sequence as a feature vector of dimension $n$ where $n$ is the number of onsets or beats - has many flaws:

What if you want to classify sequences which are longer/shorter than your training samples?
What if the onset detector skips an event or detects a spurious event, and the entire vector is shifted left or right?
What about timing information, which is essential in the perception of rhythm!
A single real number computed from the spectrogram is unlikely to represent well the musicological dimension that matters to your problem (which "stroke" or "bol" has been played on the tabla).

Strings/sequence matching methods are more likely to be what you need. What matters to your classification problem is "what follows what" rather than "what happens at the $n$-th position" - because the concept of "$n$-th" position is not very robust, and also very relative.

Here's how I would approach the problem:

Detect onsets on your audio signal, and for each onset, compute a vector of timbral features (MFCCs will be fine). At this point you have a sequence of onset times and a sequence of feature vectors.
On your entire training database, train a vector-quantizer with a small-ish codebook (say 16 or 32 vectors) on the feature vectors. There's a good chance that it will map well to the tabla bols ('dha', 'ti', 're', 'tin'...) - each codeword will correspond to a distinct tabla stroke.
Quantize your data using this vector quantizer. Your sequences will now be represented as a list of onset times and codeword id [(0.0, codeword_10), (0.51, codeword_7), (1.02, codeword_7), (1.27, codeword_16), (1.52, codeword_17), (2.03, codeword_10)...] etc.
Use a classic tactus/tatum detection algorithm (they use inter-onset times as their input) to recover the temporal structure of the rhythm. This will allow you to group events into matras. The example above will become: codeword_10, codeword_7, codeword_7, codeword_16|codeword_17, codeword_10 (note the grouping - the fourth and fifth events are played in sequence within the same matra).
Use any string/sequence modelling technique on the representation above:

Transforming the sequence obtained at step 4 into a sparse vector of n-gram counts is the most flexible method, since it allows you to use any text machine learning technique. The downside is that it needs lots of training data and works well on rather long sequences.
You can also use HMMs. Collect a small corpus of sequences in various taal, train a HMM on them for each taal; use the HMM to compute the likelihood of each of the taal models given the observed sequence.
Another option is to use an edit distance on the sequence obtained at step 4, in conjunction with nearest-neighbour classification. Compute the edit distance between the sequence to classify, and sequences representative of each of the different taals. Pick the taal with the smallest edit distance.

Note that if you have time to invest in manual annotation, steps 2 and 3 can be replaced by a supervised classifier: annotate a database of table sounds with the name or 'bol' of each stroke ('dha', 'dhin', 'ti', 're'...). Train a classifier from audio features (say MFCC) into bols. Use this classifier at step 3. You will now be able to work with a very meaningful representation of your sequences.

A final note: There's a good chance that step 5 can be done with a handful of rules, rather than through machine learning. I know the entire purpose of your project is to do tala recognition with machine-learning; but what is the point if a bit of musicological knowledge can do the same job? What machine learning can do is transcribe the audio signal into a symbolic representation (a sequence of bols) - but further processing of this representation to extract meaningful musicological properties/categories does not necessarily require machine learning techniques...


What is normalized angular frequency? I am new to DSP, and I am self-studying. I came upon this question, and I am stuck at it. Could someone please help me? Here is the question:

  In order to digitally create a sinusoid with frequency $f$ = 440 Hz, you
  will need the normalized angular frequency $\omega$. How large is $\omega$
  for a sampling frequency $f_s$ = 16 kHz?

HINT:

Determine the argument of a continuous-time sinusoid with frequency $f=440$ Hz.
Sample that continuous-time sinusoid by replacing $t$ by $nT$, $n\in\mathbb{Z}$, with $T=1/f_s$ being the interval between samples.

This gives you an expression $\sin(n\omega_0)$, where $\omega_0$ depends on $f$ and $f_s$. $\omega_0$ is the normalized angular frequency.


Compressing a time-domain signal in the frequency domain Matlab I'm new to frequency domain analysis. Let's say I have a frequency domain representation from a time domain signal. What I want to do is taking frequencies with a low amplitude out of this frequency domain representation and then convert this frequency domain representation back to the time domain such that the result is almost the original time domain representation.

In Matlab what I want to do is something like this:

timeDomain(1:50) -> apply fft(timeDomain) -> inspect frequency domain and filter out low amplitude frequencies yielding a frequency domain vector x with length(x) probably (way) smaller than 50. (so x should contain only the frequencies with the highest amplitudes of the signal) -> apply ifft on this filtered complex vector -> Get an approximation of the original time domain signal (timeDomain(1:50) that has also length 50.

So my aim is to represent the time domain signal timeDomain(1:50) with less dimensions in the frequency domain without too much information about the signal getting lost, thus being able to reconstruct the time domain signal, so basically compression.

Thanks for your time.

Everything depends on what you call compressing. What you are proposing is Fourier or spectrum thresholding. You will generate a $50$-length vector with some zeroes, hopefully. Let us ignore "complexness" that requires two floats for now,  you can deal with it using Hermitian symmetry.

Your main problem is actual compression. Forget about quantization first. The main issue is index location. Assume I have a transformed domain vector $T=[5, 1, 3 ,4, 0.5]$. Suppose I use a threshold on $2$. Now I give you the compressed version  $T_c=[5,  3 ,4]$, a ratio of $3/5$.

Can you retrieve $T_d=[5, 0, 3 ,4, 0]$? Probably not, you do not know where the zeroes were. This problem somehow prevented matching pursuit to be used for image compression, because coding the index penalized the final compressed file too much.


How do I detect distinct objects when their edges touch each other? I need to find all the contours in an image retrieved from the camera. So I first use the canny edge detector to find the edges and then find the contours. Pretty simple.

However, my contours get "merged". For example, in the image below, I clearly have 4 different objects. The edges are slightly touching at some points though, so I get one big contour instead of four separate ones. I tried changing around the thresholds, eroding, morphology operations and similar things but the edges keep slightly touching. Does anyone have any suggestions on how to get separate contours in images similar to the one below? (The image below is obviously just an example, my actual images are much more complex, but have the same basic problem).

Detecting different components:
If you're trying to detect the different components, there probably are other approaches to do them than detecting the contours. Here's an example in Mathematica. An erosion followed by dilation is used to close the gap in the second component before detection (if you don't do this, it won't detect it).
img = Binarize@Import["http://i.stack.imgur.com/yqDyu.png"];
Colorize[MorphologicalComponents[Dilation[Erosion[img,1],1]]]

The figure on the left below, shows imperfect object detection (without closing the gap) and on the right, shows the correct detection (running the above code).

Detecting the different contours:
However, if you indeed want to separate only the contours, here's an example. The erosion and dilation are performed as before to close the gap and the resulting image is run through a Canny edge detector. I've made the default options explicit, so that you can see what's being used.
img2 = EdgeDetect[Dilation[Erosion[img, 1], 1], Method -> "Canny"]

This will give you both the inside and the outside edge (see figure on the left below), since the pixel width is greater than 1 all around. I haven't had much luck trying to get it thinner, as the performance degrades (may be different for your other images). The inner contours are the ones you want, and the outer contour is just the combined contour of all 4 components. Now all we need to do is drop the outermost one with:
SelectComponents[img2, "EnclosingComponentCount", # > 0 &]

which gives you just the inner contours (see bottom right). In other words, it only picks those contours that are enclosed by at least one other contour, which automatically disqualifies the outermost one. I do not know the equivalent of these commands/operations in openCV.

Note that the apparent breaks in the figure are due to saving to jpeg in a smaller size. It does not look that way on my screen.


How to calculate wave length Given that 30GHz to 300GHz is considered mm wave spectrum which can be used for high-speed wireless communication such as 5G.
My question is how to calculate the wavelength if we know the frequency? I suppose we have to know the wavespeed as $ \lambda = v / f $
But do you know where to get wavespeed reference in telecommunications?

The "wavespeed" is the speed of light, which is 299,792,458 m/s if in free-space, or slowed by the root of the relative permittivity (dielectric constant for the material where the wave is propagating).  So $v = c/\sqrt{\epsilon_r}$ and the formula for wavelength is $\lambda = c/(f\sqrt{\epsilon_r})$, where $c$ is the speed of light (299,792,458 m/s), $f$ is the frequency of the signal in Hz, and $\epsilon_r$ is the dielectric constant.

Example dielectric constants:

Air at 0°C: 1.000594
Air at 20°C: 1.000528
Glass: 5 to 10
Paper: 3.85
Silicon: 11.68

For example, a 30 GHz signal in air would have a wavelength of nearly 10 mm: $ 299,792,458/30\times10^{9}$.


Can someone show the details of how to apply AIC for sinusoidal models to specific data? NOTE: This is a question that another user has been trying (unsuccessfully) to ask. Because the multiple questions asking, essentially, the same thing have either been deleted by me (because they were duplicates of each other) or by the user, I have tried to ask the question in such a way that will help the user and is a good fit for the site.

Suppose I have a vector of data measurements, $x(n)$, and I want to fit a model to that data:
$$
x(n) =  \sum_{p=1}^{P} a_p \cos(2\pi f_p + \phi_p) + \epsilon(n)
$$
where each of the $P$ sinusoids is parameterized by $(a_p, f_p, \phi_p)$ and $\epsilon$ is an i.i.d. Gaussian, white noise process with zero mean and variance $\sigma^2$.

My question is: How do we use AIC (or BIC or MDL) to estimate the model order, $P$ ?

The Wikipedia entry just says:
$$
AIC = 2P - 2\ln(L)
$$
where $L$ is the likelihood function.

How do I calculate $L$ given $x(n)$?

The PDF (likelihood function) of the multi-tone estimation problem is:
$$
{\cal L}({\bf x}; {\bf a}, {\bf f},{\bf \phi}, {\sigma}^2) = \frac{1}{(2\pi \sigma^2)^{\frac{N}{2}}}
\exp\left( - \frac{1}{2\sigma^2} \sum_{n=0}^{N-1} ( x(n) - \sum_{p=1}^{P} {a}_p \cos(2\pi {f}_p + {\phi}_p) )^2   \right)
$$

The paper by Djuric quotes the AIC and MDL values as:

The code at the end attempts to:

Generate a 5-tone noisy signal.
Do the ML estimates of the amplitudes, frequencies and phases for a given value of $p$.
Find the (log) likelihood values for each value of $p$.
Calculate the AIC and MDL values from Djuric.

The plot for $p = 1 \ldots 60$ is:

Usually, the AIC over-estimates the model order. In this particular case, AIC gets 7 and MDL gets 5 (the correct value). The minima of each curve are plotted (as a o on the AIC curve and as a + on the MDL curve).

Djuric, P.M., "A model selection rule for sinusoids in white Gaussian noise," Signal Processing, IEEE Transactions on  (Volume:44 ,  Issue: 7 ), pp. 1744 - 1751.

scilab CODE ONLY BELOW

P = 5;
f = [ 0.1 0.12 0.231 0.333 0.478634];
phi = [0.329847 5.90324 3.09834 4.907983 1.984];
a = [ 1 1 2 2 3];
sigma = 1;

N = 100;

x = cos(2*%pi*[0:N-1]'*f + ones(N,1)* phi )*a';

xn = x + sigma*rand(N,1,'normal');

function [L,LL] = likelihood(data,ahat,fhat,phihat,sigmahat)
    N = length(data);
    datahat = cos(2*%pi*[0:N-1]'*fhat + ones(N,1)* phihat )*ahat';
    L = 1/(2*%pi*sigmahat*sigmahat)^(N/2) * exp(-1/(2*sigmahat*sigmahat)*sum((data - datahat).^2));
    LL = - log(2*%pi*sigmahat*sigmahat)*N/2 - -1/(2*sigmahat*sigmahat)*sum((data - datahat).^2);
endfunction

function [a,m] = aic(p,L,N)
    a = 3*p - log(L);
    m = 3*p/2*log(N) - log(L);
endfunction

function [ahat,fhat,phihat,regr] = sin_estimate(data)
    PAD = 10;
    N = length(data);
    DATA = fft([data; zeros((PAD-1)*N,1)],-1);
    [mx,ix] = max(abs(DATA(1:N*PAD/2,1)));
    ahat = mx/(N/2);
    fhat = (ix-1)/N/PAD;
    phihat = atan(imag(DATA(ix)),real(DATA(ix)));

    regr = data - cos(2*%pi*[0:N-1]'*fhat + ones(N,1)* phihat )*ahat
endfunction

pValues = [1:60];
likeValues = [];
loglikeValues = [];
aicValues = [];
mdlValues = [];
for p=pValues,
    signal = xn;
    a_est = zeros(1,p);
    f_est = zeros(1,p);
    phi_est = zeros(1,p);
    for pp=1:p
        [a_est(pp),f_est(pp),phi_est(pp),signal] =  sin_estimate(signal);
    end

    [l,ll] = likelihood(xn,a_est,f_est,phi_est,sigma);

    likeValues = [likeValues l];
    loglikeValues = [loglikeValues ll];

    [aicv,mdlv] = aic(p,likeValues(p));
    aicValues = [aicValues; aicv]
    mdlValues = [mdlValues; mdlv]
end

clf;
plot(aicValues);
plot(mdlValues,'r');
title('AIC and MDL plots');


what method for thresholding? I want to do one threshold method on high pass components of a wavelet frame decomposition approch,for edge detection in image. I have 14 high pass components.

for thresholding, I need keep the lowest 50% of the coefficients in absolute value of the high pass parts or in other words keep only the lowest 50% largest coefficients.I examine otsu and Hysteresis thresholding, but I can't find all edges. I also wrote one  matlab code, that does for each high pass part:

% x1 and x2 are largest pixel value of each of high pass parts}
% there is negative and positive pixels in high pass parts

function Xt=Threshold(X,x1,x2)

[m,n]=size(X);
Xt=zeros(m,n);
for k=1:m
    for l=1:n
        if(X(k,l)> min(min(X)) && X(k,l)< x1 && mod(X(k,l),2)~=0)
            Xt(k,l)=1;
        elseif(X(k,l)> x2  && X(k,l)< max(max(X)) && mod(X(k,l),2)==0)
            Xt(k,l)=1;
        end
    end
end

I guess you had some misunderstandings.

First, strong edges correspond to wavelet coefficients with large magnitudes in high pass bands. Thus, instead of keeping small coefficients, you should keep large coefficients.

Secondly, in a high pass band, there should be a large number of coefficients with magnitude close to zero. This means edge pixels would almost impossible account for more than 50% of all pixels. Thus, the usage of 50% in your assumption is also doubtful.

Finally, the difficulty of edge detection often lies in differentiating weak edges from noise, but not the detection of strong edges.

I suggest to you to take look at the MATLAB source code of Canny Edge Detector, whose parameter adaption part ( when you do not provide the two canny thresholds to matlab ) might give you some thoughts on to automatically pick some good thresholds. Its idea is based on examining the edge histogram and determines threhsolds according to this distribution dynamically.


Smooth (not) equally-spaced signal strength data I'm currently trying to reduce the noise in/smooth sampled signal strength data. The context is that I collected data from a transmitting device using bluetooth - the available data is basically a series of timestamps and the receiving signal strength.

The problem: The signal strength fluctuates quite strongly, which is why I want to smooth my data. Actually the data is not equally spaced - some signals get lost, the transmitter has some quirks leading to irregular periods in which signals are sent etc. - but I could aggregate them on a fixed time-period, if that's necessary / recommended.

I read some articles but are still not sure how to approach this. What kind of filter would probably lead to the best results? Kalman, Savitzky-Golay, adaptive, exponential...many to choose from and not enough experience to do so.

What would you recommend?

Thanks for your time!

Between these I only have experience with Kalman filters. They don't need equally spaced measurements but you can make it a lot easier with it. Missing measurements are really easy to cope with.

If you want a simple solution and you think that the most likely outcome of your next measurement is the same as the last one (there are no 'dynamics') you could also interpolate your missing values and use a linear filter.


Where to begin? - Multi source polyphonic pitch detection I'm quite interested in real time pitch tracking. I have used mcleod pitch method, Yin and other methods to perform real time pitch tracking for a monophonic audio track.

I know that Polyphonic pitch detection is at a whole different level and I want to know where to start. If the experts(masters) could point me at the right direction, It would be greatly appreciated.
Papers, books, Topics or any suggestions that you could give me would be of great help.

It's a research topic. So try the music-ir.org MIREX web site for a vast number of academic papers and transcription competition entries.

If you want to try a newer direction, try gathering many 100's of hours of snippets of tagged music, fuzz it, and feeding that as training data to a DNN that includes convolutional+LSTM layers.


Question on Levinson's proposed discrete form of Wiener filter The whole foundation of Levinson's discrete version of Wiener filter is based on the assumption of stationarity of a time series, and aims to predict a value based on the past observed values. Now, if I try to use the theory to interpolate/predict values in a time series comprising financial data or other practical time series like the GDP of a country over some years, is it justified, since the financial data series is not necessarily stationary? Does the whole theory crumble due to this when applied for financial data, or do I still get adequate results, depending on certain properties of the data?

Any insight in this regard would be helpful.

Not specific to statistical Wiener filtering but also true in general, any algorithm will fail to work, if its assumptions are not met. However, it's also known that most algorithms will perform quite acceptable when their assumptions are sufficiently holding.

In this sense, if a statistical algorithm relies on the stationarity of its inputs, then it will be acceptably performing fine when the inputs are sufficiently stationary. Hence the keyword is sufficiency...

One might therefore argue that, a slowly changing (in statistical characteristics) non-stationary input would still allow the short term predictions possible , albeit with more erros than the ideal expectation of a strictly stationary input.

Taken the other way, many (if not all) nonstationary processes can be considered to be approximately stationary over a short enough period of time. Indeed this is the whole basis of applicability of most statistical signal processing algorithms that makes any use of real world data (which is, as you have outlined, nonstationary by its very nature, except the exceptional cases). Divide the input into pieces, assume stationarity over that period, apply the algorithm.

That being said for certain engineering signals of the type: speech, music, image, video, radar, sonar, seismology, you should make it clear by yourself that those financial data are also approximately stationary over some short periods of time. Otherwise you will have trouble in applying those algorithms that rely on stationary inputs to work.