### Archive

Archive for the ‘youtube demo code’ Category

## plotting the frequency content of a signal

A ‘web friendly’ version of this code with figures is available at https://dadorran.wordpress.com/2014/02/20/plotting-frequency-spectrum-using-matlab/

%% Demonstration of how to plot the Frequency Spectrum of a Signal
%https://www.dropbox.com/s/n4kpd4pp2u3v9nf/tremor_analysis.txt
N = length(signal);
fs = 62.5; % 62.5 samples per second
fnyquist = fs/2; %Nyquist frequency

%% Quick view of double-sided (two-sided) magnitude spectrum
% When roughly interpreting this data half way along x-axis
% corresponds to half the sampling frequency
plot(abs(fft(signal)))
xlabel('Frequency (Bins - almost!)')
ylabel('Magnitude');
title('Double-sided Magnitude spectrum');
axis tight

%% Double-sided magnitude spectrum with frequency axis (in bins)
fax_bins = [0 : N-1]; %N is the number of samples in the signal
plot(fax_bins, abs(fft(signal)))
xlabel('Frequency (Bins)')
ylabel('Magnitude');
title('Double-sided Magnitude spectrum (bins)');
axis tight

%% Single-sided magnitude spectrum with frequency axis in bins
X_mags = abs(fft(signal));
fax_bins = [0 : N-1]; %frequency axis in bins
N_2 = ceil(N/2);
plot(fax_bins(1:N_2), X_mags(1:N_2))
xlabel('Frequency (Bins)')
ylabel('Magnitude');
title('Single-sided Magnitude spectrum (bins)');
axis tight

%% Single-sided magnitude spectrum with frequency axis in Hertz
% Each bin frequency is separated by fs/N Hertz.
X_mags = abs(fft(signal));
bin_vals = [0 : N-1];
fax_Hz = bin_vals*fs/N;
N_2 = ceil(N/2);
plot(fax_Hz(1:N_2), X_mags(1:N_2))
xlabel('Frequency (Hz)')
ylabel('Magnitude');
title('Single-sided Magnitude spectrum (Hertz)');
axis tight

%% Single-sided magnitude spectrum with frequency axis normalised
% Normalised to Nyquist frequency.
% Very common to use this method of normalisation in matlab
X_mags = abs(fft(signal));
bin_vals = [0 : N-1];
fax_norm = (bin_vals*fs/N)/fnyquist; % same as bin_vals/(N/2)
N_2 = ceil(N/2);
plot(fax_norm(1:N_2), X_mags(1:N_2))
xlabel({'Frequency (Normalised to Nyquist Frequency. ' ...
'1=Nyquist frequency)'})
ylabel('Magnitude');
title('Single-sided Magnitude spectrum (Normalised to Nyquist)');
axis tight

%% Single-sided magnitude spectrum - frequency in rads per sample
X_mags = abs(fft(signal));
bin_vals = [0 : N-1];
N_2 = ceil(N/2);
ylabel('Magnitude');
axis tight

%% Double-sided magnitude spectrum showing negative frequencies
% See http://youtu.be/M1bLPZdNCRA for an explanation of negative
% frequencies
X_mags = abs(fftshift(fft(signal)));
bin_vals = [0 : N-1];
N_2 = ceil(N/2);
fax_Hz = (bin_vals-N_2)*fs/N;
plot(fax_Hz, X_mags)
xlabel('Frequency (Hz)')
ylabel('Magnitude');
title('Double-sided Magnitude spectrum (Hertz)');
axis tight

%% Single-sided magnitiude spectrum in decibels and Hertz
X_mags = abs(fft(signal));
bin_vals = [0 : N-1];
fax_Hz = bin_vals*fs/N;
N_2 = ceil(N/2);
plot(fax_Hz(1:N_2), 10*log10(X_mags(1:N_2)))
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)');
title('Single-sided Magnitude spectrum (Hertz)');
axis tight

%% Single-sided power spectrum in decibels and Hertz
X_mags = abs(fft(signal));
bin_vals = [0 : N-1];
fax_Hz = bin_vals*fs/N;
N_2 = ceil(N/2);
plot(fax_Hz(1:N_2), 20*log10(X_mags(1:N_2)))
xlabel('Frequency (Hz)')
ylabel('Power (dB)');
title('Single-sided Power spectrum (Hertz)');
axis tight

%% Single-sided power spectrum in dB and frequency on a log scale
X_mags = abs(fft(signal));
bin_vals = [0 : N-1];
fax_Hz = bin_vals*fs/N;
N_2 = ceil(N/2);
semilogx(fax_Hz(1:N_2), 20*log10(X_mags(1:N_2)))
xlabel('Frequency (Hz)')
ylabel('Power (dB)');
title({'Single-sided Power spectrum' ...
' (Frequency in shown on a log scale)'});
axis tight



## Basic Matlab GUI

Here is the code I used to create a basic matlab GUI without using GUIDE in the following youtube videos:
First Video: GUI Layout http://youtu.be/kDcn83becS0
Second Video: Object handles and Object Manipulation http://youtu.be/XGRkqgDwWRA
Third Video: Callback Functions http://youtu.be/JT8gAbY2Ooc
Fourth Video: Finishing Touches http://youtu.be/VCyKQAMUpmU


% Basic Graphical User Interface (GUI) without using GUIDE

% There are three basic areas to understand:
%   1. Layout    (how to position objects on the GUI)
%   2. Handles to Objects (used to modify object properties)
%   3. Callback functions (used by User Interface Objects)

%create a figure to house the GUI
figure

%create an annotation object
ellipse_position = [0.4 0.6 0.1 0.2];
ellipse_h = annotation('ellipse',ellipse_position,...
'facecolor', [1 0 0]);

%create an editable textbox object
edit_box_h = uicontrol('style','edit',...
'units', 'normalized',...
'position', [0.3 0.4 0.4 0.1]);

%create a "push button" user interface (UI) object
but_h = uicontrol('style', 'pushbutton',...
'string', 'Update Color',...
'units', 'normalized',...
'position', [0.3 0 0.4 0.2],...
'callback', {@eg_fun,edit_box_h, ellipse_h });

%Slider object to control ellipse size
uicontrol('style','Slider',...
'Min',0.5,'Max',2,'Value',1,...
'units','normalized',...
'position',[0.1    0.2    0.08    0.25],...
'callback',{@change_size,ellipse_h,ellipse_position });

uicontrol('Style','text',...
'units','normalized',...
'position',[0    0.45    0.2    0.1],...
'String','Ellipse Size')

%eg_fun code ---------------------------------------------------
%copy paste this code into a file called eg_fun.m
function eg_fun(object_handle, event)
disp('hi')

%updated eg_fun used to demonstrate passing variables
%copy paste this code into a file called eg_fun.m
function eg_fun(object_handle, event, edit_handle, ellipse_handle)
str_entered = get(edit_handle, 'string');

if strcmp(str_entered, 'red')
col_val = [1 0 0];
elseif strcmp(str_entered, 'green')
col_val = [0 1 0];
elseif strcmp(str_entered, 'blue')
col_val = [0 0 1];
else
col_val = [0 0  0];
end
set(ellipse_handle, 'facecolor', col_val);

%change_size code --------------------------------------------------
%copy paste this code into a file called change_size.m
function  change_size(objHandel, evt, annotation_handle, orig_pos)
slider_value = get(objHandel,'Value');
new_pos = orig_pos;
new_pos(3:4) = orig_pos(3:4)*slider_value;
set(annotation_handle, 'position', new_pos)


## DFT Windowing

This is code I used in a youtube video on DFT windowing – see http://youtu.be/PYd-pzaZpYE

% Tutorial on why DFT windowing works.
% See http://youtu.be/V_dxWuWw8yM for intro video

fs = 1000;
T = 1/fs;
N = 1000; % Create a signal with N samples = 1 second
n = 0:N-1; % sample numbers
t = n*T; % sampling times over a one second duration.
x = cos(2*pi*6.5*t);

subplot(2,1,1)
plot(n,x);
xlabel('Samples');
ylabel('Amplitude')
title('Time-Domain Signal to be analysed');

subplot(2,1,2)
xlabel('Bins');
ylabel('Magnitude')
title('First 1500 Bins of Magnitude Spectrum of Zero-Padded Signal');

% Plot magnitude spectrum against frequency in Hertz
fax = [0:100000-1]*fs/100000;
xlabel('Frequency (Hertz)');
ylabel('Magnitude')

%window using a hann window function
han_win = hanning(N)';
hold on

%The zero padded signal can be considered as being a sinusoid
% which is multiplied by a rectangular winddow. The folloiwing
% code illustrates this.
Ndisp = 5000; %number of time-domain samples to display
subplot(3,1,1)
title(['Zero Padded Signal (first ' num2str(Ndisp) ' samples)'])
ylabel('Amplitude')
xlabel('Samples');
t = [0:100000-1]*T; % sampling times for 100000 samples.
sinusoid = cos(2*pi*6.5*t);
subplot(3,1,2)
plot(sinusoid(1:Ndisp))
title(['Sinusoidal Signal (first ' num2str(Ndisp) ' samples)'])
ylabel('Amplitude')
xlabel('Samples');
subplot(3,1,3)
win = [ones(1, N) zeros(1, 99000)];
%win = [hanning(N)' zeros(1, 99000)];
plot(win(1:Ndisp))
title(['Rectangular Window (first ' num2str(Ndisp) ' samples)'])
ylabel('Amplitude')
xlabel('Samples');
ylim([-0.5 1.5])

%When you multiply signals in the time-domain it is equivalent to
% convolution in the frequency domain. Let's take a look at each of
% these three waveforms in the frequency domain
%update fax to  show negative frequencies
fax = [0:100000-1]*fs/100000;
fax(end-49998:end) = fax(end-49998:end)-1000;
fax(50001) = NaN;
subplot(3,1,1)
title(['Double Sided Magnitude Spectrum of  Zero Padded Signal'])
ylabel('Magnitude')
xlabel('Frequency (Hz)');

subplot(3,1,2)
plot(fax,abs(fft(sinusoid)))
title(['Double Sided Magnitude Spectrum of  Sinusoidal Signal'])
ylabel('Magnitude')
xlabel('Frequency (Hz)');

subplot(3,1,3)
plot(fax, abs(fft(win)))
title(['Double Sided Magnitude Spectrum of' ...
' Rectangular Window function '])
ylabel('Magnitude')
xlabel('Frequency (Hz)');

% Determine the magnitude spectrums with less zero padding
% Plot magnitude spectrum against frequency in Hz rather than bins
% Observe that the lower res spectrums are a sampled version of the
% higher res ones
X_mags4000 = abs(fft([x zeros(1, 3000)]));
X_mags2000 = abs(fft([x zeros(1, 1000)]));
X_mags1000 = abs(fft(x));

figure
fax = [0:4000-1]*fs/4000;
plot(fax(1:60), X_mags4000(1:60));
hold on
plot(fax(1:60), X_mags4000(1:60),'.');
set(gca,'XTick', fax(1:4:60))
xlabel('Frequency in Hertz')
ylabel('Magnitude')

fax = [0:2000-1]*fs/2000;
plot(fax(1:30), X_mags2000(1:30),'rx');

fax = [0:1000-1]*fs/1000;
plot(fax(1:15), X_mags1000(1:15),'go');

fax = [0:100000-1]*fs/100000;

%Each of the other magnitude spectrums are a subset or sampled
% version of the high resolution one.
%It should be appreciated from these examples that all DFT's can
% be thought of as being a sampled version of the DTFT (which
% produces a continuous spectrum that could be sampled at any
% frequency)
%Understanding the effects of windowing therefore requires an
% undertstanding of what happening for the high res case.

% Comparison of some common window functions
figure
N = 1000;
plot([-5 :.01: 5],hamming_freq(L/2-500:L/2+500),'r')
hold on
plot([-5 :.01: 5],rect_freq(L/2-500:L/2+500),'g')
plot([-5 :.01: 5],hann_freq(L/2-500:L/2+500),'c')
plot([-5 :.01: 5],blackman_freq(L/2-500:L/2+500),'k')
legend('Hamming', 'Rectangular','Hann','Blackman')
ylabel('Magnitude');



This is the code I used in the youtube tutorial on zero padding (see http://youtu.be/7yYJZfc5q-I).

% This tutorial attempts to explain how zero-padding works.
% You should already have some exposure to how the DFT works and
% be have seen the effects of zero paddding and windowing.
%

fs = 1000;
T = 1/fs;
N = 1000; % Create a signal with N samples = 1 second
n = 0:N-1; % sample numbers
t = n*T; % sampling times over a one second duration.
x = cos(2*pi*3*t);

%plot the magnitude spectrum
X_mags = abs(fft(x));

subplot(2,1,1)
plot(x)
xlabel('Samples');
ylabel('Amplitude')
title('Time-Domain Signal');

num_disp_bins = 15;
subplot(2,1,2)
plot([0:num_disp_bins-1], X_mags(1:num_disp_bins));
hold on
plot([0:num_disp_bins-1], X_mags(1:num_disp_bins),'k.');
hold off
xlabel('Frequency Bins');
ylabel('Magnitude');
title('Frequency Content Magnitudes');

%Illustrate Basis functions
plot(n,x,'k')
xlabel('Samples');
ylabel('Amplitude')
for k = 0 : 20
plot(x,'kx')
hold on
cos_basis = cos(2*pi*k*n/N);
sin_basis = sin(2*pi*k*n/N);

plot(n,cos_basis,'r')
plot(n,sin_basis,'g')
title({['Signal being anlaysed (black) with basis functions'...
' that have ' num2str(k) ' cyles over '  ...
num2str(N) ' samples']...
['Cosine basis function shown in red, Sine in green']})
hold off
pause
end



## DFT Basis Function Plots

Just posting this as I had request for code I used in creating a youtube video on the DFT

%This a script I used when creating the images for the DFT explained youtube video http://youtu.be/B2iUDBZzBpY
%
% It shows the the result of multiplying a signal being analysed by a set of sinusoidal basis functions (which is how the DFT works)
% The figure produced was edit further before using it in the video
close all;
N = 8; %signals are only 8 samples long
n = 0:N-1;

%specify the magnitudes and phases of the sinusoids (integer cyles over N samples)
mags = [ 0 0 5 0 0 0 ]; %change this to include more sinusoids i.e. change 0 values to non zero
phases = [0 0 0 0 0 0]; %phase values for each of the sinusoids

y_range = [-5 5];
interp_val = 100;
n_vis = 0 : (N)*interp_val-1;
x_vis = zeros(1, (N)*interp_val);
x = zeros(1, N);

for k = 0 : N/2
x = x + mags(k+1)*cos(2*pi*(k)*n/N + phases(k+1)); %samples of signal
x_vis = x_vis + mags(k+1)*cos(2*pi*(k)*n_vis/(N*interp_val) + phases(k+1)); %interpolated sinusoid to aid with visualisation
end

num_rows = N/2+2 +N/2; %number of rows in the subplot
subplot(num_rows,2,1)
plot(x_vis);
hold on
plot(n*interp_val, x,'r.')
hold off

set(gca, 'XTick', 0:interp_val:(N)*interp_val)
set(gca, 'XTickLabel', '')
set(gca, 'Box', 'off')
set(gca, 'YTick', y_range)
ylim(y_range)
%for each of the basis functions
for k = 0 : N/2
basis_function = cos(2*pi*k*n/N); %samples of the cosine basis function
basis_function_vis = cos(2*pi*k*n_vis/(N*interp_val)); %an interpolated version to aid with visualisation

mult_sig = basis_function.*x;
mult_sig_vis = basis_function_vis.*x_vis;

basis_function2 = sin(2*pi*k*n/N); %sine wave basis function
basis_function_vis2 = sin(2*pi*k*n_vis/(N*interp_val));

mult_sig2 = basis_function2.*x;
mult_sig_vis2 = basis_function_vis2.*x_vis;

subplot(num_rows,2,(2*k+2)*2-1)
plot(basis_function_vis)
hold on
plot(n*interp_val,basis_function ,'r.')
hold off
set(gca, 'XTick', 0:interp_val:(N)*interp_val)
label = '';
set(gca, 'XTickLabel', label)
set(gca, 'Box', 'off')
if(k > 0)
subplot(num_rows,2,(2*k+1)*2-1)
plot(basis_function_vis2)
hold on
plot(n*interp_val,basis_function2 ,'r.')
hold off
set(gca, 'XTick', 0:interp_val:(N)*interp_val)
label = '';
set(gca, 'XTickLabel', label)
set(gca, 'Box', 'off')
end

subplot(num_rows,2,(2*k+2)*2)
hold on
plot(mult_sig_vis)
plot(zeros(1, N*interp_val),'k')
plot(n*interp_val,mult_sig ,'r.')
hold off
set(gca, 'XTick', 0:interp_val:(N)*interp_val)
set(gca, 'YTick', y_range)
set(gca, 'XTickLabel', label)
set(gca, 'Box', 'off')
ylim(y_range)
xlim([0 N*interp_val])
ylabel(num2str(round(sum(mult_sig)*100)/100))
if(k > 0)
subplot(num_rows,2,(2*k+1)*2)
hold on
plot(mult_sig_vis2)
plot(zeros(1, N*interp_val),'k')
plot(n*interp_val,mult_sig2 ,'r.')
hold off
set(gca, 'XTick', 0:interp_val:(N)*interp_val)
set(gca, 'XTickLabel', label)
set(gca, 'Box', 'off')
set(gca, 'YTick', y_range)
ylim(y_range)
xlim([0 N*interp_val])
ylabel(num2str(round(sum(mult_sig2)*100)/100))
end
end



## Filter Design Using Matlab Demo

This is the code I use in my youtube videos on filtering (ddorran channel)

%generate the noisy signal which will be filtered
x= cos(2*pi*12*[0:0.001:1.23]);
x(end) = [];
[b a] = butter(2,[0.6 0.7],'bandpass');
filtered_noise = filter(b,a,randn(1, length(x)*2));
x = (x + 0.5*filtered_noise(500:500+length(x)-1))/length(x)*2;
% This script is available at https://dadorran.wordpress.com search for
% filtering matlab demo
plot(x)
title('Noisy signal')
xlabel('Samples');
ylabel('Amplitude')

%plot magnitude spectrum of the signal
X_mags = abs(fft(x));
figure(10)
plot(X_mags)
xlabel('DFT Bins')
ylabel('Magnitude')

%plot first half of DFT (normalised frequency)
num_bins = length(X_mags);
plot([0:1/(num_bins/2 -1):1], X_mags(1:num_bins/2))
ylabel('Magnitude')

%design a second order filter using a butterworth design technique
[b a] = butter(2, 0.3, 'low')

%Show the signal flow diagram of the filter/system
create_signal_flow(b,a)

%plot the frequency response (normalised frequency)
H = freqz(b,a, floor(num_bins/2));
hold on
plot([0:1/(num_bins/2 -1):1], abs(H),'r');

%filter the signal using the b and a coefficients obtained from
%the butter filter design function
x_filtered = filter(b,a,x);

%plot the filtered signal
figure(2)
plot(x_filtered,'r')
title('Filtered Signal - Using Second Order Butterworth')
xlabel('Samples');
ylabel('Amplitude')

%Redesign the filter using a higher order filter
[b2 a2] = butter(20, 0.3, 'low')

%Plot the magnitude spectrum and compare with lower order filter
H2 = freqz(b2,a2, floor(num_bins/2));
figure(10)
hold on
plot([0:1/(num_bins/2 -1):1], abs(H2),'g');

%filter the noisy signal and plot the result
x_filtered2 = filter(b2,a2,x);
figure(3)
plot(x_filtered2,'g')
title('Filtered Signal - Using 20th Order Butterworth')
xlabel('Samples');
ylabel('Amplitude')

%Use a band reject filter in place of a low pass filter
[b_stop a_stop] = butter(20, [0.5 0.8], 'stop');

%plot the magnitude spectrum
H_stopband = freqz(b_stop,a_stop, floor(num_bins/2));
figure(10)
hold on
plot([0:1/(num_bins/2 -1):1], abs(H_stopband),'c');

%plot filtered signal
x_filtered_stop = filter(b_stop,a_stop,x);
figure(4);
plot(x_filtered_stop,'c')
title('Filtered Signal - Using Stopband')
xlabel('Samples');
ylabel('Amplitude')

%Use matlabs built-in buttord function to get the optimum order to meet a specification
[N Wn] = buttord(0.1, 0.5, 5, 40)

%use the N and Wn values obtained above to design the filter in the usual way
[b3 a3] = butter(N, Wn, 'low');

%plot the magnitude spectrum
H3 = freqz(b3,a3, floor(num_bins/2));
figure(10);
hold on
plot([0:1/(num_bins/2 -1):1], abs(H2),'k');
figure(10)

%filter the signal and plot the ouput of the filter
x_filtered3 = filter(b3,a3,x);
figure(5);
plot(x_filtered3,'k')
title(['Filtered Signal - Using ' num2str(N) ' th Order Butterworth'])
xlabel('Samples');
ylabel('Amplitude')
% comparison with other filter design techniques (chebyshev and elliptical)
[b_butter a_butter] = butter(4, 0.2, 'low');
H_butter = freqz(b_butter, a_butter);

[b_cheby a_cheby] = cheby1(4, 0.5, 0.2, 'low');
H_cheby = freqz(b_cheby, a_cheby);

[b_ellip a_ellip] = ellip(4, 0.5, 40, 0.2, 'low');
H_ellip = freqz(b_ellip, a_ellip);

%plot each filter to compare
figure(11)
norm_freq_axis = [0:1/(512 -1):1];
plot(norm_freq_axis, abs(H_butter))
hold on
plot(norm_freq_axis, abs(H_cheby),'r')
plot(norm_freq_axis, abs(H_ellip),'g')
legend('Butterworth', 'Chebyshev', 'Elliptical')
xlabel('Normalised Frequency');
ylabel('Magnitude')

%plot in dB for verification that spec is met
figure(12);
plot(norm_freq_axis, 20*log10(abs(H_butter)))
hold on
plot(norm_freq_axis, 20*log10(abs(H_cheby)),'r')
plot(norm_freq_axis, 20*log10(abs(H_ellip)),'g')
legend('Butterworth', 'Chebyshev', 'Elliptical')
xlabel('Normalised Frequency ');
ylabel('Magnitude (dB)')



## matlab complex exponential visualisation

% Used to help visualise complex exponentials as shown in my youtube video (http://youtu.be/K_C7htSXORY)
% by david.dorran@gmail.com

function complex_exponential_gui(varargin)

global x;global t;global w;global real_circle_pos; global imag_circle_pos;
global argand_circle_pos; global T; global freq;

diameter = 0.15; %diameter of argand circle
border = 0.1; %gap between objects

action = '';
if(nargin > 0)
action = varargin{1}
end
if(~freq)
freq = 'positive'
end
if(strcmp(action, 'toggle'))
if(strcmp(freq,'positive'))
freq = 'negative';
else
freq = 'positive';
end
end
if(nargin ==0 || strcmp(action, 'toggle'))
close all
w = pi/3;

T = 0.1;
t = 0:T:20.3;

if(~strcmp(freq, 'negative'))
x = exp(j*w*t);
else
x = exp(-j*w*t);
end

uicontrol('style','pushbutton',...
'string','start',...
'fontsize',8,...
'units','normalized',...
'position',[.15 .46 .18 .04],...
'callback','complex_exponential_gui(''play'')');
uicontrol('style','pushbutton',...
'string',' ',...
'fontsize',8,...
'units','normalized',...
'position',[.1 .1 .01 .04],...
'callback','complex_exponential_gui(''toggle'')');

if(~strcmp(freq, 'negative'))
annotation('textbox',[0.3 0.75 0.1 0.1],'string','e^{j{\omega}t} = cos({\omega}t)+jsin({\omega}t)','linestyle','none','fontsize',14)

annotation('textbox',[0.7 0.3 0.1 0.1],'string','jsin({\omega}t)','linestyle','none','fontsize',10)
else
annotation('textbox',[0.3 0.75 0.1 0.1],'string','e^{j-{\omega}t} = cos(-{\omega}t)+jsin(-{\omega}t)','linestyle','none','fontsize',14)
annotation('textbox',[0.3 0.68 0.1 0.1],'string','e^{-j{\omega}t} = cos({\omega}t)-jsin({\omega}t)','linestyle','none','fontsize',14)
annotation('textbox',[0.7 0.3 0.1 0.1],'string','-jsin({\omega}t)','linestyle','none','fontsize',10)
end
annotation('textbox',[0.3 0.1 0.1 0.1],'string','cos({\omega}t)','linestyle','none','fontsize',10)
annotation('textbox',[0.64 0.35 0.1 0.1],'string','-j','linestyle','none','fontsize',10)
annotation('textbox',[0.64 0.5 0.1 0.1],'string','j','linestyle','none','fontsize',10)
annotation('textbox',[0.4 0.25 0.1 0.1],'string','-1','linestyle','none','fontsize',10)
annotation('textbox',[0.56 0.25 0.1 0.1],'string','1','linestyle','none','fontsize',10)
subplot(3,1,3)
plot(real(x));
hold on
real_circle_pos = plot(1, real(x(1)),'r.', 'markersize',20);
plot([0 length(x)*1.3], [ 0 0 ] ,'k','linewidth',2)
%plot([0 1], [0 1.5] ,'k','linewidth',2)

plot([length(x)*1.2 length(x)*1.3], [ 0.1 0 ] ,'k','linewidth',2)
plot([length(x)*1.2 length(x)*1.3], [ -0.1 0 ] ,'k','linewidth',2)
set(gca, 'visible','off', 'position',[0.5-diameter/2 border-.1 diameter 0.5-diameter-border+0.05],'CameraUpVector',[-1 0 0])
hold off

subplot(3,1,2)
plot(imag(x));
hold on
imag_circle_pos = plot(1, imag(x(1)),'r.', 'markersize',20);
plot([0 length(x)*1.3], [ 0 0 ] ,'k','linewidth',2)
plot([length(x)*1.2 length(x)*1.3], [ 0.1 0 ] ,'k','linewidth',2)
plot([length(x)*1.2 length(x)*1.3], [ -0.1 0 ] ,'k','linewidth',2)
set(gca, 'visible','off', 'position',[0.5+diameter/2+border 0.5-diameter/2  0.5-diameter-border diameter])

subplot(3,1,1)
plot3(t, -real(x), imag(x));
hold on
plot3([0 T], [1.4 -1.4], [0 0],'k','linewidth',2)
plot3([0 T], [1.4 1.2], [0 0.1],'k','linewidth',2)
plot3([0 T], [1.4 1.2], [0 -0.1],'k','linewidth',2)
plot3([0 T], [-1.4 -1.2], [0 0.1],'k','linewidth',2)
plot3([0 T], [-1.4 -1.2], [0 -0.1],'k','linewidth',2)

plot3([0 T], [0 0] , [1.4 -1.4],'k','linewidth',2)
plot3([0 T], [0 0.1],[1.4 1.2], 'k','linewidth',2)
plot3([0 T], [0 -0.1],[1.4 1.2], 'k','linewidth',2)
plot3([0 T], [0 0.1],[-1.4 -1.2], 'k','linewidth',2)
plot3([0 T],  [0 -0.1],[-1.4 -1.2],'k','linewidth',2)

plot3([0 t(end)*1.4],  [0 0],[0 0 ],'k','linewidth',2)
plot3([t(end)*1.3 t(end)*1.4],  [0 0],[0.1 0 ],'k','linewidth',2)
plot3([t(end)*1.3 t(end)*1.4],  [0 0],[-0.1 0 ],'k','linewidth',2)

argand_circle_pos = plot3(0, -real(x(1)), imag(x(1)),'g.', 'markersize',20);

set(gca,'visible','off',  'position',[0.5-diameter 0.5-diameter  diameter*2 diameter*2],'cameraposition',[1 0 0])

set(gca,'ytick', [-1 0 1])
set(gca,'yticklabel', [-1 0 1])
set(gca,'xtick', [-1 0 1])
set(gcf, 'ToolBar', 'figure');
end
if(strcmp(action,'play') || strcmp(action,'pause'))
offset = 0;
one_pause = 0;
two_pause = 0;
for k = 1:length(t)-1
set(real_circle_pos, 'Xdata', k,'Ydata', real(x(k)));
set(imag_circle_pos, 'Xdata', k,'Ydata', imag(x(k)));
set(argand_circle_pos, 'Xdata',t(k) ,'Ydata', -real(x(k)), 'Zdata',imag(x(k)));

if(imag(x(k)) > imag(x(k+1)) && one_pause == 0 && strcmp(freq,'positive'))
pause(10)
offset = offset+ 10;
one_pause = 1;
end
if(imag(x(k)) > -1/sqrt(2) && imag(x(k)) < 0 && real(x(k)) > 0 && two_pause == 0 && strcmp(freq,'positive'))
pause(7)
offset = offset+ 7;
two_pause = 1;
end
drawnow;
pause(T)
end
end


## real-time audio visualisation/plotting in matlab


% Plot audio from soundcard captured at Fs Hz with frames of length N (with thanks to Richard Hayse)

clear all;
Fs=11025;
N=2048;t=(0:N-1)/Fs;

while 1
y = wavrecord(N,Fs,2);

plot(y);
axis( [0 N -1 1]);
drawnow;
end