function zweus(N, fmin, fmax)
%
%   zweus(N, fmin fmax) original stimulus and its spectrum after a compressive
%   nonlinearity according to Y(t) = X(t)^0.3
%
if nargin<3 fmax = 3000; end
if nargin<2 fmin = 2000; end
if nargin<1 N = 7; end
   
    dt = 1/(2*fmax);
    t=0:dt:1.0-dt;  D = length(t);
    Stim = zeros(1, D); Env = zeros(1,D);
    
    dF = fmax - fmin;
    prims = fmin + floor(dF*rand(1,N));
    Prims = sort(prims, 'ascend')
     
     %Amps = 0.1 + rand(1,N); 
     Amps = 5 + 5*rand(1,N); 
     Phis = 2*pi*(rand(1,N)-0.5);
 
     for n=1:N
         Stim = Stim + Amps(n)*cos(2*pi*Prims(n)*t + Phis(n));
     end
     
     figure(1);
     clf
     T=floor(D);
     plot(t(1:T), Stim(1:T), 'k-','linewidth', 1.5);
     axis([0,1,-1.1*max(abs(Stim)), 1.1*max(abs(Stim))]);
     set(gca, 'fontsize',13,'fontweight','bold');
     xlabel('Time (s)', 'fontsize',13,'fontweight','bold');
     ylabel('Amplitude (a.u.)', 'fontsize',13,'fontweight','bold');
     text(0.01, max(abs(Stim)), 'A', 'fontsize', 24, 'fontweight', 'bold');
     
     M = faculty(N)/(faculty(N-2)*2);    % number of components
     Amn = zeros(1,M);
     Pmn = zeros(1,M);
     Fmn = zeros(1,M);
     k=0;
     for n=1:N
         for m=n+1:N
           k=k+1;
           Fmn(k) = Prims(m)-Prims(n);
           Amn(k) = 2*Amps(n)*Amps(m);
           Pmn(k) = Phis(m)-Phis(n);
         end
     end
     
     AMP =0;
     for k=1:N
         AMP = AMP + Amps(k)^2;
     end
     
     figure(2); clf;
     plotspec(Prims, Amps, Phis);
     text(min(Prims)-80, 0.95, 'B', 'fontsize', 24, 'fontweight', 'bold');
     
     figure(3); clf;
     plotspec(Fmn, Amn, Pmn);
     text(min(Fmn)-80, 0.95, 'C', 'fontsize', 24, 'fontweight', 'bold');
     
     figure(4); clf
     for n=1:M
         Env = Env + Amn(n)*cos(2*pi*Fmn(n)*t + Pmn(n));
     end
     Env = AMP + Env;
     T=floor(D/10);
     plot(t(1:T), Env(1:T), 'r-','linewidth', 2);
     hold on
     plot(t(1:T), Stim(1:T).^2, 'k-','linewidth', 1);
     axis([0 t(T) -0.1 1.1*max(Env(1:T))]);
     set(gca, 'fontsize',13,'fontweight','bold');
     xlabel('Time (s)', 'fontsize',13,'fontweight','bold');
     ylabel('Amplitude (a.u.)', 'fontsize',13,'fontweight','bold');
     text(0.01*t(T), 1.05*max(Env(1:T)), 'D', 'fontsize', 24, 'fontweight', 'bold'); 
     
     %
     %  part 2: compressive nonlinearity N(t) = Env(t)^0.3 and the spectrum of its envelope
     %
     NonLin = Env.^0.3;
     
     figure(5); clf;
     amp = max(Env(1:T));
     plot(t(1:T), Env(1:T)/amp, 'k-','linewidth', 1);
     hold on
     amp2 = max(NonLin(1:T).^2);
     plot(t(1:T), (NonLin(1:T).^2)/amp2, 'r-','linewidth', 2);
     axis([0 t(T) -0.1 1.1]);
     set(gca, 'fontsize',13,'fontweight','bold');
     xlabel('Time (s)', 'fontsize',13,'fontweight','bold');
     ylabel('Amplitude (a.u.)', 'fontsize',13,'fontweight','bold');
     text(0.05*T, 1.05*max(Env(1:T)), 'Compressive Nonlinearity: x^0.3','fontsize', 12, 'fontweight', 'bold'); 
      
     %
     %  spectrum of the envelope and of it's nonlinear deformed version
     %
     figure(6); clf; 
     powerspectrum(Env, 2^14, 1/dt, [20 1000 20 70], 'k-', 0);
     hold on
     powerspectrum(NonLin, 2^14, 1/dt, [20 1000 20 70], 'r-', 37);
     plot([25, 35], [68, 68], 'k-');
     text(37, 68, 'Envelope spectrum', 'fontsize', 11, 'fontweight', 'bold');
     plot([25, 35], [66, 66], 'r-');
     text(37, 66, 'Spectrum (Env)^{0.3}', 'fontsize', 11, 'fontweight', 'bold');
     return
     
    function plotspec(F, A, P)
        N=length(F);
        mx = max(F); mn = min(F);
        ax = max(A); am = min(A);
        px = max(abs(P)); pm = min(abs(P));
        
        plot([mn, mx], [0.55,0.55],'k-');
        hold on
        axis([mn - 100 , mx + 100 , 0, 1]);
        
        for n=1:N
            plot( [F(n), F(n)], [0.55, 0.55+0.35*A(n)/ax],'k-','linewidth',2);
        end
        
        plot( [mx+30, mx+30], [0.55,0.9],'k-','linewidth',3);
        plot( [mx+28, mx+32], [0.55,0.55],'k-','linewidth',3);
        plot( [mx+28, mx+32], [0.9,0.9],'k-','linewidth',3);
        text(mx+35, 0.54, '0', 'fontsize',11,'fontweight','bold');
        a=sprintf('%3.1f',ax);
        text(mx+35, 0.90, a, 'fontsize',11,'fontweight','bold');
        
        plot( [mx+30, mx+30], [0.10,0.4],'k-','linewidth',3);
        plot( [mx+28, mx+32], [0.4,0.4],'k-','linewidth',3);
        plot( [mx+28, mx+32], [0.25,0.25],'k-','linewidth',3);
        plot( [mx+28, mx+32], [0.10,0.10],'k-','linewidth',3);
        text(mx+35, 0.09, '-pi', 'fontsize',11,'fontweight','bold');
        text(mx+35, 0.40, '+pi', 'fontsize',11,'fontweight','bold');
        text(mx+35, 0.24, '  0', 'fontsize',11,'fontweight','bold');
        
        plot([mn, mx], [0.25,0.25],'k-');
        for n=1:N
            plot( [F(n), F(n)], [0.25, 0.25+0.15*P(n)/pi],'k-','linewidth',2);
        end
        
        set(gca, 'fontsize',12,'fontweight','bold');
        set(gca, 'ytick',[]);
        xlabel('Frequency (Hz)', 'fontsize',13,'fontweight','bold');
        ylabel('Phase (rad)              Amplitude (a.u.)', 'fontsize',13,'fontweight','bold');
        
        return
        
function [P,F]=powerspectrum(X,N,Fs,scales, col, add)   % adapted from Jeroen Goossens

Noverlap = 0;
Window   = hanning(N);
Pr       = 0.95;

[P,F] = spectrum(X,N,Noverlap,Window,Fs,Pr);

% decibel scale
Spt=10*log10(P(:,1))+add;

semilogx(F,Spt,col);
xlabel('Frequency (Hz)','fontsize',13,'fontweight','bold');
ylabel('Signal Power (dB)','fontsize',13,'fontweight','bold');
title('power spectral density');
set(gca, 'fontsize',12,'fontweight','bold');
axis(scales);
        