Convoluzione di due segnali discreti in C

[simo85]

Un saluto a tutti

Tra poco comincerò il semestre di Segnali e Sistemi, quel poco che so l'ho studiato tempo fa da autodidatta..
Ci sono concetti che mi ricordo ed altri no, alcune cose che ho tralasciato e altre che penso di non aver appreso a pieno (forse la convoluzione è uno di questi) quindi sto ripassando e studiando in anticipo.

Sto ripassando sul libro Signals & Systems 2nd Edition by Oppenheim, Willsky e Nawab, e sono arrivato al capitolo 2 "Linear Time Invariant Systems", nel quale viene esposto un esempio di una convoluzione tra due segnali discreti:



La rappresentazione grafica dei segnali e della convoluzione è questa:

Sorgente FidoCadJ | Ingrandisci | Cos'è Fidocad

Dove:



Per chi possiede questo famoso (e stupendo a mio modesto parere) libro, l'esempio è il 2.1 a pag. 80
E quindi mi sono messo di mio conto ed ho voluto farlo in C.

Lo so che potrei andare su Google e scrivere "C convolution" o spulciare il codice sorgente di qualche audio player con effetti speciali, e questo sarebbe quanto basta a scopiazzare e non imparare, ma questo non è il mio obiettivo.

Ecco che ho scritto questo programmino, che è un po' fatto su misura per questo esempio, ma l'importante è rassicurarmi di aver capito bene l'implementazione.

Codice: Seleziona tutto

#include <stdio.h>

#define NX 2
#define NY 3

int main(void)
{
   int n = 0;
   int k = 0;

   double h[NY] = {1.0, 1.0, 1.0};
   double x[NX] = {0.5, 2.0};
   double y[NY + 1] = {0};   

   for(n = 0; n < NY + 1; n++)
   {   
      for(k = 0; k < NX; k++)
      {
         if(n >= k)
         {
            /* printf("x[%d]:%.2f\th[%d]:%.2f\ta:%.2f\n", k, x[k], n-k, h[n-k]); */
            y[n] += x[k]*h[n - k];
         }
      }

      printf("y[%d] = %.2f\n", n, y[n]);
   }
   
   return 0;
}

E questo è l'output, che coincide con l'esempio del libro.

Codice: Seleziona tutto

y[0] = 0.50
y[1] = 2.50
y[2] = 2.50
y[3] = 2.00


Ovviamente più che interessarmi il risultato dei valori ed esultare tutto contento, mi interessa essere sicuro di aver fatto bene, e quindi di aver capito bene.

Chiedo scusa per eventuali castronerie varie.
Ringrazio in anticipo.



Simo

[DarwinNE]

Io nel codice delle convoluzioni mi sbaglio sempre semprissimo

Di solito, uno dei tanti errori che faccio è non calcolare bene i limiti degli array che contengono i segnali. Nel tuo codice, non sarebbe meglio scegliere separatamente la lunghezza dell'array del segnale x e quello di h e poi ricavarsi la lunghezza dell'array che dovrà contenere y (se mi ricordo bene, dovrebbe essere legata alla somma delle due lunghezze)?

Prova semmai con un segnale x molto lungo e h con tre o quattro punti. Prova anche a fare calcolare al programma un filtro con una media mobile, cose del genere di cui puoi calcolare la soluzione in altro modo. Leggere file WAV non è difficile, potresti provare a filtrare un segnale e vedere cosa viene fuori con Audacity.

In certi casi, magari, potrebbe essere più efficiente passare per la FFT, ma quella che otterresti è una convoluzione circolare e non una convoluzione tout-court.

[simo85]

Ciao DarwinNE grazie per i consigli e gli spunti !

Per il momento sto provando ad implementare un esercizio successivo..
Questo era giusto per cominciare a recuperare...

Posso anche provare ad effettuare delle simulazioni con matlab in parallelo !
Per il momento è la via più veloce.

[simo85]

DarwinNE, per errore ho cancellato il messaggio precedente.

La dimensione è , e questo è un esempio (spartano) che ho scritto:

Codice: Seleziona tutto

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

void Conv(double x[], int sx, double h[], int sh)
{
   int n = 0;
   int k = 0;

   int sy = (sh + sx - 1);
   double *y = NULL;
   
   y = (double *) malloc(sy * sizeof(double));
   
   for(n = 0; n < sy;  n++)
   {
      y[n] = 0.0;
      for(k = 0; k < sx; k++)
      {   
         if((n >= k) && (n - k) < sh)
            y[n] += x[k]*h[n - k];
      }
         
      printf("y[%d] = %.2f\n", n, y[n]);
   }
   
   free(y);
}

int main(void)
{
   double h[5] = {0.5, 2, 2.5, 1};
   double x[6] = {1, 1, 1, 0, 0, 0};
   
   //double h[4] = {1.0, 2.0, 1.0, 2.0};
   //double x[8] = {2.0, 3.0, 4.0, 1.0, 9.0, 7.0, 8.0, 9.0};

   //double h[4] = {1.0, 2.0, 3.0, 4.0};
   //double x[10] = {1.2, 0.3, 5.7, 3, 2, 8.6, 3.2, 8.8, 3.2, 7.7};
   
   //double x[6] = {3.2, 4.5, 8.6, 2.1, 9.8, 2.2};
   //double h[4] = {2.2, 4.4, 1.2, 3.3};

   //double x[6] = {1.2, 2.5, 8.6, 2.4, 9.8, 2.5};
   //double h[4] = {2.2, 4.4, 1.2, 3.3};

   ///double x[8] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
   //double h[4] = {0.5, 2, 2.5, 1};

   Conv(x, 6, h, 5);
   return 0;
}


Il codice del messaggio # 1 è sbagliato perché non tiene conto della dimensione di e dei suoi indici..

Ho fatto varie prove coln matlab in parallelo e dovrebbe essere OK
Ho poi trovato anche questo semplice esempio grafico che è utile.

Convolution_example.pdf

(391.99 KiB) Scaricato 134 volte

Ciao,
Simo

[dimaios]

Esistono vari algoritmi per la convoluzione e vengono utilizzati in base alla natura dei segnali.
Qui una piccola rassegna non esaustiva ... solo per curiosità.
Se devi per esempio implementare l'algoritmo su un microcontrollore e non hai a disposizione la FFT oppure non cerchi la convoluzione ciclica puoi utilizzare altri algoritmi.

[simo85]

Grazie dimaios

[dimaios]

Di nulla. Ricorda che in questo campo bisogna scegliere il tool giusto per l'occasione giusta.
Infatti, soprattutto sui sistemi con scarse risorse di calcolo, esistono tecniche ad-hoc non molto diffuse che però funzionano per segnali con certe caratteristiche ( es. con pochi campioni oppure e brevi ecc. ).