LA CONVOLUZIONE 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 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 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
function
vet=rect1; tempo=20e-3; T=1/40000; campioni=(tempo/T)+1; k=(1:1:campioni); yf(length(k))=0; xf=(k-1)*T; for
i=1:round(campioni/2)-1
yf(i)=1; end 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'); zoom function
vet=rect2; tempo=20e-3; T=1/40000; campioni=(tempo/T)+1; k=(1:1:campioni); yf(length(k))=0; xf=(k-1)*T; for
i=1:campioni
yf(i)=1; end 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'); zoom
|