Power spectral density calculation coming negative Topic is solved

If you have a simple question and just want an answer.
Post Reply
shaileshwankhede
Experienced Member
Posts: 65
Joined: Fri Dec 02, 2016 1:30 pm

Power spectral density calculation coming negative

Post by shaileshwankhede »

Hi,
I am having xCore array microphone kit with 7 mics and using lib_dsp api for calculation of power spectral density of voice signal arriving at each mic.
In one thread I am continuously storing buffer and once half of the new buffer data is available, I am applying algorithm to calculate power spectral density.
Since I want to later perform cross-correlation between adjacent mic data, I am using split spectrum technique to save time and memory.
Keeping in mind, fourier transform of auto-correlation function is Power spectral density, I am calculating auto-correlation of mic Data. But when I trace or analyze this data, I am getting some negative values in between.

How can PSD be negative? Is there any calculation problem? However when I observe particular frequency region of PSD, values are raised properly when speak.
Shall I ignore these negative spikes or do I need to take absolute values of auto-correlation?
Half buffer data size is 500 samples and number of FFT points is 1024. Even I am zero padding signal properly. I have attached PSD trace (in purple color) for some of the frequency bins (5 to 130). (Sampling rate 16KHz).
Below is code snippet for calculation. Please let me know what's going wrong.

Code: Select all

#define BUFFER_SIZE           1000
#define HALF_BUFFER           500

#define N_FFT_POINTS          1024

/*Below part in C file for use in interfaces*/
dsp_complex_t audioBuffer[3][BUFFER_SIZE];						  //This will store signal in pair. mic1 & mic2; 3&4; 5&6
void fillMicData(mic_array_frame_time_domain* audio, uint16_t wrIdx)
{
    /*Ignore centre mic 0 data*/
    for(unsigned i=0; i<3; i++) {
        audioBuffer[i][wrIdx].re = audio->data[2*i + 1][0];       //1st signal
        audioBuffer[i][wrIdx].im = audio->data[2*i + 2][0];       //2nd signal
    }
}

void getMicData(dsp_complex_t micBuff[N_FFT_POINTS], uint8_t buffIndex, uint16_t rdIdx)
{
    dsp_complex_t zero_padded_data[N_FFT_POINTS - HALF_BUFFER] = {{0}};
    memcpy(micBuff, (dsp_complex_t*)&audioBuffer[buffIndex][rdIdx], sizeof(dsp_complex_t)*HALF_BUFFER);
    memcpy((dsp_complex_t*)&micBuff[HALF_BUFFER], zero_padded_data, sizeof(dsp_complex_t)*(N_FFT_POINTS - HALF_BUFFER));            //padding zeros to make length of N_FFT_POINTS.
}
/*C file part ends*/

/*In xc file. This buffer Data is called continously on every sampling of mic data*/
unsafe void bufferData(mic_array_frame_time_domain * unsafe audio, server interface algorithm_if algo_if)
{
    static uint16_t wrIdx = 0;
    static unsigned buffReadPosition = 0;

    select {
        case algo_if.getBufferReadPosition(unsigned& pos) :
            pos = buffReadPosition;
            break;
        default:break;
    }

    fillMicData(audio, wrIdx);			//Filling of mic buffer happens here

    if(wrIdx < BUFFER_SIZE-1)
        wrIdx++;
    else {
        wrIdx = 0;
        buffReadPosition = BUFFER_SIZE >> 1;
        algo_if.newBufferAvailable();	//triggers calculation of PSD
    }
    if(wrIdx == BUFFER_SIZE/2) {
        buffReadPosition = 0;
        algo_if.newBufferAvailable();	//triggers calculation of PSD
    }
}

int32_t temp_re[2][N_FFT_POINTS/2], temp_im[2][N_FFT_POINTS/2];
int32_t result_re[N_FFT_POINTS/2], result_im[N_FFT_POINTS/2];
#define Q_FORMAT 24
void soundLocalizationAndBeamforming(client interface algorithm_if algo_if)
{
	while(1) {
        select {
            case algo_if.newBufferAvailable() :
                {
                    dsp_complex_t micChData[N_FFT_POINTS];
                    unsigned rdIdx;
                    algo_if.getBufferReadPosition(rdIdx);				//get which half part of circular buffer to be read
                    for(unsigned i=0; i<3; i++){
                        getMicData(micChData, i, rdIdx);
                        dsp_fft_bit_reverse(micChData, N_FFT_POINTS);
                        dsp_fft_forward(micChData, N_FFT_POINTS, FFT_SINE(N_FFT_POINTS));

                        dsp_fft_split_spectrum(micChData, N_FFT_POINTS);
                        dsp_complex_t (*w)[2] = (dsp_complex_t(*)[2])micChData;				//w[0][N_FFT_POINTS/2] -> 1st signal FTT, w[1][N_FFT_POINTS/2] -> 2nd signal FFT

                        /*Separate real and imaginary data and Calculate complex conjugate of 2nd*/
                        for(unsigned j=0; j<(N_FFT_POINTS/2); j++) {
                            temp_re[0][j] = (*w)[0].re;
                            temp_im[0][j] = (*w)[0].im;
                            temp_re[1][j] = temp_re[0][j];				
                            temp_re[1][j] = -temp_im[0][j];			//Complex conjugate
                            w++;
                        }
                        /*PSD calculation*/
                        dsp_vector_mulv_complex(temp_re[0], temp_im[0], temp_re[1], temp_im[1], result_re, result_im, N_FFT_POINTS/2, Q_FORMAT);    //cross correlation
						
						/*result_re should have auto-correlation ie. PSD of micChData
						Why some of its values negative? What is significance of result_im?*/
                }
                break;
            default:
                break;
        }
    }
}

Thanks,
Shailesh
Image


View Solution
User avatar
johned
XCore Addict
Posts: 185
Joined: Tue Mar 26, 2013 12:10 pm
Contact:

Post by johned »

Absolutely nothing wrong with -ve values.

Feel free to try it in Python :

Code: Select all

import numpy as np
import pylab as pl
rate = 8000.0
t = np.arange(0, 10, 1/rate)
x = np.sin(2*np.pi*500*t) + np.sin(2*np.pi*1100*t) + np.random.randn(len(t))*0.2
acor = np.correlate(x, x, "same")
p = np.fft.rfft(acor)
f = np.linspace(0, rate/2, len(p))
pl.plot(f, p)
Are you thinking about dB ? Even then, all the numbers will be negative if you scale the peak to 0dBFS.

John
shaileshwankhede
Experienced Member
Posts: 65
Joined: Fri Dec 02, 2016 1:30 pm

Post by shaileshwankhede »

Hi John,

So to calculate total power in spectra of interest can I use Absolute vector sum of Real result and thus ignoring imaginary part?

Thanks,
Shailesh
Post Reply