Archive

Archive for the ‘youtube demo code’ Category

plotting the frequency content of a signal

February 17, 2014 3 comments

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
% download data from
%https://www.dropbox.com/s/n4kpd4pp2u3v9nf/tremor_analysis.txt
signal = load('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];
fax_rads_sample = (bin_vals/N)*2*pi;
N_2 = ceil(N/2);
plot(fax_rads_sample(1:N_2), X_mags(1:N_2))
xlabel('Frequency (radians per sample)')
ylabel('Magnitude');
title('Single-sided Magnitude spectrum (rads/sample)');
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

January 14, 2014 2 comments

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
% Available at https://dadorran.wordpress.com

% 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

January 7, 2014 Leave a comment

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
% Code available at https://dadorran.wordpress.com. 

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);

x_zpad = [x zeros(1, 99000)]; %zero pad
X_ZPAD_mags= abs(fft(x_zpad));

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

subplot(2,1,2)
bins = 0:length(x_zpad)-1;
plot(bins(1:1500), X_ZPAD_mags(1:1500))
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;
plot(fax(1:1500), X_ZPAD_mags(1:1500));
title('Magnitude Spectrum of Zero-Padded Signal');
xlabel('Frequency (Hertz)');
ylabel('Magnitude')

%window using a hann window function
han_win = hanning(N)';
x_zpad_win = [x.*han_win zeros(1, 99000)]; %zero pad
X_ZPAD_WIN_mags= abs(fft(x_zpad_win));
hold on
plot(fax(1:1500), X_ZPAD_WIN_mags(1:1500),'r');


%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)
plot(x_zpad(1:Ndisp))
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)
plot(fax,abs(fft(x_zpad)))
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;
plot(fax(1:1500), X_ZPAD_mags(1:1500),'k');

%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;
zpad_len = 99000;
L = N+zpad_len;
zpad = zeros(1,zpad_len);
hamming_freq = fftshift(abs(fft([hamming(N)', zpad])))/(N);
plot([-5 :.01: 5],hamming_freq(L/2-500:L/2+500),'r')
hold on
rect_freq = fftshift(abs(fft([rectwin(N)', zpad])))/(N);
plot([-5 :.01: 5],rect_freq(L/2-500:L/2+500),'g')
hann_freq = fftshift(abs(fft([hanning(N)', zpad])))/(N);
plot([-5 :.01: 5],hann_freq(L/2-500:L/2+500),'c')
blackman_freq = fftshift(abs(fft([blackman(N)', zpad])))/(N);
plot([-5 :.01: 5],blackman_freq(L/2-500:L/2+500),'k')
legend('Hamming', 'Rectangular','Hann','Blackman')
ylabel('Magnitude');
xlabel('Bin width (non zero padded)');


Zero Padding Tutorial M-file

December 20, 2013 Leave a comment

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

October 23, 2013 Leave a comment

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

October 18, 2013 8 comments

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))
xlabel('Normalised frequency (\pi rads/sample)')
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
%this function available from https://dadorran.wordpress.com
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')
% script available from https://dadorran.wordpress.com
% 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

September 13, 2013 Leave a comment
% 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

August 31, 2013 Leave a comment

% 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