function EN40_HW1_2016
%  EN40 2016 Homework 1  Solutions
%  17 POINTS TOTAL + 6 EXTRA CREDIT
close all  % This closes any open figure windows

%%
%   Problem 7:  Vector of equally spaced points starting at 0 and ending at 8192
%   You can also solve this problem with a loop, eg for i=1:2*8193 x(i) = 0.5*(i-1)/8192; end 
    x = 0:0.5:8192;
%   [1 POINT]
%%
%   Problem 8: Vectors y1-y5 creating sin waves with different frequencies    
    y1 = sin(4*pi*440*x/8192);
    y2 = sin(4*pi*550*x/8192);
    y3 = sin(4*pi*660*x/8192);
    y4 = sin(4*pi*880*x/8192);
    y5 = sin(4*pi*1760*x/8192);
%     sound(y1)
%     pause(2)
%     sound(y2)
%     pause(2)
%     sound(y3)
%     pause(2)
%     sound(y4)
%     pause(2)
    
%   [1 POINT]    
%%  
%   Problem 9 Plot of y1+y5 as a function of x (first 200 points only)
    figure1 = figure;
    % Create axes
    axes1 = axes('Parent',figure1,'YGrid','on','XGrid','on','FontSize',16);
    xlim(axes1,[0 6.2]);
    ylim(axes1,[-2 3]);
    box(axes1,'on');
    hold(axes1,'on');
    
    plot1 = plot(x(1:201)*0.5/8.192,y1(1:201)+y5(1:201),'LineWidth',2);
%    sound(y1+y5)
%    pause(2)

%   [2 POINTS, plus 1 extra credit for labels etc]    
    
%%  Problem 10: Filter the signal from problem 9 (see below for the   function   
    ynew = filter(y1+y5,0.2);
    sound(ynew)
    pause(2)
    plot2 = plot(x(1:201)*0.5/8.192,ynew(1:201),'LineWidth',2);
    ylabel('Signal (arbitrary units)');
    xlabel('Time (milliseconds)');
    title('Problems 10 and 11: Low pass filter demo');
    set(plot1(1),'DisplayName','Raw signal');
    set(plot2(1),'DisplayName','Filtered signal');
    legend1 = legend(axes1,'show');
    set(legend1,...
    'Position',[0.587500160695629 0.775803898024184 0.193856342540406 0.116949725164222],...
    'FontSize',14.4);    
    
%   [3 POINTS]    
    
%%  Problem 11 (optional, no credit) - this plays part of JSBach A minor sonata for unaccompanied flute     
%    [notes,accidental,octave,beats] = bach; 
%    for i = 1:length(notes)
%        singnote(notes(i),accidental(i),octave(i),beats(i),1,4)
%    end

%% Problem 12 Filtering a discontinuous function

    for i=1:length(x)
       yval = sin(880*pi*x(i)/8192);
       if (yval>0)
           r(i) = yval;
       else
           r(i) = 0;
       end
    end
    figure2 = figure;
    % Create axes
    axes1 = axes('Parent',figure2,'YGrid','on','XGrid','on','FontSize',16);
%    xlim(axes1,[0 12.4]);
    ylim(axes1,[-2 3]);
    box(axes1,'on');
    hold(axes1,'on');
    
    plot1 = plot(x(1:801)*0.5/8.192,r(1:801),'LineWidth',2);
    
    z = filter(r,0.01);
    plot2 = plot(x(1:801)*0.5/8.192,z(1:801),'LineWidth',2);
    ylabel('Signal (arbitrary units)');
    xlabel('Time (milliseconds)');
    title('Problems 10 and 11: Low pass filter demo');
    set(plot1(1),'DisplayName','Raw signal');
    set(plot2(1),'DisplayName','Filtered signal');
    legend1 = legend(axes1,'show');
    set(legend1,...
    'Position',[0.587500160695629 0.775803898024184 0.193856342540406 0.116949725164222],...
    'FontSize',14.4);    
    

%   [3 POINTS + 1 extra credit for pretty plots]
%%        
% Problem 13: Vidale-Wolfe advertising model
    
    K = 1;                      % Time constant, months^-1
    A = 10;                     % Spending on advertising, 1000$/mo
    Smax = 200;                  % Max revenue, 1000$/mo
    D = 0.3;                    % Time constant, months^-1
    initial_s = 0.;             % Initial sales rate
    
    [times1,sols1] = ode45(@(t,s) vidalewolfe(t,s,K,A,Smax,D),[0,36],initial_s);
    A = 20;
    [times2,sols2] = ode45(@(t,s) vidalewolfe(t,s,K,A,Smax,D),[0,36],initial_s);
    A = 60;
    [times3,sols3] = ode45(@(t,s) vidalewolfe(t,s,K,A,Smax,D),[0,36],initial_s);   
    
    figure3 = figure;
    axes1 = axes('Parent',figure3,'YGrid','on','XGrid','on','FontSize',16);
    box(axes1,'on');
    hold(axes1,'on');
    ylabel('Revenue (1000$/month)');
    xlabel('Time (months)');
    title('Problem 13: Predicted sales revenue rate');
    plot(times1,sols1,'DisplayName','A=10 000$/month','LineWidth',2);
    plot(times2,sols2,'DisplayName','A=20 000$/month','LineWidth',2);
    plot(times3,sols3,'DisplayName','A=60 000$/month','LineWidth',2);
    legend1 = legend(axes1,'show');
    set(legend1,...
    'Position',[0.435875590885719 0.650045242101745 0.439550771130598 0.224355152516907],...
    'FontSize',18);   

%    [3 POINTS + 2 extra credit for labeling plots, etc]

%%
%   Problem 14: Synthetic biology problem
    alpha = 4000;
    beta = 0.1;
    T = 5;
    
    initial_solution = [10,15,20,25,30,30];
    
    [times,sols] = ode45(@(t,w) genes(t,w,alpha,beta,T),[0,800],initial_solution);
    
    figure4 = figure;
    axes1 = axes('Parent',figure4,'YGrid','on','XGrid','on','FontSize',16);
    ylim(axes1,[0 4000]);
    hold(axes1,'on');
    ylabel('Protein concentration (#/bacterium)');
    xlabel('Time (seconds)');
    title('Problem 14: Synthetic Biology');
    plot1 = plot(times,sols(:,2:2:6),'LineWidth',2,'Parent',axes1);
    set(plot1(1),'DisplayName','p_1');
    set(plot1(2),'DisplayName','p_2');
    set(plot1(3),'DisplayName','p_3');
    legend1 = legend(axes1,'show');
    set(legend1,...
    'Position',[0.807044358436043 0.703592056965547 0.085799542194287 0.216666660641415],...
    'FontSize',14.4);

    figure4 = figure;
    axes1 = axes('Parent',figure4,'YGrid','on','XGrid','on','FontSize',16);
    ylim(axes1,[0 4500]);
    hold(axes1,'on');
    ylabel('Messenger concentration (#/bacterium)');
    xlabel('Time (seconds)');
    title('Problem 14: Synthetic Biology');
    plot1 = plot(times,sols(:,1:2:5),'LineWidth',2,'Parent',axes1);
    set(plot1(1),'DisplayName','m_1');
    set(plot1(2),'DisplayName','m_2');
    set(plot1(3),'DisplayName','m_3');
    legend1 = legend(axes1,'show');
    set(legend1,...
    'Position',[0.807044358436043 0.703592056965547 0.085799542194287 0.216666660641415],...
    'FontSize',14.4);    

 

%    [4 POINTS + 2 extra credit for labeling plots, etc]
end
    function z = filter(y,alpha)
%      Simple first order digital filter: y - vector to be filtered;
%      0<alpha<1  constant related to cutoff frequency of filter;
%      z - filtered signal.
        z = zeros(1,length(y));
        z(1) = y(1);
        for i = 2:length(y)
          z(i) = alpha*y(i) + (1.-alpha)*z(i-1);
        end
    end
    
    function [dsdt] = vidalewolfe(t,s,K,A,Smax,D)
%      Vidale Wolfe differential equation: t=time; s=revenue rate;
%      A,K,Smax,D are constants.
       dsdt = K*A*(1-s/Smax) - D*s;
    end
    
    function dwdt = genes(t,w,alpha,beta,T)
%      Differential equations for synthetic biology problem
%      w stores solution [m1,p1,m2,p2,m3,p3] at time t
%      dwdt is the vector of time derivatives
       m1 = w(1); p1 = w(2); m2 = w(3); p2 = w(4); m3 = w(5); p3 = w(6);
       dm1dt = (alpha/(1 + p3^2) - m1)/T;
       dp1dt = -beta*(p1-m1);
       dm2dt = (alpha/(1 + p1^2) - m2)/T;
       dp2dt = -beta*(p2-m2);
       dm3dt = (alpha/(1 + p2^2) - m3)/T;
       dp3dt = -beta*(p3-m3);
       dwdt = [dm1dt;dp1dt;dm2dt;dp2dt;dm3dt;dp3dt];      
    end
    
function singnote(tone,semitone,octave,beats,amplitude,beats_per_sec)
 %  singnote('A','b',4,1,0.5,2) sounds A4b quarter note at 2 quarter notes per sec 
 %  Frequencies use equal temper chromatic scale starting at C4
       scale = ['C','D','E','F','G','A','B','R'];  % 'R' is a rest
%      Naturals C       D      E     F      G     A     B 
       nats = [261.63,293.66,329.63,349.23,392.0,440,493.88,0,];
%      Sharps   C#       D#     E#     F#     G#    A#     B#
       sharps = [277.18,311.13,349.23,369.99,415.3,466.16,523.25];
%      Flats    Cb        Db      Eb    Fb      Gb    Ab     Bb
       flats =  [493.88/2,277.18,311.13,329.63,349.23,415.3,466.16];

       index = find(scale==tone);     
       freq = nats(index);
       if (semitone=='b') freq = flats(index); end
       if (semitone=='#') freq = sharps(index); end
       freq = freq * 2^(octave-4);
       x = 0:1:8192*beats/beats_per_sec;
       y = amplitude*sin(2*pi*freq*x/8192);
       sound(y);
       pause(0.9*beats/beats_per_sec)
        
        
end
function [notes,accidentals,octave,beats] = bach
% First movement of A minor sonata for unaccompanied flute, JS Bach
% Each 'notes' gives one measure.
  notes1 = ['R','E','A','G','A','C','A','E','A','E','A','G','A','C','A','E'];
  for i = 1:length(notes1)
      accidentals1(i) = 'n';
      beats1(i) = 1;
      octave1(i) = 5;
  end
  accidentals1(4) = '#';
  accidentals1(12) = '#';
  octave1(6) = 6;
  octave1(9) = 4;
  octave1(14) = 6;
  notes2 = ['A','C','E','F','G','F','E','D','C','E','G','A','E','D','C','B'];
  for i = 1:length(notes2)
      accidentals2(i) = 'n';
      beats2(i) = 1;
      octave2(i) = 5;
  end
  accidentals2(5) = '#';
  accidentals2(11) = '#';
  octave2(1) = 4;
  octave2(5) = 4;
  octave2(13) = 4;
  octave2(16) = 4;
  notes = [notes1,notes2,notes2];
  accidentals = [accidentals1,accidentals2,accidentals2];
  octave = [octave1,octave2,octave2];
  beats = [beats1,beats2,beats2];
  notes4 = ['A','C','G','A','F','D','G','A','E','C','G','A','D','C','B','A'];
  for i = 1:length(notes4)
      accidentals4(i) = 'n';
      beats4(i) = 1;
      octave4(i) = 4;
  end
  octave4(2) = 5;
  octave4(6) = 5;
  octave4(10) = 5;
  octave4(14) = 5;
  accidentals4(3) = '#';
  accidentals4(7) = '#';
  accidentals4(11) = '#';
  
  notes5 = ['G','E','G','B','E','G','B','D','F','D','B','G','E','D','C','B'];
   for i = 1:length(notes4)
      accidentals5(i) = 'n';
      beats5(i) = 1;
      octave5(i) = 4;
   end
  octave5(5:16) = 5;
  octave5(8) = 6;
  octave5(10) = 6;
  octave5(16) = 4;
  accidentals5(1) = '#';
  accidentals5(3) = '#';
  accidentals5(6) = '#';
  accidentals5(12) = '#';
  notes6 = ['C','A','C','E','A','C','E','F','G','B','G','E','C','B','A','G'];
   for i = 1:length(notes6)
      accidentals6(i) = 'n';
      beats6(i) = 1;
      octave6(i) = 5;
   end
   octave6(2) = 4;
   octave6(6:7) = 6;
   octave6(14:16) = 4;
   accidentals6(10) = 'b';
   accidentals6(14) = 'b';
   notes7 = ['A','F','A','C','D','F','A','C','B','G','B','D','E','G','B','D'];
   for i = 1:length(notes7)
      accidentals7(i) = 'n';
      beats7(i) = 1;
      octave7(i) = 5;
   end
   octave7(1:3) = 4;
   octave7(9:11) = 4; 
   notes8 = ['C','A','C','E','F','A','C','E','D','B','D','F','G','B','D','F'];
   accidentals8(1:16) = 'n';
   beats8(1:16) = 1;
   accidentals8(12) = '#';
   octave8(1:16) = 5;
   octave8(2) = 4;
   octave8(7) = 6;
   octave8(10) = 4;
   octave8(15) = 6;
   notes9 = ['E','G','C','B','C','E','C','G','C','G','C','B','C','E','C','G'];
   accidentals9(1:16) = 'n';
   beats9(1:16) = 1;
   octave9(1:16) = 5;
   octave9(3) = 6;
   octave9(5:7) = 6;
   octave9(11) = 6;
   octave9(13:15) = 6;
   notes10 = ['C','G','B','A','B','E','B','G','C','G','B','A','B','E','B','G'];
   accidentals10(1:16) = 'n';
   accidentals10(3) = 'b';
   accidentals10(5) = 'b';
   accidentals10(7) = 'b';
   accidentals10(11) = 'b';
   accidentals10(13) = 'b';
   accidentals10(15) = 'b';
   beats10(1:16) = 1;
   octave10(1:16) = 5;
   octave10(6) = 6;
   octave10(14) = 6;
   notes11 = ['C','G','B','A','B','E','B','G','C','G','E','C','A','G','F','E'];
   accidentals11(1:16) = 'n';
   accidentals11(1) = '#';
   accidentals11(3:2:7) = 'b';
   accidentals11(9) = '#';
   accidentals11(12) = '#';
   beats11(1:16) = 1;
   octave11(1:16) = 5;
   octave11(6) = 6;
   octave11(13:16) = 4;
   notes12 = ['D','G','F','E','F','A','F','D','B','G','B','D','E','D','C','B'];
   accidentals12(1:16) = 'n';
   accidentals12(9) = '#';
   beats12(1:16) = 1;
   octave12(1:16) = 5;
   octave12(1) = 4;
   octave12(9:11) = 4;
   octave12(16) = 4;
   notes13 = ['C','E','A','G','A','C','A','F','D','B','D','F','B','A','G','F'];
   accidentals13(1:16) = 'n';
   accidentals13(4) = '#';
   accidentals13(8) = '#';
   accidentals13(9) = '#';
   accidentals13(11:12) = '#';
   accidentals13(16) = '#';
   beats13(1:16) = 1;
   octave13(1:16) = 5;
   octave13(6) = 6;
   octave13(10) = 4;
   notes14 = ['G','F','E','G','A','E','G','C','F','E','D','F','G','D','F','B'];
   accidentals14(1:16) = 'n';
   accidentals14(2) = '#';
   accidentals14(9) = '#';
   accidentals14(12) = '#';
   accidentals14(15) = '#';
   beats14(1:16) = 1;
   octave14(1:16) = 5;
   octave14(5) = 4;
   octave14(8) = 6;
   octave14(13) = 4;
   notes15 = ['E','D','C','E','F','E','F','A','D','B','C','A','D','C','B','A'];
   accidentals15(1:16) = 'n';
   accidentals15(5) = '#';
   accidentals15(7) = '#';
   accidentals15(9) = '#';
   accidentals15(13) = '#';
   beats15(1:16) = 1;
   octave15(1:16) = 5;
   octave15(5) = 4;
   octave15(10) = 4;
   octave15(12) = 4;
   octave15(13) = 4;
   octave15(15:16) = 4;
   notes16 = ['G','E','F','D','G','F','E','D','A','F','G','E','A','G','F','E'];
   accidentals16(1:16) = 'n';
   accidentals16(5) = '#';
   accidentals16(10) = '#';
   accidentals16(13) = '#';
   accidentals16(15) = '#';
   beats16(1:16) = 1;
   octave16(1:16) = 5;
   octave16(1) = 4;
   octave16(5) = 4;
   octave16(9) = 4;
   octave16(13) = 4;
   notes17 = ['D','B','G','D','C','A','F','C','B','G','E','B','A','F','D','A'];
   accidentals17(1:16) = 'n';
   accidentals17(1) = '#';
   accidentals17(3) = '#';
   accidentals17(5) = '#';
   accidentals17(7) = '#';
   accidentals17(12) = 'b';
   accidentals17(14:15) = '#';
   beats17(1:16) = 1;
   octave17(1:16) = 5;
   octave17(9) = 4;
   octave17(12:13) = 4;
   octave17(16) = 4;
   notes18 = ['G','E','C','A','F','C','A','F','D','F','A','C','B','A','G','F'];
   accidentals18(1:16) = 'n';
   accidentals18(5) = '#';
   accidentals18(8:10) = '#';
   accidentals18(16) = '#';
   beats18(1:16) = 1;
   octave18(1:16) = 5;
   octave18(1) = 4;
   octave18(4:5) = 4;
   octave18(7:11) = 4;
   octave18(13) = 4;
   notes19 = ['G','F','E','B','C','F','E','E'];
   accidentals19(1:8) = 'n';
   accidentals19(2) = '#';
   accidentals19(6) = '#';
   beats19(1:8) = 1;
   beats19(2:3) = 0.5;
   beats19(6) = 3;
   beats19(8) = 8;
   octave19(1:8) = 5;
   octave19(5) = 6;
   notes = [notes1,notes2,notes2,notes4,notes5,notes6,notes7,notes8,notes9,notes10,notes11,notes12,notes13,...
            notes14,notes15,notes16,notes17,notes18,notes19];
  accidentals = [accidentals1,accidentals2,accidentals2,accidentals4,accidentals5,accidentals6,...
                 accidentals7,accidentals8,accidentals9,accidentals10,accidentals11,accidentals12,accidentals13,...
                 accidentals14,accidentals15,accidentals16,accidentals17,accidentals18,accidentals19];
  beats = [beats1,beats2,beats2,beats4,beats5,beats6,beats7,beats8,beats9,beats10,beats11,beats12,beats13,...
           beats14,beats15,beats16,beats17,beats18,beats19];
  octave = [octave1,octave2,octave2,octave4,octave5,octave6,octave7,octave8,octave9,octave10,octave11,octave12,...
           octave13,octave14,octave15,octave16,octave17,octave18,octave19];
end

