Costruiamo adesso una spline bicubica mettendo in evidenza l’effetto del ricampionamento di un fattore non intero e dell’interpolazione su dati non sufficientemente sovracampionati come fatto per l’interpolatore lineare e la spline.

 

bicubica.m

 

% COSTRUZIONE DI UNA SPLINE BICUBICA

 

% Considereremo un segnale generico da introdurre dall'esterno

 

% può essere inserito anche un fattore non intero di ricampionamento

% e il periodo di campionamento può essere generico

 

clear all

tempo=input('inserisci l''istante finale di rappresentazione del segnale (in sec) : ');

f1=input('inserisci la frequenza di campionamento (in Hz) : ');

T1=1/f1;

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

T=1/(100*f1);

k_=(tempo/T)+1;

k=(1:1:k_);

x=(k-1)*T;

y=feval(F,x);

x=(k-1)*T*1000;

figure

plot(x,y,'y');

xlabel('ms');

hold on

 

T=T1;

k_=(tempo/T)+1;

k=(1:1:k_);

x=(k-1)*T;

y=feval(F,x);

x=(k-1)*T*1000;

if length(x)<40

   plot(x,y,'r*');

else

   plot(x,y,'r');

end

 

disp('OPERAZIONE DI RICAMPIONAMENTO ');

valore=input('introduci il fattore di ricampionamento : ');

a=input('inserire il valore di a : ');

n=length(y);

valori=valore*(n-1)+1;

T1=x(n)-x(1);

T2=T1/(valori-1);

xi(1)=x(1);

for i=1:valori-1

   xi(i+1)=xi(i)+T2;

end

m=length(xi);

for j=1:m

   yi(j)=0;

   for i=1:n-2

      if xi(j)>=x(i) & xi(j)<=x(i+1)

         xin=(xi(j)-x(i))*f1/1000;

         yi(j)=yi(j)+(1-(1+a)*xin^2+a*xin^3)*y(i);

      end

      if xi(j)>=x(i+1) & xi(j)<=x(i+2)

         xin=(xi(j)-x(i))*f1/1000;

         yi(j)=yi(j)+(a-2)*(xin^3-5*xin^2+8*xin-4)*y(i);

      end

   end

   for i=3:n

      if xi(j)>x(i-1) & xi(j)<x(i)

         xin=-(xi(j)-x(i))*f1/1000;

         yi(j)=yi(j)+(1-(1+a)*xin^2+a*xin^3)*y(i);

      end

      if xi(j)>x(i-2) & xi(j)<x(i-1)

         xin=-(xi(j)-x(i))*f1/1000;

         yi(j)=yi(j)+(a-2)*(xin^3-5*xin^2+8*xin-4)*y(i);

      end

   end

   if xi(j)>x(1) & xi(j)<x(2)

      xin=-(xi(j)-x(2))*f1/1000;

      yi(j)=yi(j)+(1-(1+a)*xin^2+a*xin^3)*y(2);

   end

   if xi(j)>x(n-1) & xi(j)<x(n)

      xin=(xi(j)-x(n-1))*f1/1000;

      yi(j)=yi(j)+(1-(1+a)*xin^2+a*xin^3)*y(n-1);

   end

end

yi(m)=y(n);

if m<55

   plot(xi,yi,'b+');

else

   plot(xi,yi,'b');

end

hold off

zoom

ff=fft(yi);

freq=[0:(f1*valore)/m:f1*valore*(1-1/m)];

figure

plot(freq,abs(ff))

xlabel('Hz');

zoom

 

 segnale.m

 

function f=segnale(x);

f=sin(2*pi*200*x)+sin(2*pi*350*x+pi/7)+sin(2*pi*420*x-pi/4);     % la frequenza è in Hz

 

impulso.m

 

function i=impulso(x)

n=length(x);

i(n)=0;

i(round(n/2))=1;

 

rect1.m

 

function r=rect1(x)

n=length(x);

r(n)=0;

for i=round(2*n/5) : round(3*n/5)

   r(i)=1;

end

 

rect2.m

 

function r=rect2(x)

n=length(x);

r(n)=0;

for i=round(49*n/100) : round(51*n/100)

   r(i)=1;

end