Embedded Music for low freq signals #1

PAk-CatchFire · 2019-11-05T23:24:17Z

Hello.
I am trying to implement a low memory Music algorithm to separate breathing and heart rate signals, which are 0.5Hz and 1-2Hz, on an ARM processor.

After some filtering I arrived to your implementation of Music, but I have seen that it is mostly prepared for 8kHz and 16 sample signals, and I am using 56 Hz as Fs and 200-500 sample points FFTs.

Would it be too difficult or memory dependant to adapt this library?
Thank you

spinlockirqsave · 2019-11-07T13:01:55Z

Hi, and many thanks for your question.

The libmusic will work with 56 Hz and 200-500 point FFT. The primary need addressed by this library was detection of DTMF frequencies in audio range, and this is why you can find lm_dtmf_detect function and other DTMF related things there, but library is meant to work with signals of other natures as well. The number of samples in a probe vector determines the accuracy of the detection, and depends on a number of frequencies (signal components) to be detected. It should be at least 8 times number of components, so for your 0.5 Hz and 1.2 Hz components it's at least 16 samples. The number of components should be known in advance for algorithm to give good results in the current implementation (limitation of MUSIC algorithm).

RAM

In your case of 2 signal components there will need to be at least this memory allocated(*):

Array 'A': 9 * 8 * sizeof(double) = 9 * 8 * 8 (on x86 and x86_64) = 576 bytes
Array A2: 9 * 8 * sizeof(double) = 9 * 8 * 8 (on x86 and x86_64) = 576 bytes
Array S: 9 * 8 * sizeof(double) = 9 * 8 * 8 (on x86 and x86_64) = 576 bytes
Array U: 9 * 9 * sizeof(double) = 9 * 9 * 8 (on x86 and x86_64) = 648 bytes
Array VT: 8 * 8 * sizeof(double) = 8 * 8 * 8 (on x86 and x86_64) = 512 bytes
Detections array: 0.5Fs * 10 = 40 kB
Re Array: 8 * 8 = 64 bytes
Im Array: 8 * 8 = 64 bytes
LAPACK's output: 4544 bytes
Other internal LAPACK memory: ? Test needed
Other internal structures = ~300 bytes
TOTAL RAM: (43.016 + 4.544 + ~0.3) kB + other LAPACK's internal mem = ~50.818 kB + other LAPACK's internal mem

(*) - on x86 and x86_64, but size of double may be different on your CPU/ALU

Library will reuse buffers (they are allocated once per all calls to lm_detect for a given detection session).

Please shout if you have any other questions.

PAk-CatchFire · 2019-11-07T23:57:36Z

Thank you!! That answer is quite helpful.
Right now we are working with a 57.6Hz sampling freq (Actually is a matrix of <4k points at 250kHz and then a sencond 2D FFT of 56.7Hz). The number of samples is not an issue, and we are using more than 200 hundred points for every MUSIC sliding window, we are searching for two frequencies.

What are the changes you would recommend to make the code work on this environment?

On the other hand, since I am trying to use it on a DSP, would it be possible to change the double variables to float?

Do you suggest an example for my test?

Thank you again

spinlockirqsave · 2019-11-08T12:46:08Z

Floats are possible. Which DSP you want to port this code for? Do you have SVD available on that processor?

PAk-CatchFire · 2019-11-08T13:42:10Z

Probably F28379D, or something in the family.

You can get a low cost board here:
http://www.ti.com/tool/LAUNCHXL-F28379D

BTW, why are you using the names: peak_697, etc?
As in:

if (f == 697) {
			d->info.peak_697 = peak;
		} else if (f == 770) {
			d->info.peak_770 = peak;
		} else if (f == 852) {
			d->info.peak_852 = peak;
		} else if (f == 941) {
			d->info.peak_941 = peak;
		}

```
if (f == 1209) {
		d->info.peak_1209 = peak;
	} else if (f == 1336) {
		d->info.peak_1336 = peak;
	} else if (f == 1477) {
		d->info.peak_1477 = peak;
	} else if (f == 1633) {
		d->info.peak_1633 = peak;
	}

``
Thank you again

PAk-CatchFire · 2019-11-08T13:46:23Z

Do you have SVD available on that processor?

You have some DSP and Math libraries, but probably not SVD per se. That's why I want to implement Music (or Burg) on it.

Both seem to work for this application (not for Yule Walker)

spinlockirqsave · 2019-11-13T13:19:09Z

Those peak_* variables are only for convenience of DTMF detection, they correspond to DTMF frequencies and are only used for quick access to frequencey detection results corresponding to DTMF.

Thank you for the link to the DSP board, I am having a look.

You have some DSP and Math libraries, but probably not SVD per se. That's why I want to implement Music (or Burg) on it

For MUSIC you need SVD to find eigenvectors corresponding to nullspace of correlation matrix.

PAk-CatchFire · 2019-11-13T14:32:21Z

For MUSIC you need SVD to find eigenvectors corresponding to nullspace of correlation matrix.

Where is the source code of your SVD function dgesvd?

Thank you

piotrgregor · 2019-11-13T14:35:59Z

dgesvd comes with LAPACK. It's definition can be found here

PAk-CatchFire · 2019-11-16T22:43:26Z

What are the minimum files I could use from the Laplace library to make it work on the DSP?
I am worried about memory size.
Regards

spinlockirqsave · 2019-11-21T18:01:29Z

Sorry, I am not sure what you mean by Laplace library? Did you mean LAPACK?

PAk-CatchFire · 2020-02-08T22:33:46Z

Sorry for the delay answering. Yes, I meant LAPACK....using a low memory solution would be the target.

PAk-CatchFire · 2020-02-09T11:18:05Z

Actually for SVD I was planning to use some other compilable library like:

https://github.com/lucasmaystre/svdlibc
or this:
https://github.com/kaushikb258/SVD_C

What do you think?
Is it compatible?

Regards

spinlockirqsave · 2020-02-15T14:27:30Z

Could be. Need to make sure you get the same results as in LAPACK tests (getting same eigenvectors is critical). Porting libmusic to use new implementation for SVD, having tests with LAPACK as a reference for new tests would be an option.

But since you care about memory footprint that part of main.c from https://github.com/lucasmaystre/svdlibc may signal that you will have more stuff to do than just making sure you get right vectors:

    if (transpose) {
      if (SVDVerbosity > 0) printf("  Transposing the matrix...\n");
      S = svdTransposeS(S); /* This leaks memory. */
    }
    svdWriteSparseMatrix(S, argv[optind], writeFormat);
  } else {
    DMat D = svdLoadDenseMatrix(optarg, readFormat);
    if (!D) fatalError("failed to read dense matrix");
    if (transpose) {
      if (SVDVerbosity > 0) printf("  Transposing the matrix...\n");
      D = svdTransposeD(D); /* This leaks memory. */

spinlockirqsave · 2020-02-15T14:30:07Z

As to the RAM requirements of libmusic they would be as follows:

RAM

In your case of 2 signal components there will need to be at least this memory allocated(*):

Array 'A': 9 * 8 * sizeof(double) = 9 * 8 * 8 (on x86 and x86_64) = 576 bytes

Array A2: 9 * 8 * sizeof(double) = 9 * 8 * 8 (on x86 and x86_64) = 576 bytes

Array S: 9 * 8 * sizeof(double) = 9 * 8 * 8 (on x86 and x86_64) = 576 bytes

Array U: 9 * 9 * sizeof(double) = 9 * 9 * 8 (on x86 and x86_64) = 648 bytes

Array VT: 8 * 8 * sizeof(double) = 8 * 8 * 8 (on x86 and x86_64) = 512 bytes

Detections array: 0.5Fs * 10 = 40 kB

Re Array: 8 * 8 = 64 bytes

Im Array: 8 * 8 = 64 bytes

LAPACK's output: 4544 bytes

Other internal LAPACK memory: ? Test needed

Other internal structures = ~300 bytes

TOTAL RAM: (43.016 + 4.544 + ~0.3) kB + other LAPACK's internal mem = ~50.818 kB + other LAPACK's internal mem

(*) - on x86 and x86_64, but size of double may be different on your CPU/ALU

Library will reuse buffers (they are allocated once per all calls to lm_detect for a given detection session).

Library has been tested against memory leaks, it is safe.

PAk-CatchFire · 2020-03-25T00:49:22Z

I am reviewing the library, but Detection array is too high.
Why is *10?

I would need to detect only two main frequencies (fundamental and first harmonic), using an array of 500 points of data (max window) and detecting frequencies betwee 0.8 and 4 Hz.

PAk-CatchFire · 2020-03-25T23:20:50Z

On the other hand, I would like to give it a test under Visual Studio, at least to test the results on real time on my environment.
Could you givev me a hand about porting it?
Reagrds

spinlockirqsave · 2020-04-07T07:05:28Z

Yes, please reach out to me at: piotr@dataandsignal.com

PAk-CatchFire · 2020-04-07T09:44:07Z

Done....email sent.

I am sharing my piece of code for Music here:

                                 //////////////////////////////////////
				//MUSIC
				//////////////////////////////////////
				cnt_music++;
				//if (cnt_music >= 200)//3seconds
				if (music_ahora)
				{
					cnt_music = 0;


					// Zero out Rxx prior to filling it using sger
					zeros(m, m, Rxx);

					// Create covariance matrix by doing outer product of v with
					// itself.
				
					cblas_sger(CblasRowMajor,    //Row-major storage //
						m,                // Row count for Rxx  //
						m,                // Col count for Rxx  /
						1.0f,             // scale factor to apply to v*v' //
						Fase_Latido_PC,
						1,                // stride between elements of v. //
						Fase_Latido_PC,
						1,                // stride between elements of v. //
						Rxx,
						m);               // leading dimension of matrix Rxx. //
										 
		  		//	printf("\nMatrix Rxx (%d x %d) =\n", m, m);
			   	//	print_matrix(Rxx, m, m);

					// Now compute SVD of Rxx.
					
					//info = sgesvd_( 'A', 'A',
					info = LAPACKE_sgesvd(LAPACK_ROW_MAJOR, 'A', 'A',
						m, m, Rxx,
						m, S, U, m,
						VT, m, superb);

					if (info != 0) {
						fprintf(stderr, "Error: dgesvd returned with a non-zero status (info = %d)\n", info);
						//return();
					}

				

        	//			printf("\nMatrix U (%d x %d) is:\n", m, m);
        	//			print_matrix(U, m, m);

        	//			printf("\nVector S (%d x %d) is:\n", m, 1);
        	//			print_matrix(S, m, 1);
        
        	//			printf("\nMatrix VT (%d x %d) is:\n", m, m);
        	//			print_matrix(VT, m, m);

				// Extract noise vectors here.  The noise vectors are held in Nu
				extract_noise_vectors(U, m, m, PSIG, Nu);
		
    		//printf("\nMatrix Nu (%d x %d) is:\n", m, NUMPTS-PSIG);
			  //	print_matrix(Nu, m, NUMPTS-PSIG);

				// Now find max freq.  Set up initial grid endpoints.  Freqs are
				// in units of Hz.
				fleft = 0.0f;
				fright = FSAMP / 2.0f;

				for (int j = 0; j<MAXRECURSIONS; j++) {
					 printf("j:%d fleft = %f, fright = %f\n", j, fleft, fright);

					// Set up search grid
					linspace(fleft, fright, NGRID, f);
					
				//	printf("\nVector f =\n");
				//	print_matrix(f, NGRID, 1);

					// Compute vector of amplitudes Pmu on grid.  music_sum wants normalized
					// frequencies 
					for (int i = 0; i < NGRID; i++) {
						Pmu[i] = music_sum(f[i] / FSAMP, Nu, NUMPTS, NUMPTS - PSIG);
					}
					printf("\nVector Pmu =\n");
					print_matrix_scale(Pmu, NGRID, 1,100000.0f);

					find_bracket(NGRID, Pmu, &ileft, &iright);
					fleft = f[ileft];
					fright = f[iright];
				}
				fpeak = (fleft + fright) / 2.0f;   // Assume peak is average of fleft and fright

				BPM_MUSIC = fpeak * 60;
				//printf("Peak frequency found at f = %f Hz\n", fpeak);
				}
				// usleep(500000);   // delay 1/2 sec.
/**/
				//////////////////////////////////////////////////
				// END MUSIC
			      //////////////////////////////////////////////////

PAk-CatchFire changed the title ~~Music~~ Embedded Music for low freq signals Nov 5, 2019

spinlockirqsave assigned spinlockirqsave and piotrgregor Nov 7, 2019

spinlockirqsave added question Asking for comments Feedback Request for feedback labels Nov 7, 2019

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Embedded Music for low freq signals #1

Embedded Music for low freq signals #1

PAk-CatchFire commented Nov 5, 2019 •

edited

Loading

spinlockirqsave commented Nov 7, 2019 •

edited

Loading

PAk-CatchFire commented Nov 7, 2019

spinlockirqsave commented Nov 8, 2019 •

edited

Loading

PAk-CatchFire commented Nov 8, 2019 •

edited

Loading

PAk-CatchFire commented Nov 8, 2019 •

edited

Loading

spinlockirqsave commented Nov 13, 2019

PAk-CatchFire commented Nov 13, 2019

piotrgregor commented Nov 13, 2019

PAk-CatchFire commented Nov 16, 2019

spinlockirqsave commented Nov 21, 2019

PAk-CatchFire commented Feb 8, 2020

PAk-CatchFire commented Feb 9, 2020

spinlockirqsave commented Feb 15, 2020

spinlockirqsave commented Feb 15, 2020 •

edited

Loading

PAk-CatchFire commented Mar 25, 2020

PAk-CatchFire commented Mar 25, 2020

spinlockirqsave commented Apr 7, 2020

PAk-CatchFire commented Apr 7, 2020 •

edited

Loading

Embedded Music for low freq signals #1

Embedded Music for low freq signals #1

Comments

PAk-CatchFire commented Nov 5, 2019 • edited Loading

spinlockirqsave commented Nov 7, 2019 • edited Loading

PAk-CatchFire commented Nov 7, 2019

spinlockirqsave commented Nov 8, 2019 • edited Loading

PAk-CatchFire commented Nov 8, 2019 • edited Loading

PAk-CatchFire commented Nov 8, 2019 • edited Loading

spinlockirqsave commented Nov 13, 2019

PAk-CatchFire commented Nov 13, 2019

piotrgregor commented Nov 13, 2019

PAk-CatchFire commented Nov 16, 2019

spinlockirqsave commented Nov 21, 2019

PAk-CatchFire commented Feb 8, 2020

PAk-CatchFire commented Feb 9, 2020

spinlockirqsave commented Feb 15, 2020

spinlockirqsave commented Feb 15, 2020 • edited Loading

PAk-CatchFire commented Mar 25, 2020

PAk-CatchFire commented Mar 25, 2020

spinlockirqsave commented Apr 7, 2020

PAk-CatchFire commented Apr 7, 2020 • edited Loading

PAk-CatchFire commented Nov 5, 2019 •

edited

Loading

spinlockirqsave commented Nov 7, 2019 •

edited

Loading

spinlockirqsave commented Nov 8, 2019 •

edited

Loading

PAk-CatchFire commented Nov 8, 2019 •

edited

Loading

PAk-CatchFire commented Nov 8, 2019 •

edited

Loading

spinlockirqsave commented Feb 15, 2020 •

edited

Loading

PAk-CatchFire commented Apr 7, 2020 •

edited

Loading