### Archive

Archive for the ‘matlab code’ Category

## Heartrate (BPM) Example Matlab Code

May 22, 2014 Leave a comment

This is the code I used in my youtube video at http://youtu.be/3tdumuwHgxc

% program to determine the BPM of an ECG signal

% count the dominant peaks in the signal (these correspond to heart beats)
% - peaks are defined to be sampels greater than their two nearest neighbours and greater than 1

beat_count = 0;
for k = 2 : length(sig)-1
if(sig(k) &gt; sig(k-1) &amp; sig(k) &gt; sig(k+1) &amp; sig(k) &gt; 1)
%k
%disp('Prominant peak found');
beat_count = beat_count + 1;
end
end

% Divide the beats counted by the signal duration (in minutes)
fs = 100;
N = length(sig);
duration_in_seconds = N/fs;
duration_in_minutes = duration_in_seconds/60;
BPM_avg = beat_count/duration_in_minutes;


## cross correlation demo

April 25, 2014 2 comments
% A demonstration of cross correlation in action.
% Hit the space bar during the demo to execute
%
% https://dadorran.wordpress.com/2014/04/25/cross-correlation-demo/
clc;close all
a = [0.1 0.2  -0.1 4.1 -2 1.5 0 ];
b = [0.1 4 -2.2 1.6 0.1 0.1 0.2];

len = length(a);
if(len ~= length(b))
error('vectors supplied must be the same length');
end
figure
set(gcf, 'position', [ 285   347   642   367]);

max_amp = max([max(a) max(b)]);
min_amp = min([min(a) min(b)]);

plot_h = 0.25;
text_h = 0.1;
ax1 = subplot(2,1,1);
pl1_line = plot(a);
labels1 = text([1:len], a , num2str(a'), 'VerticalAlignment','bottom', ...
'HorizontalAlignment','right','fontsize',8);
hold on; pl1_dot = plot(a,'r.');
xlim([1 len])
ylim([min_amp max_amp])

set(ax1,'position', [(1/3) 0.95-plot_h (1/3) plot_h])
set(ax1,'visible','off')

ax2 = subplot(2,1,2);
pl2_line = plot(b);
labels2 = text([1:len], b , num2str(b'), 'VerticalAlignment','bottom', ...
'HorizontalAlignment','right','fontsize',8);
hold on; pl2_dot = plot(b,'r.');
xlim([1 len])
ylim([min_amp max_amp])
set(ax2,'visible','off')
set(ax2,'position', [(1/3) 0.9-plot_h*2 (1/3) plot_h])

str = '';
for k = 1: len
str = [str '(' num2str(a(k)) ')(' num2str(b(k)) ') + '];
end
str(end-1)  = '=';
str = [str num2str(sum(a.*b))];
r_ba = xcorr(a,b);

corr_calc_text = annotation('textbox',  [0 0.85-plot_h*2-text_h 1 text_h], 'linestyle','none','horizontalalignment','center' ,'string', {'correlation at zero lag is ' str}, 'fontsize', 8);

annotation('textbox',  [0.5 0.8-plot_h*2-text_h*2 1 text_h], 'linestyle','none','horizontalalignment','left' ,'string', sprintf('%.2f',r_ba(len)),'color','red', 'fontsize', 8);

pause

x_inc= (1/3)/(len-1);
for k = 1:len-1

str = '';
for m = 1: len-k
str = [str '(' num2str(a(m+k)) ')(' num2str(b(m)) ') + '];
end
str(end-1)  = '=';
str = [str num2str(r_ba(len+k))];

set(corr_calc_text,'string', {['correlation at lag of ' num2str(k) ' is '] str}, 'fontsize', 8);

set(ax2,'position', [(1/3)+k*x_inc 0.9-plot_h*2 (1/3) plot_h])

annotation('textbox',  [0.5+x_inc*k 0.8-plot_h*2-text_h*2 1 text_h], 'linestyle','none','horizontalalignment','left' ,'string', sprintf('%.2f',r_ba(len+k)),'color','red', 'fontsize', 8);
if(k ==1)
pause
annotation('textbox',  [0.5 0.01 1 text_h], 'linestyle','none','horizontalalignment','left' ,'string', ['   0'] ,'color','blue', 'fontsize', 8);
annotation('textbox',  [0.001 0.01 1 text_h], 'linestyle','none','horizontalalignment','left' ,'string', ['Lag:'] ,'color','blue');
annotation('textbox',  [0.5+x_inc*(len) 0.8-plot_h*2-text_h*2 1 text_h], 'linestyle','none','horizontalalignment','left' ,'string', ']' ,'color','red');
annotation('textbox',  [0.5-x_inc*(len-1)-x_inc/2 0.8-plot_h*2-text_h*2 1 text_h], 'linestyle','none','horizontalalignment','left' ,'string', '[' ,'color','red');
annotation('textbox',  [0.001 0.8-plot_h*2-text_h*2 1 text_h], 'linestyle','none','verticalalignment','middle','horizontalalignment','left' ,'string', {'Correlation' 'Sequence:'} ,'color','red');

end
annotation('textbox',  [0.5+x_inc*k 0.01 1 text_h], 'linestyle','none','horizontalalignment','left' ,'string', ['   ' num2str(k)] ,'color','blue', 'fontsize', 8);

pause
end

for k = 1:len-1

str = '';
for m = 1: len-k
str = [str '(' num2str(a(m)) ')(' num2str(b(m+k)) ') + '];
end
str(end-1)  = '=';
str = [str num2str(r_ba(len-k))];

set(corr_calc_text,'string', {['correlation at lag of ' num2str(-1*k) ' is '] str}, 'fontsize', 8);

set(ax2,'position', [(1/3)-k*x_inc 0.9-plot_h*2 (1/3) plot_h])
annotation('textbox',  [0.5-x_inc*k 0.8-plot_h*2-text_h*2 1 text_h], 'linestyle','none','horizontalalignment','left' ,'string', sprintf('%.2f',r_ba(len-k)),'color','red', 'fontsize', 8);
annotation('textbox',  [0.5-x_inc*k 0.01 1 text_h], 'linestyle','none','horizontalalignment','left' ,'string', ['   ' num2str(k*-1)] ,'color','blue', 'fontsize', 8);

pause
end

% Uncomment the next two lines if you would like to see a plot of the
% correlation sequence
% [corr_seq lags] =  xcorr(a,b);
% plot(lags,corr_seq)
% xlabel('lags');ylabel('correlation measure');



## Matlab Fourier Demo

March 27, 2014 Leave a comment
% illustration of Fourier Theory using plots
%
% usage : fourier_demonstration(1) % a well defined signal
%          fourier_demonstration(2)  % a segment of speech signal (download from https://www.dropbox.com/s/bw4dpf93xxz1lyb/speech_seg.wav)
%          fourier_demonstration(3)  % a SQUARE WAVE
%          fourier_demonstration(4)  % an impulse
function fourier_demonstration(num)
if(num==1)

%sum of sinusoids as input
t = 0:1/16000:1-1/16000;
s1 = cos(2*pi*1*t);%cosw(16000,1, 16000, 0);
s2 = cos(2*pi*5*t)*5;%cosw(16000,5, 16000, 0)*5;
s3 = cos(2*pi*10*t + pi/2+0.56)*3;%cosw(16000,10, 16000, pi/2+0.56)*3;
s4 = cos(2*pi*8*t+pi)*3;%cosw(16000,8, 16000, pi)*3;
s5 = cos(2*pi*12*t+pi/2.2)*3;%cosw(16000,12, 16000, pi/2.2)*3;
s6 = cos(2*pi*3*t+pi/2+0.4);%cosw(16000,3, 16000, pi/2+0.4)*3;
sig = s1+s2+s3+s2+s3+s4;
elseif(num==2)
sig = wavread('speech_seg.wav')';
elseif(num==3)
sig = [zeros(1,100) ones(1,100) zeros(1,100) ones(1,100) zeros(1,100) ones(1,100) zeros(1,100) ones(1,100) ];
sig = (sig-0.5)*2;
else
sig = [zeros(1, 200) 1 zeros(1,200)];

end
sig = sig - mean(sig);

N = length(sig);
ft = fft(sig);
ft(round(length(ft)/2)-2:end) = [];
mags = abs(ft);
phases = angle(ft);
dc_mag = mags(1);
mags(1) = [];
phases(1) = [];

[sorted_mags sorted_indices] = sort(mags,2);
% sort in decending order
sorted_mags = fliplr(sorted_mags);
sorted_indices = fliplr(sorted_indices);
sorted_phases = phases(sorted_indices);

synth_op = zeros(1, N);
sig = sig -dc_mag;
n = [0:N-1]; % sample numbers

ylim_max = max([ sorted_mags(1)/N*2 sig])*1.05;
ylim_min = min([ -sorted_mags(1)/N*2 sig])*1.05;
ylims = [ylim_min ylim_max ];
subplot(3,1,1)
plot(sig)
set(gca,'Xticklabel','', 'YLim', ylims)
set(gca,'Xticklabel','')
ylabel('Amplitude')
xlabel('Time')
title('Example Signal');
pause
colors = 'rgkbm';
significant_freqs = find(sorted_mags > max(sorted_mags/100));
freq_vals_to_display = max(sorted_indices(1:length(significant_freqs))) +1 ;
length(mags)
for k  = 1: length(mags)
omega = 2*pi*(sorted_indices(k))/N;
sinusoid = cos(n*omega+sorted_phases(k))*sorted_mags(k)/N*2 ;
if(sorted_mags(k) < 10^-6)
break
end
synth_op = synth_op + sinusoid;
subplot(3,1,1)
plot(sig)
set(gca,'Xticklabel','', 'YLim', ylims)
hold on
plot(synth_op,'r')
set(gca,'Xticklabel','')
ylabel('Amplitude')
xlabel('Time')
if k ==1
title('Example signal and sinusoid shown in lower plot');
else
title(['Example signal and ' num2str(k) ' sinusoids shown in middle plot added together.']);
end
hold off
subplot(3,1,2)
hold on
plot(sinusoid,colors(rem(k,5)+1))
set(gca,'Xticklabel','')
%set(gca, 'Ylim', [-max(mags)/N*2 max(mags)/N*2])
set(gca, 'Ylim', ylims)
ylabel('Amplitude')
xlabel('Time')
subplot(3,1,3);
ft_mag_vals = ones(1, freq_vals_to_display)*NaN;
sorted_indices(k)
if( sorted_indices(k) < freq_vals_to_display)
ft_mag_vals(sorted_indices(k)) = sorted_mags(k)/N*2;
end
hold on

%stem([0:freq_vals_to_display],[NaN ft_mag_vals],'^')
hold off
ylabel('Magnitude')
xlabel('Normalised Frequency (1/(periods displayed))')
pause
%plot(sinusoid,colors(rem(k,6)+1))
end
close all


## Correlation implemented in matlab

February 24, 2014 Leave a comment

This is the code shown at the end of http://youtu.be/_r_fDlM0Dx0

% Code to calculate the correlation measurement between to
% two signals of length N samples
% Available at https://dadorran.wordpress.com

a = [2.1 3.8 -3.6 4.1 -2.9];
b = [-1.1 3.2 -3.6 2.2 -4.2];
N = length(a); % must be the same as length(b)

corr_measure1 = 0;
for n = 1:N
mult_result = a(n)*b(n);
corr_measure1 = corr_measure1 + mult_result;
end

disp(['The correlation measurement is ' num2str(corr_measure1)])

% Alternative concise implementation
corrmeasure2 = sum(a.*b)



## Plotting Frequency Spectrum using Matlab

February 20, 2014 8 comments

This post shows a variety of ways of how to plot the magnitude frequency content of a discrete signal using matlab.

### Load Example Data

% 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

## 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



## create_tf_equation

February 3, 2014 Leave a comment

A function which creates an annotation based plot that shows the transfer function of a system given the systems b and a coefficients.

% This function creates a transfer function given a set of b and a
% coefficients
%
% example usage:
%             create_tf_equation([0.23 0.5], [1 -0.34 -0.23432]);
%
function create_tf_equation(b,a, varargin)

%To Do: cater for s-domain
if(~length(b))
error('b must have some values')
end
if(~length(a))
error('a must have some values')
end

b_str = '';
for k = 1:length(b)
if( k ==1)
z_var = '';
sign_var = '';
else
z_var = ['z^{' num2str((k-1)*-1) '}'];
if(b(k) < 0)
sign_var = ' - ';
b(k) = b(k)*-1;
else
sign_var = ' + ';
end
end
b_str = [ b_str sign_var num2str(b(k)) z_var];
end

a_str = '';
for k = 1:length(a)
if( k ==1)
z_var = '';
sign_var = '';
else
z_var = ['z^{' num2str((k-1)*-1 ) '}'];
if(a(k) < 0)
sign_var = ' - ';
a(k) = a(k)*-1;
else
sign_var = ' + ';
end
end
a_str = [ a_str sign_var num2str(a(k)) z_var];
end

if(length(a) == 1)
H_str = b_str;
else

H_str = ['\frac{' b_str   '}{' a_str '}'];
end
fh = figure
text('Interpreter','latex',...
'String',['$$H(z) = ' H_str '$$'],...
'Position',[0,0.2],...
'FontSize',14);
set(gca, 'visible', 'off')
fwidth = round(100 + max([length(a_str) length(b_str)])*8.5);
set(fh, 'position', [232   246 fwidth      54]);


Categories: matlab code, Uncategorized

## freq_plot.m

January 21, 2014 Leave a comment
%a function which plots the magnitude spectrum of a signal from DC (0Hz) up to
%max_freq.
%
%   example:
%           fs = 10000;
%           t = 0:1/fs: 4;
%           sig = 3*cos(2*pi*100*t) + 3*cos(2*pi*20*t);
%           freq_plot(sig, 120, fs)
%
% by david.dorran@gmail.com, Feb 2012
function freq_plot(sig, max_freq, fs)
if(max_freq>fs/2)
error('Max freq must be less than Nyquist frequency');
end
X = abs(fft(sig));
N = length(sig);

f = 0:fs/length(sig):fs;
num_bins = find(f> max_freq);
X = X(1:num_bins(1));
f = f(1:num_bins(1));
plot(f, X/(N)*2+0.0001)
xlabel('frequency (Hz)')
ylabel('Magnitude');

Categories: matlab code, Uncategorized

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