CPS 2015-03-31

dft2d.m

N=6;
M=4;
s = rand(N,M);
S = mat_dft(N) * s * mat_dft(M);
%backs = inv(mat_dft(N)) * S * inv(mat_dft(M));
backs = 1/N * conj(mat_dft(N)) * S * 1/M * conj(mat_dft(M))
err = max(abs(backs(:) - s(:)))

korohoda1.m

% m-plik skryptowy: DFT_przyklad_sygnalu.m
%
% Laboratorium:
% Cyfrowe Przetwarzanie Sygnalow
% Katedra Elektroniki, WIEiT, AGH
%
% Dane "realistyczne"
%
% Opracowanie: P.Korohoda, 13/03/2015;
% Modyfikacje: P.Korohoda, 23/03/2015;

clc; clear; close all;

N=2^10;
     n=0:N-1;  k=0:N-1; k=k(:);
     kn=k*n;
     
     w_N=exp(-j*2*pi/N);
     W=w_N.^kn;

     fp=1200; 
     Dt=1/fp;  df=fp/N;   t=0:Dt:(N-1)*Dt; 
     f=k*df;
     
        k0=300; % wybieramy indeks wskazujacy na numer wartosci "f";
        f1=k(end-k0)*df, % trafiamy w punkt na osi "f";
        %f1=((k(end-k0)+k(end-k0+1))/2)*df, % dokladnie w srodku miedzy punktami;
        
        dP=pi/4; % dla dP=0 otrzymujemy nieprzesuniety cosinus, dla dP=-pi/2 jest sinus;
     
     x=cos(2*pi*f1*t+dP); x=x(:);
     
     ta=0:Dt/100:t(end)/2;
     xa=cos(2*pi*f1*ta+dP);
     
     X=W*x;
     IW=inv(W);
     x1=real(IW*X); % nie dzielimy przez N, poniewaz wyznaczylismy macierz odwrotna do W;
     
     err=max(abs(x1-x)),

             k1=find(abs(X)>0.001);
             
        figure(1);
            subplot(1,2,1);
                plot(real(W(:)),imag(W(:)),'b.');    grid on; axis equal;
            subplot(1,2,2);
                plot(real(IW(:)),imag(IW(:)),'b.');  grid on; axis equal;
        figure(2);
                plot(ta,xa,'r-');  grid on; hold on;
                plot(t,x,'b.-');
        
        figure(3);       
            subplot(2,2,1); 
                plot(f,real(X),'b.');   grid on; hold on;
                plot([fp/2,fp/2],[min(real(X)),max(real(X))],'k','linewidth',2);
            subplot(2,2,3);
                plot(f,imag(X),'b.');   grid on; hold on;
                plot([fp/2,fp/2],[min(imag(X)),max(imag(X))],'k','linewidth',2);
            subplot(2,2,2);
                plot(f,abs(X),'b.');    grid on; hold on;
                plot([fp/2,fp/2],[min(abs(X)),max(abs(X))],'k','linewidth',2);
            subplot(2,2,4);
                plot(f,angle(X)/pi,'b.');  grid on; hold on;
                plot([fp/2,fp/2],[min(angle(X)/pi),max(angle(X)/pi)],'k','linewidth',2);
                plot(f(k1),angle(X(k1))/pi,'ro');

% KONIEC PLIKU;

korohoda3.m

% m-plik skryptowy: czas_obliczen_fft_vs_dft.m
%
% Laboratorium z Cyfrowego przetwarzania sygnalow;
%
% Autor: dr inż. P. Korohoda, Katedra Elektroniki, AGH;
% data utworzenia: 31 marca 2014 r.;

clc; clear; close all;

K=2^10;   
N0=2^12;
Nv=N0+(-99:0); 
m=0;
for N=Nv;
    m=m+1;
    x=rand(N,1);
    czas(m)=0;
    tic;
    for k=1:K,
           X=fft(x);
    end;
    czas(m)=toc;
end;
calkowity_czas_obliczen_w_minutach=sum(czas)/60, % w minutach;
czas=czas/K; % czas obliczen przeliczony na pojedyncze wyliczenie za pomoca FFT;

K1=16;
m=0;
for N=Nv;
    m=m+1;
    x=rand(N,1);
    A=fft(eye(N));
    czas1(m)=0;
    tic;
    for k=1:K1,
           X=A*x;
    end;
    czas1(m)=toc;
end;
calkowity_czas_obliczen_w_minutach1=sum(czas1)/60, % w minutach;
czas1=czas1/K1; % czas obliczen przeliczony na pojedyncze wyliczenie za pomoca FFT;


    figure(1);
        plot(Nv,czas/czas(end),'b.-');    grid on;    hold on;
        plot(Nv,czas1/czas(end)/2^10,'r.-');
        plot([Nv(1),Nv(end)],[1,1],'k');
            xlim([Nv(1),Nv(end)]);
                xlabel('N');  ylabel(['Ile razy dluzej niz dla N=',num2str(N0)]);
                title('Porownanie szybkosci obliczen FFT dla roznych dlugosci ciagu (N)');
        
%save czas_1.mat czas;

% KONIEC PLIKU;

korohoda4.m

% m-plik skryptowy: praca_z_plikiem_wav.m
%
% Laboratorium z Cyfrowego przetwarzania sygnalow;
%
% Autor: dr inż. P. Korohoda, Katedra Elektroniki, AGH;
% data utworzenia: 25 marca 2014 r.;

clc; clear; close all;

[x,fp]=wavread('utwor_2b.wav');
[N,kanaly]=size(x),
xL=x(:,1); N=length(xL); dt=1/fp; t=0:dt:(N-1)*dt;
XL=fft(xL); dfk=fp/N; fk=0:dfk:(N-1)*dfk;
        
            figure(1); 
                plot(ones(N,1),xL,'b.'); grid on;
            figure(2); 
                plot(t,xL,'b.-'); grid on;
            figure(3); n=20000+(1:5000);
                plot(t(n),xL(n),'b.-'); grid on;
            figure(4); k=1:N/2;
                subplot(2,1,1); plot(fk(k),abs(XL(k)),'b.-'); grid on;
                subplot(2,1,2); plot(fk(k),angle(XL(k))/pi,'b.-'); grid on;
                
% KONIEC PLIKU;

mat_dft,.m

function A=mat_dft(N)
% A=mat_dft(N);
%
% Laboratorium:
% Cyfrowe Przetwarzanie Sygnalow
% Katedra Elektroniki, WIEiT, AGH
%
% Funkcja wuznacza macierz DFT dla zadanej dlugosci ciagu: N;
% W matlabie mozna to wprawdzie zrobic proscie: A=fft(eye(N)) - ale tak nie cwiczymy
% odpowiedniej teorii;
%
% Opracowanie: P.Korohoda, 13/03/2015;

n=0:N-1; k=n; k=k(:);

A=exp(-j*2*pi/N).^(k*n);
end
% KONIEC FUNKCJI;