Realizziamo adesso un programma con il quale è possibile calcolare la convoluzione sia lineare che circolare tra due generici segnali introdotti dall’esterno. E’ possibile anche calcolare la convoluzione utilizzando i metodi : overlap & save e overlap & add.

 

conv.m

 

%                     LA CONVOLUZIONE

 

% Considereremo un segnale generico da introdurre dall'esterno

 

% così come la risposta impulsiva del filtro

 

clear all

 

err=0;

while err==0

   F=input('inserisci il nome della function contenente il segnale:','s');

   if exist(F)==0

      disp('LA FUNZIONE INTRODOTTA NON ESISTE');

   else

      err=1;

   end

end

err=0;

while err==0

   R=input('inserisci il nome della function contenente la risposta impulsiva del filtro:','s');

   if exist(R)==0

      disp('LA FUNZIONE INTRODOTTA NON ESISTE');

   else

      err=1;

   end

end

disp('Scegli il problema da risolvere : ');

err=0;

while err==0

   disp('1. Convoluzione lineare');

   disp('2. Convoluzione circolare');

   disp('3. Convoluzione circolare con l''uso della FFT');

   disp('4. Convoluzione utilizzando il metodo overlap & save');

   disp('5. Convoluzione utilizzando il metodo overlap & add');

   risp=input('Introduci il numero corrispondente : ');

   if risp~=1 & risp~=2 & risp~=3 & risp~=4 & risp~=5

      disp('Il valore introdotto é errato.Ricomincia');

   else

      err=1;

   end

end

f=feval(F);

r=feval(R);

N=length(f);

M=length(r);

if M<N

   c1=r;

   m1=f;

else

   c1=f;

   m1=r;

   pass=M;

   M=N;

   N=pass;

end

if risp==1

  

   % CONVOLUZIONE LINEARE

  

   c2=[zeros(1,N-1),c1,zeros(1,N-1)];

   for i=1:length(c2)

      c3(i)=c2(length(c2)-i+1);

   end

   for i=1:N+M-1

      y1(i)=sum(c3([N+M-i:2*N+M-i-1]).*m1);

   end

   figure

   plot(y1);

   title('convoluzione lineare');

   zoom

end

 

if risp==2

  

   %CONVOLUZIONE CIRCOLARE

   

   c2=[c1,zeros(1,N-M)];

   for i=1:length(c2)

      c3(i)=c2(length(c2)-i+1);

   end

   for i=1:N

      y2(i)=0;

      for j=1:N

         k=i-j+1;

         if k<=0

            k=k+N;

         end

         y2(i)=y2(i)+c3(N+1-j)*m1(k);

      end

    end

    figure

    plot(y2);

    title('convoluzione circolare');

    zoom

end

 

if risp==3

  

   %CONVOLUZIONE CIRCOLARE CON DFT

  

   c2=[c1,zeros(1,N-M)];

   ff1=fft(c2);

   ff2=fft(m1);

   ff=ff1.*ff2;

   y3=ifft(ff);

   figure

   plot(real(y3));

   title('convoluzione circolare con DFT');

   zoom

end

 

if risp==4

  

   % CONVOLUZIONE UTILIZZANDO IL METODO OVERLAP & SAVE

  

   a=N/M;

   if a>1

      m=ceil(2*(M-1)*log10(N/M));

      N1=M-1+m;

      cicli=fix(((N-N1)/m)+1);     % manca il calcolo per l'ultimo pezzo del vettore

   else

      m=0;

   end

   if m==0

      disp('Il metodo overlap & save non si può effettuare per questi segnali.');

      disp('Usare la convoluzione normale');

   else

      c2=[c1,zeros(1,N1-M)];

      C2=fft(c2);

      for i=1:cicli

         n=(i-1)*m+1;

         m0=m1(n:n+N1-1);

         M0=fft(m0);

         Y4i=M0.*C2;

         y4i=ifft(Y4i);

         if n==1

            y4(n:n+N1-1)=y4i(1:length(y4i));

         else

            n=n+M-1;

            y4(n:n+N1-M)=y4i(M:length(y4i));

         end

      end

      disp('Il vettore di dimensione maggiore è stato suddiviso in : ');

      if length(y4)==N

         cicli

      else

         cicli=cicli+1

         N1=N-length(y4)+M-1;

         m0=m1(N-N1+1:N);

         M0=fft(m0);

         c2=[c1,zeros(1,N1-M)];

         C2=fft(c2);

         Y4i=M0.*C2;

         y4i=ifft(Y4i);

         y4(N-N1+M:N)=y4i(M:length(y4i));

      end

      figure

      plot(real(y4));

      title('convoluzione con il metodo overlap & save');

      zoom

   end

end

 

if risp==5

  

   % CONVOLUZIONE UTILIZZANDO IL METODO OVERLAP & ADD

  

   a=N/(M-1);

   if a>=2

      m=fix(2*(M-1)*(log10(N/(M-1))-log10(2)));

      N1=M-1+m;

      cicli=fix(N/N1);      % manca il calcolo per l'ultimo pezzo del vettore

   else

      m=-1;

   end

   if m==-1

      disp('Il metodo overlap & add non si può effettuare per questi segnali.');

      disp('Usare la convoluzione normale');

   else

      c2=[c1,zeros(1,N1-1)];

      C2=fft(c2);

      y5j=zeros(1,M-1);

      for i=1:cicli

         n=(i-1)*N1+1;

         m0=[zeros(1,M-1),m1(n:n+N1-1)];

         M0=fft(m0);

         Y5i=M0.*C2;

         y5i=ifft(Y5i);

         y5(n:n+N1-1)=y5i(M:length(y5i));

         y5(length(y5)-N1+1:length(y5)-N1+M-1)=y5(length(y5)-N1+1:length(y5)-N1+M-1)+y5j;

         y5j=y5i(1:M-1);

      end

      disp('Il vettore di dimensione maggiore è stato suddiviso in : ');

      if length(y5)==N

         cicli

      else

         cicli=cicli+1

         N1=N-length(y5);

         m0=[zeros(1,M-1),m1(N-N1+1:N)];

         M0=fft(m0);

         c2=[c1,zeros(1,N1-1)];

         C2=fft(c2);

         Y5i=M0.*C2;

         y5i=ifft(Y5i);

         y5(N-N1+1:N)=y5i(M:length(y5i));

         if M-1>N1

            y5(N-N1+1:N)=y5(N-N1+1:N)+y5j(1:N1);

         else

            y5(N-N1+1:N-N1+M-1)=y5(N-N1+1:N-N1+M-1)+y5j;

         end

      end

      figure

      plot(real(y5));

      title('convoluzione con il metodo overlap & add');

      zoom

   end

end

segnale.m

function vet=segnale;

 

tempo=50e-3;

T=1/10000;

f1=200;

f2=420;

campioni=(tempo/T)+1;

k=(1:1:campioni);

xf=(k-1)*T;

yf=sin(2*pi*f1*xf)+sin(2*pi*f2*xf);

vet=yf;

figure

plot(xf,yf);

xlabel('sec');

title('grafico del segnale nel tempo');

ff=fft(yf);

zoom

figure

freq=[0:1/(T*campioni):(1/T)*(1-1/campioni)];

plot(freq,abs(ff));

xlabel('Hz');

title('spettro del segnale');

axis([0 500,min(abs(ff)) max(abs(ff))]);

zoom

banda.m

function r=banda;

 

f1=400;         %frequenza di taglio inferiore

f2=440;         %frequenza di taglio superiore

tempo=10e-3;

T=1/10000;

campioni=(tempo/T)+1;

k=(1:campioni);

x=(k-(length(k)+1)/2)*T;

for i=1:length(x)

   if i~=(length(k)+1)/2

      y(i)=(sin(2*pi*f2*x(i))-sin(2*pi*f1*x(i)))./(pi*x(i));

   end

end

y((length(k)+1)/2)=2*f2-2*f1;

r=y;

figure

plot(x*1000,y);

xlabel('ms');

title('grafico del segnale nel tempo');

zoom

ff=fft(r);

figure

freq=[0:1/(T*campioni):(1/T)*(1-1/campioni)];

plot(freq,abs(ff));

xlabel('Hz');

title('spettro del segnale');

axis([0 500,min(abs(ff)) max(abs(ff))]);

zoom