%  Abgabe: 04.01.2012


%% ANMERKUNG:
%%           --> geplottetter Fehler gibt fälschlicherweise die Ausgangsfunktion wieder
%%           --> Aufgabe 4 wurde nicht vollständig programmiert 
%%           --> Die Umsetzung von 3 Perioden im vorgegebenen Intervall ist nicht für alle Aufgaben umgesetzt worden
%%           
%%              

%% Aufgabenstellung
%Leiten Sie die Bestimmungsgleichungen der Fourierkoeffizienten 
%einer Funktion mit beliebiger Periode T her, und erläutern Sie 
%dabei die Begriffe Basisfunktion, Orthogonalbasis und Frequenzplot. 
%Klären Sie die Fragen nach der Existenz und der Konvergenz. 
%Erklären Sie das Gibbssche Phänomen. 
%Berechnen Sie die Fourierreihenentwicklung für:
% (1) f(x)=x^2; -1<x<1
% (2) f(x)=x^2; 0<x<2
% (3) f(x)=abs(x)-1; -1<x<0
%     f(x)=1-(x-1)^2; 0<x<2
% (4) diskretes Signal (Datei pv.mat)
%Gibt es MATLAB-Funktionen für die Fourierreihenentwicklung? 
%Wenn ja, welche und wie kann man sie für ihre Probleme nutzen?

%% Programmlisting

clear all; clc;       % WorkspaceVariablen)löschen. Löscht alle Ein- und Ausgaben im Comand-Window
    syms x AK BK      % mit der Anweisung syms wird x als symbolische Variable definiert
    
%% Benutzereingabe:  
    Eingabe= input('Gewünschte Aufgabe eingeben und mit Eingabetaste bestätigen, Nr. 1,2 oder 3:');
    
    

    
%% Aufgabenstellung 1 bis 3 als switch-Befehl: 
    Auswahl =1;  %Variable zum plotten von stetigen und unstetigen Funktionen (1 für Aufgabe 1 und 2, 2 für Aufgabe 3) 
    
    switch Eingabe
   
        
        case 1
                fx=x^2;          % Funktion der Aufgabe 1 
                a=-1; d=1;
           
                %p=-1+d-a/3
                %q=p+d-a/3
                
                
                       % Intervall der Aufgabe 1
                xmin=-1;xmax=1;ymin=0.9;ymax=1;% Definierte Achsenwerte für die grafische Darstellung
        case 2
                fx=x^2+x^2+x^2;        % Funktion der Aufgabe 2
                a=0;d=2;       % Intervall der Aufgabe 2
                xmin=0;xmax=2;ymin=1.5;ymax=4; % Definierte Achsenwerte für die grafische Darstellung
                
        case 3
                Auswahl =2;
                fxa=abs(x)-1;  % 1.Funktion der Aufgabe 3
                a=-1; b=0;     % 1.Intervall der Aufgabe 3
                
                fxb=1-(x-1)^2; % 2.Funktion der Aufgabe 3
                c=0;d=2;       % 2.Intervall der Aufgabe 3
                xmin=-1;xmax=2;ymin=-1;ymax=0.5; % Definierte Achsenwerte für die grafische Darstellung
    end


%% Definition der Periode:
    T=d-a/3;                                      %  Periode der Funktion 
    h=T/2;                                        %  halbe Periode   
    
    
%% Differenzierung des Intervalls:        
 if Auswahl == 1                          % Bei den Aufgaben 1 und 2         
  b=a+h;                                   %  Mittelpunkt des Intervalls
  c=b; 
  fxa=fx;                                  %  fxa=fxb=fx für Aufgabe 1 und 2
  fxb=fx; 
   
 else if Auswahl == 2                     % Bei Aufgabe 3, als zusammengesetzte 
                                          % Funktion,fx= fxa und fxb
      
     else
     end
 end
 
%% Berechnung der Fourierkoeffizienten:
  ak=int(fxa,x,a,b)/h+int(fxb,x,b,d)/h;                % Berechnung von a0
        bk=[];                                               % b0 ist leer 
        f=ak/2;                                               %a0/2                                                        
  
        farb= ['k','y','g','b','r','m','g','k','b'];  %  Ein Feld mit Farbcodes für Kurvenfarben
        F=0;                                          %  Farb-Variable für die dargestellten Funktion mit dem Anfangswert Null 
       
       % Berechung der Fourierkoeffizienten ak,bk, und der Funktion f  
        for i=1:10 % Anzahl der Rechnschritte (i=start:ende)
        aki=int(fxa*cos(i*pi*x/h),x,a,b)/h+int(fxb*cos(i*pi*x/h),x,c,d)/h; %Berechnung der ak-Koeffizienten
        bki=int(fxa*sin(i*pi*x/h),x,a,b)/h+int(fxb*sin(i*pi*x/h),x,c,d)/h; %Berechnung der bk-Koeffizienten
        ak=[ak,aki];  % die berechneten Werte von aki werden fortlaufend ak zugewiesen
        bk=[bk,bki];  % die berechneten Werte von bki werden fortlaufend bk zugewiesen
%% Fourier-Reihe f als Summenfunktion:        
        f=f+aki*cos(i*pi*x/h)+bki*sin(i*pi*x/h); % Berechnung der Näherungsfunktion f mit den Fourierkoeffizienten
%% Visualisierung der Fourier-Reihe:      
        F=F+1;                        % Farbvariable wird um 1 erhöht und der 1. Funktionsverlauf wird mit Farbcode 'k' geplottet
        if(i<=8) ||(i==9)                % Anweisung zum plotten der ersten 10 Kurven (inkl. fx >>Zeile 101)
%% figure 1 stellt die ersten 10 kumulierten Werte der Fourierreihenentwicklung dar:       
        figure (1)                    % 1.Grafikfenster wird geöffnet zur Darstellung der Funktion
        [xx,yy]=fplot(char(f),[a,d])  % Plotten der Funktion f im Intervall a bis b
        plot(xx,yy,farb(F));
        hold on;                      % Bei Neuplot wird alte Grafik beibehalten
        grid on;                      % Gitternetz ein
%% figure 2 figure 4 zoomt auf vordefinierte Achenwerte:    
       %figure (2)                    % 2.Grafikfenster wird geöffnet zur Darstellung der Funktion
       %[xx,yy]=fplot(char(f),[a,d])  % Plotten der Funktion f im Intervall a bis b
       %plot(xx,yy,farb(F));
       %axis([xmin,xmax,ymin,ymax]);  % Grafische Darstellung wird im gewünschten Achsenbereich dargestellt
       %hold on;                      % Bei Neuplot wird alte Grafik beibehalten
       %grid on;                      % Gitternetz ein 
   
% Modulator, damit nur die 10., 20., 30., 40. und 50. Fourier-Reihe dargestellt wird:              
       else if(mod(i,10)==0)         % nur um Zahlen ohne Rest 
%% figure 3 stellt bis 50 die kumulierten Werte der Fourierreihenentwicklung dar:
              figure (3)             % 3.Grafikfenster wird geöffnet zur Darstellung der Funktion
              F=i/10;                % zur Verknüpfung der Anzahl der Schritte mit der Farbscala
              [xx,yy]=fplot(char(f),[a,d]); % plotten der 
              plot(xx,yy,farb(F));  % hier wird jeder 10. Plot dargestellt (10,20,30 ....)
              hold on               % Bei Neuplot wird alte Grafik beibehalten 
              grid on               % Gitternetz ein
              
%% figure 4 figure 4 zoomt auf vordefinierte Achenwerte:              
              %figure (4)            % 4.Grafikfenster wird geöffnet zur Darstellung der Funktion
              %F=i/10;
              %[xx,yy]=fplot(char(f),[a,d]); 
              %plot(xx,yy,farb(F));
              %axis ([xmin,xmax,ymin,ymax]); % Grafische Darstellung wird im gewünschten Achsenbereich dargestellt
              %hold on   % Bei Neuplot wird alte Grafik beibehalten
              %grid on   %Gitternetz ein
           end  
    end
        end
        
 %% Visualisierung der gewählten Ausgangsfunktion (1,2 oder 3) (s. Aufgabenstellung):          
   if Auswahl == 1 %Nachfolgender Ausdruck soll nur für die Aufgaben 1 und 2 gelten
   [xx,yy]=fplot(char(fx),[a,d]); % Ausgangsfunktion soll über vorgegebenes Intervall gezeichnet werden

   % Ausgangsfunktion wird in die figure 1 reingezeichnet
    figure(1) 
     plot(xx,yy,'k','linewidth',2);
     xlabel('x');ylabel('f(x)');                                                  
     legend('k=1','k=2','k=3','k=4','k=5','k=6','k=7','k=8','k=9','f(x)','Location',['North']); % Einfügen einer Legende.
     title ('[fx]')
     hold off   % Bei Neuplot wird alte Grafik löschen
   
    % Ausgangsfunktion wird in die figure 2 reingezeichnet
    %figure(2)
     %plot(xx,yy,'k','linewidth',2);
     %axis ([xmin,xmax,ymin,ymax])
     %xlabel('x');ylabel('f(x)');                                                  
     %title ('[fx]')
     %legend('k=1','k=2','k=3','k=4','k=5','k=6','k=7','k=8','k=9','f(x)','Location',['North']); % Einfügen einer Legende.
     %hold off   
     
   % Ausgangsfunktion wird in die figure 3 reingezeichnet
     figure(3)
     plot(xx,yy,'k','linewidth',2);
     xlabel('x');ylabel('f(x)');
     legend('k=10','k=20','k=30','k=40','k=50','f(x)','Location',['North']); % Einfügen einer Legende.
     hold off  % Gitternetz aus  
    
     %% figure 4 zoomt auf vordefinierte Achenwerte
     %figure(4)
     %plot(xx,yy,'k','linewidth',2);
     %axis ([xmin,xmax,ymin,ymax])
     %xlabel('x');ylabel('f(x)');
     %legend('k=10','k=20','k=30','k=40','k=50','f(x)','Location',['North']); %Einfügen einer Legende.
     %hold off  % Gitternetz aus

     
%% figure 7 visualisiert den Fehler:     
     figure(7)
     
     ff=f-fx;
     [xf,yf]=fplot(char(ff),[a,d]);
     plot(xf,yf,'r','linewidth',2);
     plot(xx,yy,'b','linewidth',2);
     %surf(xx,yy, f-fx) %Geht nicht
     xlabel('x');ylabel('f(x)');                                                  
     legend('ff','f(x)','Location',['North']);        % Einfügen einer Legende.
    hold off;
     
     
   else if Auswahl == 2
      
     [xx1,yy1]=fplot(char(fxa),[a,b]);
     [xx2,yy2]=fplot(char(fxb),[c,d]);
     
     figure(1)                     

     plot(xx1,yy1,'c','linewidth',2);
     plot(xx2,yy2,'c','linewidth',2);
     axis ([xmin,xmax,ymin,ymax])
     xlabel('x');ylabel('f(x)');                                                  
     legend('k=1','k=2','k=3','k=4','k=5','k=9','Location',['North']);        
     hold off                                                                    
     grid on                                                                             
   
     figure(2)

     plot(xx1,yy1,'c','linewidth',2);
     plot(xx2,yy2,'c','linewidth',2);
     axis ([xmin,xmax,ymin,ymax])
     xlabel('x');ylabel('f(x)');
     legend('k=10','k=20','k=30','k=40','k=50','Location',['North']);
     hold off
     grid on
      
      figure(3)
     plot(xx1,yy1,'c','linewidth',2);
     plot(xx2,yy2,'c','linewidth',2);
     xlabel('x');ylabel('f(x)');
     legend('k=10','k=20','k=30','k=40','k=50','f(x)','Location',['North']);
     hold off   
    
     figure(4)
     plot(xx1,yy1,'c','linewidth',2);
     plot(xx2,yy2,'c','linewidth',2);
     axis ([xmin,xmax,ymin,ymax])
     xlabel('x');ylabel('f(x)');
     legend('k=10','k=20','k=30','k=40','k=50','f(x)','Location',['North']);
     hold off    
     
%% figure 5 visualisiert den Frequenzplot der Fourierkoeffizienten Ak:     
     figure(5)
     
     if a<=x<b;
     ff1=fxa-f;
     else if c<=x<=d;
         
         ff2=fxb-f;
         end
         end
     [xf1,yf1]=fplot(char(ff1),[a,b]);
     [xf2,yf2]=fplot(char(ff2),[c,d]);
     plot(xf,yf,'r','linewidth',2);
     plot(xx1,yy1,'c','linewidth',2);
     plot(xx2,yy2,'c','linewidth',2);
     xlabel('x');ylabel('f(x)');                                                  
     legend('ff','f(x)','Location',['North']);        % Einfügen einer Legende.
     title ('[fx]')
     
       end
   end
   
   for i=1:2
      
        if i==1
         
         Frequenz_Plot=ak;
         Frequenz_Plot_Name=char(AK);
          
        else
            
         Frequenz_Plot=bk;        
         Frequenz_Plot_Name=char(BK);
       
        end


        figure(4+i)
     
     subplot(2,1,1)
    
     if Frequenz_Plot_Name==char(AK);
        
         stem(0:1:5,double(Frequenz_Plot(1:6)),'b.');    
     else
         stem(1:1:5,double(Frequenz_Plot(1:5)),'b.'); 
         
     end
     
     if Auswahl == 1                            
       title(strcat('Frequenzplot')); 
     end
     
     xlabel('k');ylabel(strcat(Frequenz_Plot_Name)); 
     grid on
     set(gca,'xgrid','off');                          
   
     subplot(2,1,2)
    
      if Frequenz_Plot_Name==char(AK)
        
         stem(5:1:(length(Frequenz_Plot)-1),double(Frequenz_Plot(6:length(Frequenz_Plot))),'b.'); 
    
    else
    
        stem(5:1:length(Frequenz_Plot),double(Frequenz_Plot(5:length(Frequenz_Plot))),'b.');   
      end
    
    if Auswahl == 1                            
        title(strcat(' Frequenzplot')); 
    end
    
    xlabel('k');ylabel(strcat(Frequenz_Plot_Name)); 
    grid on
    set(gca,'xgrid','off');
    
    end