### Archive

Archive for the ‘youtube demo code’ Category

## FFT based FIR filter design

%% This script deals issues related with 'fft' based FIR filter design
% See youtube video at https://youtu.be/R4RpNG_Botk

% GO TO lines 115 to 133 if you want to skip the preamble

%% First show fft of finite number of impulse response samples  lie on the continuous spectrum
%Design a band reject filter using buit-in fir1 function
b = fir1(22, [0.3913 0.6522],'stop'); %cutoff frequency is (0.3, 0.6 x pi) radians/sample = 0.3,0.6 x fs/2;
h = b; % The non-zero values of the impulse response of any fir filter
%   is equal to the b coefficients
H = fft(h); % Frequency response bin values
bin_freqs = [0:length(H)-1]/length(H)*2*pi;

%% Create the continuous spectrum
% Going to do this three ways - just for the purpose of demonstration.
% The mathematical formula for the frequency response of an mth order FIR is:
%   H(w) = b0 + b1*e^-jw + b2*e^-2jw + b3*e^-2*jw + ... ... + bm*e^(-m*jw)
%
%   This equation comes from the fact that H(w) = H(z) when z= e^jw i.e. when H(z) is
%   evaluated around the unit circle (see https://youtu.be/esZ_6n-qHuU )
%  The transfer function H(z) = b0*Z^0+b1*z^-1+b2*z^-2+b3*z-3+ ... +bmz^-m

% We'll create a 'continuous' spectrum by evaluating this equation for a
% "large" number of frequencies, w.

% Two alternatives to creating a 'continuous' spectrum are to zero pad the
% b coefficients by a "large" amount prior to taking the fft; and to use the built in freqz
% function. I'll use all three approaches here for comparison purposes.

% METHOD 1 - zero-pad method
H_cont_1 = abs(fft([h zeros(1,10000)])); % the value of 10000 could be any relatively large number i.e. large enough to capture the spectral detail of the frequency response
w_cont_1 = [0:length(H_cont_1)-1]/length(H_cont_1)*2*pi;

% METHOD 2 - freqz method
[H_cont_2_complex w_cont_2] = freqz(b,1,length(H_cont_1), 'whole'); % the second parameter of freqz is the a coefficients
H_cont_2 = abs(H_cont_2_complex);

% METHOD 3 -evaluate H(z) around the unit circle method (see
%https://youtu.be/esZ_6n-qHuU)
H_cont_3_complex = [];
k = 0;
w_cont_3 = [0:length(H_cont_1)-1]/length(H_cont_1)*2*pi;
for w = w_cont_3
k = k + 1;
H_cont_3_complex(k) = b(1);
for ind = 2:length(b)
H_cont_3_complex(k) = H_cont_3_complex(k) + b(ind)*exp(-1i*w*(ind-1));
end
end

H_cont_3= abs(H_cont_3_complex);

%% plot 'continuous' frequency response and overlay 'sampled' bin values
figure
plot(w_cont_1, H_cont_1,'LineWidth', 6)
hold on
plot(w_cont_2, H_cont_2,'LineWidth', 3)
plot(w_cont_3, H_cont_3)
plot(bin_freqs, abs(H),'g.','MarkerSize', 16)
ylabel('Magnitude')
legend('Cont. Freq Resp Method 1', 'Cont. Freq Resp Method 2', 'Cont. Freq Resp Method 3', 'Bin Values')
title('Frequency response of filter using a built-in filter design technique')
set(gcf,'Position',[22 416 613 420]);

%% Now let's try and design a filter by starting with the desired frequency response
% First we'll look at the band pass filter
H_desired = [0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0];
freq_bins_desired = [0:length(H_desired)-1]/length(H_desired)*2*pi;
b_ifft = ifft(H_desired); %b coefficients of newly designed filter
h_ifft = b_ifft;

% Now compute the continuous spectrum using any of the methods above. (Zero
% padding b_ifft prior to taking the fft is the most straightforward
% appraoch - I think!)
H_cont_desired = abs(fft([b_ifft zeros(1,1000)]));
w_cont_desired = [0:length(H_cont_desired)-1]/length(H_cont_desired)*2*pi;

figure
plot(w_cont_desired, H_cont_desired);
hold on
plot(freq_bins_desired, H_desired,'r.','MarkerSize', 12)

ylabel('Magnitude')
legend('Continuous Frequency Response','DFT bin values of desired frequency response')
title('Frequency response of band pass filter using ''inverse fft'' filter design technique')
set(gcf,'Position',[649   190   610   420]);

%%
% We'll now mimic the 'desired' filter shown in the video i.e. a band reject filter
H_desired = [1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1];
freq_bins_desired = [0:length(H_desired)-1]/length(H_desired)*2*pi;
b_ifft = ifft(H_desired);

% Now compute the continuous spectrum using any of the methods above. (Zero
% padding b_ifft prior to taking the fft is the most straighforward
% appraoch I think)
H_cont_desired = abs(fft([b_ifft zeros(1,1000)]));
w_cont_desired = [0:length(H_cont_desired)-1]/length(H_cont_desired)*2*pi;

figure
plot(w_cont_desired, H_cont_desired);
hold on
plot(freq_bins_desired, H_desired,'r.','MarkerSize', 12)

ylabel('Magnitude')
legend('Continuous Frequency Response','DFT bin values of desired frequency response')
title('Frequency response of band reject filter using ''inverse fft'' filter design technique')
set(gcf,'Position',[649   190   610   420]);

%%
% Now try using a 'practical' linear phase response rather than phase values of zero

H_desired = [1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1];
% An alternative way to create H_desired to ensure first half is a mirror
% image of the second half (for odd lengths only!)
% H_desired = [1 1 1 1 1 0 0 0 1 1 1 1];
% H_desired = [H_desired fliplr(H_desired(2:end))]

phase_diff = pi/length(H_desired)-pi;
phase_desired = [0:floor(length(H_desired)/2)]*phase_diff;
phase_desired = [phase_desired fliplr(phase_desired(2:end))*-1];

freq_bins_desired = [0:length(H_desired)-1]/length(H_desired)*2*pi;
h_ifft = ifft(H_desired.*exp(j*phase_desired));

b_ifft = h_ifft; %b coefficients of newly designed filter
%b_ifft = h_ifft.*hanning(length(h_ifft))'; %windowed b coefficients of newly designed filter

% Now compute the continuous spectrum using any of the methods above. (Zero
% padding b_ifft prior to taking the fft is the most straighforward
% approach - I think!)
H_cont_desired = abs(fft([b_ifft zeros(1,1000)]));
w_cont_desired = [0:length(H_cont_desired)-1]/length(H_cont_desired)*2*pi; %frequency axis

plot(w_cont_desired, H_cont_desired, 'Linewidth', 1);
hold on
plot(freq_bins_desired, H_desired,'r.','MarkerSize', 12)

ylabel('Magnitude')
legend('Continuous Frequency Response','DFT bin values of desired frequency response')
title('Frequency response of band reject filter using ''inverse fft'' filter design technique')
set(gcf,'Position',[649   190   610   420]);

%% Demonstration of filters being a sum of prototype bandpass filters
% H_acc is the accumulation of individual prototype bandpass filter

H_prot = zeros(size(H_desired)); %initialise
H_prot(1) = H_desired(1).*exp(j*phase_desired(1));
H_acc = fft(ifft(H_prot), length(H_cont_desired));

% loop through 11 prototype band pass filters
show_plots = 0;
if(show_plots)
figure
end
for k = 1: 11
prot_bins = [k+1 24-k] %pair of bins associated with positive and negative frequencies
H_prot = zeros(size(H_desired)); %initialise
H_prot(prot_bins) = H_desired(prot_bins).*exp(j*phase_desired(prot_bins));

H_cont_prot = fft(ifft(H_prot), length(H_cont_desired));
H_acc = H_acc + H_cont_prot; %accumulate each individual band pass prototype
if(show_plots)
plot(abs(H_acc))
hold on
plot(abs(H_cont_prot))
hold off
pause(1)
end
end

sum(abs(H_acc) - abs(H_cont_desired))



## Linear Phase Filters – why they are used

%% Linear phase filters - preserve shape of a filtered signal
% This is the code used during a youtube video presentation dealing with linear phase filters
%      Search for linear phase at http://youtube.com/ddorran
%
%
close all ; clear all; clc
fs = 100;
T = 1/fs; %sampling interval
N = 2000; %length of signal being synthesised
n = 0:N-1; %samples of the signal
t = n*T;

plot_range = [N/2-100:N/2+100];
%% synthesise a signal
x = cos(2*pi*10*t) + 0.5*cos(2*pi*20*t + 1.4);
subplot(2,1,1);
plot(t(plot_range),x(plot_range))
xlabel('Time (seconds)');
ylabel('Amplitude')
title('Synthesised Signals')
axis tight

ns = randn(1,length(x)+100)*2;
%filter the noise to synthesise band limited noise
[b a] = butter(5, [0.28 0.33],'bandpass');
ns_filtered = filter(b,a,ns);
x_ns = x +ns_filtered(end-length(x)+1:end);
hold on
noisy_x = plot(t(plot_range), x_ns(plot_range),'r');
legend('clean signal', 'noisy signal')

%% Plot frequency Content of Noisy Signal
subplot(2,1,2)
X_ns = fft(x_ns);
fax = [0:N-1]/(N/2); % normalised frequency axis
plot(fax(1:N/2), abs(X_ns(1:N/2))/(N/2)) ; %plot first half of spectrum
ylabel('Magnitude')
title('Magnitude Spectrum of Noisy Signal')

%% Filter out the noise using an IIR filter (non-linear phase)
[b_iir a_iir] = cheby1(10, 0.5, [0.27 0.34], 'stop');
y_iir = filter(b_iir,a_iir, x_ns);

[H_iir w] = freqz(b_iir,a_iir); %determine frequency response
subplot(2,1,2);
hold on
plot(w/pi, abs(H_iir),'r')
legend('|X(\omega)|','|H(\omega)|')

pause
Y_iir = fft(y_iir);
plot(fax(1:N/2), abs(Y_iir(1:N/2))/(N/2),'g') ; %plot first half of spectrum
legend('|X(\omega)|','|H(\omega)|','|Y(\omega)|')

pause
subplot(2,1,1)
non_linear_y = plot(t(plot_range),y_iir(plot_range),'g')
legend('clean signal', 'noisy signal','filtered signal')
pause
set(noisy_x,'visible', 'off')

%% Examine the magnitude and phase response of the IIR filter
figure(2)
subplot(2,1,1)
plot(w/pi,abs(H_iir))
ylabel('Magnitude')
title('Magnitude Response of filter')
subplot(2,1,2)
plot(w/pi,angle(H_iir))
ylabel('Phase Shift')
title('Phase Response of filter')

%% Now filter using an FIR filter (with linear phase)
b_fir = fir1(100,  [0.27 0.34],'stop');
a_fir = 1;
y_fir = filter(b_fir,a_fir, x_ns);

figure(1)
subplot(2,1,1)
plot(t(plot_range),y_fir(plot_range),'k')
legend('clean signal', 'noisy signal','filtered signal (non-linear)','filtered signal (linear)')

[H_fir, w ]= freqz(b_fir,a_fir);
subplot(2,1,2)
plot(w/pi, abs(H_fir),'k')
legend('|X(\omega)|','|H(\omega) Non-linear|','|Y(\omega)|','|H(\omega)| linear')

%% Compare the frequency responses of the two filter design approaches
figure(2)
subplot(2,1,1)
hold on
plot(w/pi,abs(H_fir),'g')
legend('non-linear filter','linear filter')
subplot(2,1,2)
hold on
plot(w/pi,angle(H_fir),'g')
legend('non-linear filter','linear filter')
pause

%% Why does linear phase preserve the shape??
close all
clear all; clc;
fs = 1000;
t = 0:1/fs:2;
x1 = cos(2*pi*3*t-pi/2);
x2 = cos(2*pi*5*t-(pi/2)/3*5);

pause
subplot(3,1,1)
plot(t,x1)
subplot(3,1,2)
plot(t,x2)
subplot(3,1,3)
plot(t,x1+x2,'g')
hold on



## pitch/period tracking using autocorrelation

%% Using Autocorrelation to track the local period of a signal
% This code is used as part of a youtube video demonstration
%
%
%       https://www.dropbox.com/s/3y25abf1xuqpizj/speech_demo.wav
%% speech analysis example

max_expected_period = round(1/50*fs);
min_expected_period = round(1/200*fs);
frame_len = 2*max_expected_period;

for k = 1 : length(ip)/frame_len -1;
range = (k-1)*frame_len + 1:k*frame_len;
frame = ip(range);

%show the input in blue and the selected frame in red
plot(ip);
set(gca, 'xtick',[],'position',[ 0.05  0.82   0.91  0.13])
hold on;
temp_sig = ones(size(ip))*NaN;
temp_sig(range) = frame;
plot(temp_sig,'r');
hold off

%use xcorr to determine the local period of the frame
[rxx lag] = xcorr(frame, frame);
subplot(3,1,3)
plot(lag, rxx,'r')
rxx(find(rxx < 0)) = 0; %set any negative correlation values to zero
center_peak_width = find(rxx(frame_len:end) == 0 ,1); %find first zero after center
%center of rxx is located at length(frame)+1
rxx(frame_len-center_peak_width : frame_len+center_peak_width  ) = min(rxx);
%     hold on
%     plot(lag, rxx,'g');
%     hold off
[max_val loc] = max(rxx);
period = abs(loc - length(frame)+1);

title(['Period estimate = ' num2str(period) 'samples (' num2str(fs/period) 'Hz)']);
set(gca, 'position', [ 0.05  0.07    0.91  0.25])

[max_val max_loc] = max(frame);
num_cycles_in_frame = ceil(frame_len/period);
test_start_positions = max_loc-(period*[-num_cycles_in_frame:num_cycles_in_frame]);
index = find(test_start_positions > 0,1, 'last');
start_position = test_start_positions(index);
colours = 'rg';

subplot(3,1,2)
plot(frame);

set(gca, 'position',[ 0.05 0.47 0.91 0.33])
pause
for g = 1 : num_cycles_in_frame
if(start_position+period*(g) <= frame_len && period > min_expected_period)
cycle_seg = ones(1, frame_len)*NaN;
cycle_seg(start_position+period*(g-1):start_position+period*(g))  =...
frame(start_position+period*(g-1):start_position+period*(g));
hold on

plot(cycle_seg,colours(mod(g, length(colours))+1)) %plot one of the available colors
hold off
end
end
pause
end

%% synthesise a periodic signal to use as a basic demo
fs = 500;
T = 1/fs;
N = 250; % desired length of signal
t = [0:N-1]*T; %time vector
f1 = 8; f2=f1*2;
x = sin(2*pi*f1*t-pi/2) + sin(2*pi*f2*t);
plot(t, x)
ylabel('Amplitude')
xlabel('Time (seconds)')
title('Synthesised Signal');

%% Determine the autocorrelation function
[rxx lags] = xcorr(x,x);
figure
plot(lags, rxx)
xlabel('Lag')
ylabel('Correlation Measure')
title('Auto-correlation Function')

%% Illustrate the auto correlation process
illustrate_xcorr(x,x)

%% Identify most prominent peaks
% Most prominent peak will be at the center of the correlation function
first_peak_loc = length(x) + 1;

% Lots of possible ways to identify second prominent peak. Am going to use a crude approach
% relying on some assumed prior knowledge of the signal. Am going to assume
% that the signal has a minimum possible period of .06 seconds = 30 samples;
min_period_in_samples = 30;
half_min = min_period_in_samples/2 ;

seq = rxx;
seq(first_peak_loc-half_min: first_peak_loc+half_min) = min(seq);
plot(rxx,'rx');
hold on
plot(seq)

[max_val second_peak_loc] = max(seq);
period_in_samples =  abs(second_peak_loc -first_peak_loc)
period = period_in_samples*T
fundamental_frequency = 1/period

%% Autocorrelation of a noisy signal
x2 = x + randn(1, length(x))*0.2;
plot(x2)
ylabel('Amplitude')
xlabel('Time (seconds)')
title('Noisy Synthesised Signal');

[rxx2 lags] = xcorr(x2,x2);
figure
plot(lags, rxx2)
xlabel('Lag')
ylabel('Correlation Measure')
title('Auto-correlation Function')

%% Autocorrelation technique can be problematic!
% Consider the following signal
f1 = 8; f2=f1*2;
x3 = sin(2*pi*f1*t) + 5*sin(2*pi*f2*t);
plot(t, x3)
ylabel('Amplitude')
xlabel('Time (seconds)')
title('Synthesised Signal');

[rxx3 lags] = xcorr(x3,x3,'unbiased');
figure
plot(lags, rxx3)
xlabel('Lag')
ylabel('Correlation Measure')
title('Auto-correlation Function')

seq = rxx3;
seq(first_peak_loc-half_min: first_peak_loc+half_min) = min(seq);
plot(seq)

[max_val second_peak_loc] = max(seq);
period_in_samples =  abs(second_peak_loc -first_peak_loc)



## illustrate_xcorr – code for cross correlation demos

September 24, 2014 1 comment
% This function illustrates the cross correlation process in action
%
% Usage:
%           fs = 1000;
%             T = 1/fs;
%             N = 500; % desired length of signal
%             t = [0:N-1]*T; %time vector
%             f1 = 8; f2=f1*2;
%             x = sin(2*pi*f1*t) + sin(2*pi*f2*t);
%
%           % To step though each sample use the following:
%           illustrate_xcorr(x,x)
%
%           % to step through using 50 steps use:
%           illustrate_xcorr(x,x, 50)
%
function illustrate_xcorr(x, y, varargin)
if(length(x) > length(y))
y(end+1:length(x)) = 0; %zero pad so the signals are the same length
else
x(end+1:length(y)) = 0; %zero pad so the signals are the same length
end

num_steps = 2*length(x)-1;
if(nargin ==3)
arg = varargin{1};
if(isnumeric(arg))
num_steps = ceil(abs(arg));
end
end
if(nargin > 3)
error('See help on this function to see how to use it properly')
end

[rxy lags] = xcorr(x,y); %cross correlate signals

disp('The signal being autocorrelated is shown in blue (two instances)')
disp('As you hit the space bar the lower plot will move into different lag positions')
disp('The correlation function shown in red is updated for each lag position');
disp('keep pressing the space bar to step through the illustration ...');

figure
plot_width = 0.3; plot_height = 0.25;

top_ax_h = subplot(3,1,1);
plot(x)
axis tight
set(top_ax_h, 'visible','off', 'units', 'normalized')
set(top_ax_h,'position', [0.5-plot_width/2 5/6-plot_height/2 plot_width plot_height])

mid_ax_h = subplot(3,1,2);
plot(y)
axis tight
set(mid_ax_h, 'visible','off', 'units', 'normalized')
set(mid_ax_h,'position', [0.5-plot_width/2-plot_width 5/6-3*plot_height/2-0.01 plot_width plot_height])

bottom_ax_h = subplot(3,1,3);
corr_h = plot(lags,rxy,'r');
axis tight
set(bottom_ax_h,'units', 'normalized','Ytick',[])
set(bottom_ax_h,'position', [0.5-plot_width*3/2 0.2-plot_height/2 plot_width*3 plot_height])
set(corr_h, 'Ydata', ones(1, length(rxy))*NaN); %clear the correlation funciton plot once its set up

normalised_shift_size = 2*plot_width/(num_steps-1);
corr_seg_len = length(rxy)/num_steps;
for k = 1 : num_steps
if(k > 1)
new_pos = get(mid_ax_h,'position') + [normalised_shift_size 0 0 0];
set(mid_ax_h,'position', new_pos);
end
set(corr_h, 'Ydata', [rxy(1:round(corr_seg_len*k)) ones(1,length(rxy)-round(corr_seg_len*k))*NaN])
pause
end



## Heartrate (BPM) Example Matlab Code

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

% A demonstration of cross correlation in action.
% Hit the space bar during the demo to execute
%
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

% illustration of Fourier Theory using plots
%
% usage : fourier_demonstration(1) % a well defined signal
%          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)
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

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

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

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

### Contents

% download data from 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 ## 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