Спектральный анализ дискретных сигналов

Уральский технический институт связи и информатики (филиал)

Сибирского государственного университета

Телекоммуникаций и информатики

(УрТИСИ ГОУ ВПО «СибГУТИ»)

 

 

 

 

 

 

Курсовая работа.

По дисциплине

«Языки программирования».

Тема: Спектральный анализ дискретных сигналов

 

 

 

 

 

 

 

Выполнил: Астраханский Д.А.

Группа: МЕ-51

 

 

Екатеринбург, 2008

Задание для курсовой работы.

Написать программу на языке программирования Паскаль для решения следующей задачи (вариант задания индивидуальный). Результаты расчетов должны выводиться на экран и в файл. Оформление графиков и таблиц выполнять средствами математических пакетов (Maple, MathCad). Демонстрационный вариант программы подготовить в среде визуального программирования Delphi.

Проверить решение промежуточных задач средствами математических пакетов. Построить блок-схемы задачи и вспомогательных частей алгоритма. Оценить погрешность выполненных расчетов.

Напряжение U=U(t) на входе транзистора как функция времени описывается дифференциальным уравнением

с начальными условиями (1), где n – последняя цифра номера зачетной книжки, k – коэффициент усиления (см. ниже), fs(t) – периодический сигнал

(рис. 1), m – коэффициент обратной связи.

 

Указания и пояснения.

1. Дифференциальное уравнение с заданными начальными условиями (задача Коши) решается методом Рунге-Кутта второго порядка с коррекцией (3) на отрезке [0;5] с шагом  h=0.01.(в узлах  tj =jh, j=0,1,2…). Функция fs(t) в правой части представляет собой регулирующий периодический (период Т) сигнал единичной амплитуды (рис 1, номер варианта  n – последняя цифра номера зачетной книжки,). Результаты расчетов—таблица (tj,Uj) и график функции U(t) (на экран и в файл).
2. Значение коэффициента усиления k в правой части дифференциального уравнения есть наименьший положительный корень полинома (2), который вычисляется одним из методов нахождения корней уравнения (метод касательных, метод простой итерации).
3. Построить спектральные характеристики периодического сигнала fs(t), заданного в аналоговой форме и в виде дискретного сигнала. Длительность сигнала равна 1, период T=k.
4. Период функции U(t) определить с помощью функции автокорреляции.

 

Курсовая работа выполняется в ЧЕТЫРЕ  этапа.

1. Средствами математического пакета Maple  решается задача спектрального анализа аналогового и дискретного периодического сигнала fs(t). Сравниваются спектры амплитуд аналогового и дискретного представлений сигнала. (образец выполнения задания – файл вариант11.mws).
2. Создается проект в визуальной среде Delphi, решающий эту же задачу для дискретного сигнала, а результаты выполнения сравниваются визуально.
3. С помощью языка программирования системы Maple решается задача интегрирования дифференциального уравнения (задача Коши) методом (по варианту задания). Окончательные вычисления в программе зависят от результатов расчета программы в Delphi (следующий пункт). Образец выполнения задания – файл RUTTA.mws.
4. Создается проект в визуальной среде Delphi, решающий ту же задачу Коши, результаты расчета которой записываются в файл,  который используется в предыдущем пункте. Выводятся графики результатов вычислений в Maple и Delphi и сравниваются между собой. Явные несовпадения свидетельствуют об ошибке в программе на  Delphi.

Оформление:

• Формат А4.
• Титул
• Постановка задачи
• Алгоритмы решения вспомогательных задач
• Блок-схемы
• Результаты расчетов, графики
• Литература

Индивидуальное задание № 18

1. Начальные условия: U(0)=0.2
2. полином:x^4-x^3-3
3. коррекция:по средней производной
4. метод:итерации  
   

Часть 1.

Решение задачи спектрального анализа (вариант 02 )

> restart;

> with(linalg):with(plots):with(plottools):

> fun:= proc(t) local z ;z:=piecewise(t<0,0,t<1/2,2*t,t<1,1,0); evalf(z);end;

> fun(t/tau);

> p(x):=x^4-x^3-3;

> Koeff:=fsolve(p(x)=0,x,0..3);

> tau:=Koeff:

> plot(p(x),x=Koeff-0.5..Koeff+0.5,thickness=2,color=black);

> R1:=plot(fun(t),t=0..2.5,thickness=2,linestyle=3,color=blue):

> R11:=plot(fun(t/tau),t=0..2.5,thickness=2,color=black):

> display(R1,R11);

> Fourier_T:=proc(F,T0,TN,k::evaln) local T;

  global A0,Ak,Bk;

    T:=TN-T0;

   A0:=2/T*Int(F(x),x=T0..TN);

   Ak:=2/T*int(F(x)*cos(k*x*2*Pi/T),x=T0..TN):

   Bk:=2/T*int(F(x)*sin(k*x*2*Pi/T),x=T0..TN):

end proc:

> Trig_polynom:=proc(N,T,a0,ak,bk,k::evaln) local n,Pol,A0,A,B;

  global a,b,RisTrig;

  a:=array(0..N);b:=array(0..N);

    A0:=evalf(a0);a[0]:=A0;b[0]:=0;

    A:=seq(evalf(subs(k=n,ak)),n=1..N);

    B:=seq(evalf(subs(k=n,bk)),n=1..N);

     for n from 1 to N do

      a[n]:=A[n];b[n]:=B[n];

     end do;

    Pol:=A0/2+sum(A[k]*cos(2*Pi*k*x/T)+B[k]*sin(2*Pi*k*x/T),k=1..N):

    RisTrig:=plot(Pol,x=-T/2..3*T/2,color=blue,thickness=2):

  RETURN(Pol);

end proc:

> ARR:=proc(n::integer,c) local L,H,ma,mi,k::integer,

   Sim::array(0..n);

   ma:=c[0];mi:=c[0];

   L:=line([0,c[0]],[n,c[n]],thickness=2,color=red);

  for k from 1 to n do

   if c[k]>ma then ma:=c[k];end if;

   if c[k]<mi then mi:=c[k];end if;

  end do;

  H:=ma-mi;

  if H=0 then RETURN(L) end;

  for k from 0 to n do

   if abs(c[k])<H/1000 then

     Sim[k]:=ellipse([k,c[k]],0.2,0.01*H,color=blue);

   else

    Sim[k]:=plottools[arrow]([k,0],[k,c[k]],0.2,0.2,0,color=black);

   end if;

  end do;

  convert(Sim,list);

end:

> Spectr:=proc(n,a,b,c,Risphi) local k,R,phi;

   for k  from 0 to n do

     c[k]:=evalf(abs(I*a[k]+b[k])):

#    print(k,c[k]);

     phi:=evalf(argument(I*a[k]+b[k]));

     R[k]:=[eval(k),eval(phi)];

   end:;

Risphi:=plot(convert(R,list)):

end:

 

> T:=3;# величина периода

> F_for_all:=proc(t) global tau;fun(t/tau);end proc:;

> Ris1:=plot(F_for_all(t),t=0..T,color=brown,thickness=2,discont=true):display(Ris1);

> Fourier_T(F_for_all,0,T,k):

> a0:=evalf(A0);

> Nk:=80;

> Trig_polynom(Nk,T,A0,Ak,Bk,k):

> display(RisTrig,Ris1);

> Spectr(Nk,a,b,c,'Risphi1');

> display(ARR(Nk,c));

>

> Ampl:=display(ARR(Nk,c)):;

> 2: DTF:=proc (y,N,Y) local n,k,j,p,h;

n:=N-1;

h:=2*Pi/N;

2.1: for k from 0 to N do

p:=0;

   for j from 0 to n do

     p:=p+evalf(y[j]*exp(-I*k*j*h));

   end;

  Y[k]:=evalf(1/N*p);

end:

end:;

> 3: CDTF:=proc(N,Y,y) local n,k,h,p,j;

n:=N-1;

  h:=2*Pi/N;

for k from 0 to n do

p:=0;

   for j from 0 to n do

     p:=p+Y[j]*exp(I*k*j*h);

   end;

  y[k]:=evalf(Re(p));

end:

end:

> Setka_DTF:=proc(Nt,T,F,Y::array) local h,j,x,R,RL; 

   global GrafF;

   h:=T/Nt;

  for j from 0 to Nt do

   x:= evalf(j*h);

   Y[j]:= F(x);

   R[j]:=[j,eval(Y[j])];

  end:

    5.1: R[Nt]:=[j,eval(Y[0])];

     GrafF:=plot(convert(R,list),0..Nt-1,color=brown,

     style=point,symbol=circle):

end:

> Spectr_DTF:=proc(n,C,A,phi) local k,R;global Risphi;

  6.1:  for k  from 0 to n do

     A[k]:=evalf(abs(C[k])):

     phi[k]:=evalf(argument(C[k]));

     R[k]:=[eval(k),eval(phi[k])];

end:;

Risphi:=plot(convert(R,list),thickness=2,color=blue,style=point,symbol=box):

end:

Параметры задачи

> Nt:=80:`число дискретных отсчетов `:

> n:=Nt;N:=Nt-1;# параметры ДПФ

> C:=array(0..n):phi:=array(0..n):A:=array(0..n):;

Y:=array(0..N):

> Setka_DTF(N,T,F_for_all,Y);

 

> DTF(Y,Nt,C,n):

> Spectr_DTF(n,C,A,phi):

Для четных N

> display(ARR(n,A));

> display(ARR((n-1)/2,A));

> CDTF(Nt,C,F):

> display(GrafF,ARR(n-1,F));

> Setka:=proc(Nt,T,F,Y::array) local h,j,x,R,RL; 

   global GrafF;

   h:=T/Nt;

  for j from 0 to Nt do

   x:= evalf(j*h);

   Y[j]:= F(x);

   R[j]:=[x,eval(Y[j])];

  end:

    5.1: R[Nt]:=[x,eval(Y[0])];

     GrafF:=plot(convert(R,list),0..T,color=brown,

     style=point,symbol=circle):

end:

> F_Discret:=proc (Y,N,a,b,n) local k,j,p,q,h;

  h:=2*Pi/N;

for k from 0 to n do

p:=0;q:=0;

   for j from 0 to N do

     p:=p+evalf(Y[j]*cos(k*j*h));

     q:=q+evalf(Y[j]*sin(k*j*h));

   end;

  a[k]:=2/N*p;b[k]:=2/N*q;

# print(k,a[k],b[k]);

end:

if 2*n=N then b[n]:=0; end;

RETURN(n);

end:;

> 3: Trig:=proc(t,n,T,a,b) local z,k;

z:=a[0]/2+sum(a[k]*cos(k*t*2*Pi/T)+b[k]*sin(k*t*2*Pi/T),k=1..n);

end:

Для сравнения графиков спектра амплитуд и спектра фаз зададим одинаковое число коэффициентов

> M:=Nk:

> a:=array(0..M):b:=array(0..M):c:=array(0..M);

> Setka(N,T,F_for_all,Y):

> F_Discret(Y,N,a,b,M):

> Cl:=blue,red,brown:;

> 15: RT:=seq(plot(Trig(t,8*k,T,a,b),t=-0.1..T+0.1,

numpoints=500,color=Cl[k]),k=1..3):

> 16: display(RT,GrafF);

> Spectr(M,a,b,c,'Grafphi');:

 

> display(ARR(M,c));:

> display(Ampl);

>

 

 

 

 

 

 

unit koren1;

interface

uses

  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

  Dialogs, StdCtrls, ExtCtrls;

type

  TForm1 = class(TForm)

    Label1: TLabel;

    Button2: TButton;

    Label2: TLabel;

    Button3: TButton;

    ListBox1: TListBox;

    Button4: TButton;

    Button5: TButton;

    Image1: TImage;

    Button1: TButton;

    procedure Button1Click(Sender: TObject);

    procedure Button2Click(Sender: TObject);

    procedure Button3Click(Sender: TObject);

    procedure Button4Click(Sender: TObject);

    procedure FormCreate(Sender: TObject);

    procedure Button5Click(Sender: TObject);

  private

    { Private declarations }

  public

    { Public declarations }

  end;

var

  Form1: TForm1;

implementation

{$R *.dfm}

{============================}

const

N=80;

m=N div 2;

x0=0;

T=3;

hx=T/N;

type

koeff=array[0..m] of real;

dann=array[0..N] of real;

var

Tau:real;

Y:dann;

a,b:koeff;

h:real;

eps:real;

La:real;

Rb:real;

Nkoeff:integer;

function pol(t:real):real;

begin

  pol:=sqr(sqr(t))-sqr(t)*t-3;

end;

function derive(t:real):real;

begin

  derive:=4*sqr(t)*t-3*sqr(t);

end;

function root(a,b:real):real;

var

X0,X1,delta:real;

begin

X0:=(a+b)/2;

repeat

   X1:=X0-pol(X0)/derive(X0);

   delta:=abs(X1-X0);

   X0:=X1;

until    delta<0.00001;

root:=X0;

end;

function riter(a,b:real):real;

var

X0,X1,delta:real;

const

lambda=0.0001;

begin

X0:=(a+b)/2;

repeat

   X1:=X0-pol(X0)*lambda;

   delta:=abs(X1-X0);

   X0:=X1;

until    delta<0.00001;

riter:=X0;

end;

function signal(t:real):real;

var

z:real;

begin 

   if t<0 then

   z:=0

   else

   if t<1/2 then

   z:=2*t

   else

   if t<1 then

   z:=1

   else z:=0;

   signal:=z;

end;

procedure Trig(m,N:integer;Y:dann;var a,b:koeff);

    var

    j,k:integer;

    p,q:real;

    x:real;

    h:real;

    begin

      h:=2*Pi/N;

      for k := 0 to m do

        begin

          p:=0;q:=0;

          for j := 1 to N do

            begin

              x:=j*h;

              p:=p+Y[j]*cos(x*k);

              q:=p+Y[j]*sin(x*k);

            end;

            a[k]:=p*2/N;

            b[k]:=q*2/N;

        end;

    end;

function Tpol(m:integer;x:real):real;

var

z:real;

k:integer;

begin

  z:=a[0]/2;

   for k:=1 to m do

    z:=z+(a[k]*cos(k*2*Pi/T*x)+b[k]*sin(k*2*Pi/T*x));

   Tpol:=z;

end;

 

 

procedure grafik(numvar:integer);

type

   dann= array[0..N] of real;

var

   L,R,W,H: integer;

   X: dann;

   Y: dann;

   k:integer;

   ymin,ymax:real;

   Mx,My:real;

   x0,y0: integer;

   posx,posy:integer;

   Nkf:string;

   tx:real;

   ypol:real;

procedure min_max(N:integer;Y:dann; var min, max:real);

   var

     k: integer;

   begin

      min:=Y[0];max:=Y[0];

      for k := 1 to N do

      if Y[k]> max then

         max:=Y[k]

        else if Y[k]< min then

         min:=Y[k];

      {увеличим диапазон}

       max:=max+0.1;

       min:=min-0.1;

   end;

begin

   L:=20;

   R:=form1.image1.clientHeight-20;

   W:=form1.image1.Width-50;

   H:=form1.image1.clientheight-50;

   case numvar of

1: begin

      for  k:=0 to N do

        X[k]:=signal(hx*k/Tau);

      min_max(N,X,ymin,ymax);

      Mx:=W/N;

      My:=H/(ymax-ymin);

      x0:=L;

      y0:=R-abs(Round(ymin*My));

      with form1.image1.Canvas do

      begin

        pen.Color:=clblue;

        font.Name:='Tahoma';

        font.Size:=8;

        font.Color:=claqua;

        for k:=0 to N do

         begin

           posx:=x0+round(k*Mx);

           posy:=y0-round(X[k]*My);

           textout(posx-2,posy-8,'o');

           Pixels[posx,posy]:=clRed;

         end;

 

 

        pen.Width:=2;

        Moveto(L,R);lineto(L,R-H);

        moveto(x0,y0);lineto(x0+W,y0);

        font.Color:=clred;

        textout(x0+W,y0+10,'x');

        textout(x0+W,y0-20,floattostrF(T,ffFixed,3,0));

        textout(x0+round(W*Tau/T), y0-20,'tau='+ floattostrF (Tau,ffFixed, 6, 3));

        Nkf:=Inputbox('Число коэффициентов полинома','например 10','20');

        Nkoeff:=strtoint(Nkf);

        pen.Color:=clNavy;

        tx:=0;

        ypol:=Tpol(Nkoeff,tx/Tau);

         posx:=x0+round(0*Mx/2);

         posy:=y0-round(ypol*My);

        moveto(posx,posy);

        for k:=1 to 2*N do

        begin

          tx:=hx*k/2;

          ypol:=Tpol(Nkoeff,tx/Tau);

          posx:=x0+round(k*Mx/2);

          posy:=y0-round(ypol*My);

          lineto(posx,posy);

        end;

      end;

    end;

   2: begin

        for  k:=0 to m do

         Y[k]:=sqrt(sqr(a[k])+sqr(b[k]));

         min_max(m,Y,ymin,ymax);

         Mx:=W/m;

         My:=H/(ymax-ymin);

         x0:=L;

         y0:=R-abs(Round(ymin*My));

      with form1.image1.Canvas do

      begin

         pen.Width:=2;

         pen.Color:=clred;

         Moveto(L,R);lineto(L,R-H);

         moveto(x0,y0);lineto(x0+W,y0);

         pen.Width:=5;

         pen.Color:=clblue;

        for  k:=0 to m do

         begin

          posx:=x0+round(k*Mx);

          posy:=y0-round(Y[k]*My);

          moveto(posx,y0);

          lineto(posx,posy);

         end;

      end;

      end;

   end;

end;

 

 

 

 

 

 

 

 

 

{===============================}

procedure TForm1.Button1Click(Sender: TObject);

begin

Form1.Close;

end;

procedure TForm1.Button2Click(Sender: TObject);

const

a=0; b=3;

begin

Tau:=riter(a,b);

Label2.Caption:='корень='+floattostr(Tau);