% Visualization of a polynomial on the complex domain:
% A circle of exponentially increasing radius is fed into the polynomial.
% We listen to and watch the polynomial's output.
% Of particular interest are the roots.
% Jörn Loviscach, 2009-12-26

duration = 60;
samplingRate = 48000;
framesPerSecond = 30;
numPoints = 10000;
width = 1280;

% be sure to set this to a meaningful folder
folderWhereToStore = 'C:\Users\jl\Desktop\Polynomial';

myRoots = [30+4i, 10-5i, 2+3i, 1-0.5i, -0.5+0.3i];
p = (0.5-1i)*poly(myRoots);
k = 1;
titleString(k) = {'Output: '};
for i = 1:length(p)-1
    if i ~= 1
        titleString(k) = strcat(titleString(k), '+');
    end
    titleString(k) = strcat(titleString(k), '(', num2str(p(i), '%3.1f'), ')z^{', num2str(length(p)-i), '}');
    if mod(i, 3) == 2
        k = k+1;
        titleString(k) = {''};
    end
end
titleString(k) = strcat(titleString(k), '+', num2str(p(length(p)), '%3.1f'));        

% compute the sound
t = 0:1/samplingRate:duration;
r = 10.^(t/duration*7-3);
frequency = 440/8.0;
z = r.*exp(1i*2*pi*frequency*t);
zDetuned = r.*exp(1i*2*pi*(frequency+0.865)*t);
zDetuned2 = r.*exp(1i*2*pi*(frequency+1.3761)*t);
y1 = polyval(p,z);
y1Detuned = polyval(p,zDetuned);
y1Detuned2 = polyval(p,zDetuned2);
d = abs(y1);
mi = zeros(1, length(y1));
mx = zeros(1, length(y1));
for i = 1:length(y1)
    imin = max(1, int32(i-samplingRate/frequency)); % one period before
    imax = min(length(d), int32(i+samplingRate/frequency)); % one period after
    mi(i) = min(d(imin:imax));
    mx(i) = max(d(imin:imax));
end
% Three waves with 120° phase shift each:
% among other benefits, any DC offset will cancel.
% When we are close to a root (mi/mx -> 0), the output is boosted.
% The division by mx ensures that the overall amplitude remains bounded.
y = 0.2*(y1+(1i)^(1/3)*y1Detuned-(1i)^(1/3)*y1Detuned2)./(50.0*mi+mx);
% wavplay([real(y)' imag(y)'], samplingRate, 'async');
wavwrite([real(y)' imag(y)'], samplingRate, strcat(folderWhereToStore,'\audio.wav'))

% compute the video
figure('Position', [100 100 width width*9/16])
videoPart = 0;
for frame = 0:framesPerSecond*duration
    % no correct avi files >= 2 GB, so split the video
    if mod(frame, framesPerSecond*15) == 0
        if frame ~= 0
            aviobj = close(aviobj);
        end
        videoPart = videoPart+1;
        aviobj = avifile(strcat(folderWhereToStore, '\video', num2str(videoPart), '.avi'), 'fps', framesPerSecond, 'compression', 'none');
    end    
    radius = 10^(frame/(framesPerSecond*duration)*7-3);
    phi = linspace(0, 2*pi, numPoints+1);
    z = radius*exp(1i*phi);
    y = polyval(p, z);
    mi = min(abs(y));
    mx = max(abs(y));
    
    subplot(121)
    hold off
    scatter(0, 0, 30, 'k')
    hold on
    interiorRoots = myRoots(abs(myRoots)<=radius);
    if ~isempty(interiorRoots)
        q1 = (abs(interiorRoots)/radius).^10; % normally is 0, but tends to 1 as we approach a root
        scatter(real(interiorRoots), imag(interiorRoots), 30+70*q1, [1 0 0], 'filled')
    end
    scatter(real(z), imag(z), 3, phi, 'filled')
    axis equal
    set(gca, 'GridLineStyle', '-')  
    grid on
    for k = 1:numPoints/20:numPoints
        text(real(z(k)), imag(z(k)), 42, [num2str((k-1)*360/numPoints), '°'], 'Color', [0.5 0.5 0.5])
    end
    title({'Input: z', ''}, 'FontWeight', 'bold', 'FontSize', 14)
    xlabel('Re', 'FontWeight', 'bold', 'FontSize', 14)
    ylabel('Im', 'FontWeight', 'bold', 'FontSize', 14)    
    
    subplot(122)
    hold off
    scatter(0.001, -0.001, 30, 'k') % 0.001: hack to prevent this marker from disappearing sporadically
    hold on
    q2 = (1-mi/mx)^10; % normally is 0, but tends to 1 as we approach a root
    scatter(0, 0, 100*q2, [q2 0 0], 'filled')
    scatter(real(y), imag(y), 3, phi, 'filled')
    axis equal
    set(gca, 'GridLineStyle', '-')   
    grid on
    for k = 1 : numPoints/20 : numPoints
        text(real(y(k)), imag(y(k)), 42, [num2str((k-1)*360/numPoints), '°'], 'Color', [0.5 0.5 0.5])
    end
    title(titleString, 'FontWeight', 'bold', 'FontSize', 14)
    xlabel('Re', 'FontWeight', 'bold', 'FontSize', 14)
    ylabel('Im', 'FontWeight', 'bold', 'FontSize', 14)
    
    aviobj = addframe(aviobj, getframe(gcf));
end
aviobj = close(aviobj);