Archive
Archive for the ‘youtube demo code’ Category
Negative Frequency GUI Code
March 16, 2022
Leave a comment
function negfgui(varargin)
% An interactive GUI which attempts to help students visualise negative frequencies
% by David Dorran, TU Dublin, March 2022
% reuses some code to track mouse movement from the zpgui function by Tom Krauss from Perdue University
% Youtube video explanation provided at https://youtu.be/Rxc_ypVdruw
global zero_locations diff_equ_text op_text_samps op_freq_pts op_points draw_plt disp_plt op_disp_plt freq_plt xyz_pts op_xyz_pts xy_pts ptr negfqui_fh freq_pts points text_samps b a y gain pole_vals zero_vals freq_resp_pts
buf_len = 50;
tail_len = 10;
max_points = 100000;
fft_len = 1024;
win_func = hamming(buf_len);
if nargin == 0
action = 'init';
else
action = varargin{1};
end
if(strcmp(class(action), 'matlab.ui.control.UIControl'))
zero_vals = str2num(zero_locations.String);
negfgui('update_freq_response');
end
switch action
case 'init'
zero_vals = [ -j j];
pole_vals = [0 0];
if(strcmp(class(negfqui_fh),'matlab.ui.Figure'))
if(~isvalid(negfqui_fh))
negfqui_fh = figure;
end
if(nargin ==0)
close(negfqui_fh);
negfqui_fh = figure;
end
else
negfqui_fh = figure;
end
set(negfqui_fh,'position',[ 1 41 1366 651]);
ptr = 1;
points = ones(1,max_points)*NaN;
op_points = ones(1,max_points)*NaN;
freq_plt = subplot(2,3,1);
fax = (([0:fft_len-1])/fft_len*2*pi)-pi;
yyaxis left
freq_pts = plot(fax,points(1:fft_len));
hold on
op_freq_pts = plot(fax,points(1:fft_len),'r');
xlim([-pi pi])
ylim([0 buf_len/2])
ylabel('Magnitude Content')
yyaxis right
freq_resp_pts = plot(fax,points(1:fft_len));
xlim([-pi pi])
ylim([0 1])
title('Double sided-spectrum');
xlabel('frequency (radians per sample)')
ylabel('Magnitude Response')
draw_plt = subplot(2,3,2);
xy_pts = plot([]);
for m =1:tail_len
hold on
xy_pts(m) = plot(points(m)+points(m)*j,'r.','MarkerSize',tail_len-m+1);
end
ylabel('Y')
xlabel('X')
xlim([-1 1]);
ylim([-1 1]);
disp_plt = subplot(2,3,3);
xyz_pts = plot3(real(points),imag(points),[0:length(points)-1]);
hold on
surf([0 0 ],[-1.1 1.1],[-1 length(points); -1 length(points)],'FaceAlpha',0.2)
surf([-1.1 1.1], [ 0 0 ],[ length(points) length(points); -1 -1],'FaceAlpha',0.2)
xlabel('X')
ylabel('Y')
zlabel('Samples')
xlim([-1 1]);
ylim([-1 1]);
zlim([0 buf_len])
set(gca,'CameraUpVector', [0 1 0 ])
grid on
set(draw_plt,'buttondownfcn', 'negfgui(''draw'')')
op_disp_plt = subplot(2,3,6);
op_xyz_pts = plot3(real(op_points),imag(op_points),[0:length(op_points)-1]);
hold on
surf([0 0 ],[-1.1 1.1],[-1 length(op_points); -1 length(op_points)],'FaceAlpha',0.2)
surf([-1.1 1.1], [ 0 0 ],[ length(op_points) length(op_points); -1 -1],'FaceAlpha',0.2)
xlabel('X')
ylabel('Y')
zlabel('Samples')
xlim([-1 1]);
ylim([-1 1]);
zlim([0 buf_len])
set(gca,'CameraUpVector', [0 1 0 ])
grid on
text_samps = annotation('textbox','position',[.0500 0.48 0.85 0.05],'LineStyle','none','FontSize',10);
op_text_samps = annotation('textbox','position',[.0500 0 0.85 0.05],'LineStyle','none','FontSize',10);
diff_equ_text = annotation('textbox',...
[0.34 0.337941628264209 0.285237188872621 0.0844854070660522],...
'String',[''],...
'FontSize',14,...
'FitBoxToText','off',...
'EdgeColor','none');
annotation('textbox',...
[0.34 0.2 0.285237188872621 0.0844854070660522],...
'String',['Zero Locations:'],...
'FontSize',14,...
'FitBoxToText','off',...
'EdgeColor','none');
zero_locations = uicontrol('style','edit','units','normalized','FontSize',14,'position', [0.34 0.17 0.2 0.05]);
update_button = uicontrol('style','pushbutton', 'string','update locations','units','normalized','position',[0.34 0.12 0.2 0.05],'callback',{@negfgui});
negfgui('update_freq_response')
case 'update_freq_response'
b = poly(zero_vals);
a = poly(pole_vals);
H = freqz(b,a,fft_len,'whole');
gain = max(abs(H));
freq_resp_pts.YData= fftshift(abs(H)/gain); %ensure max gain is 1
b = b/gain;
subplot(2,3,4)
zplane(b,a);
title('Pole-Zero Plot')
try
str = create_difference_equation(b,1,'tex');
diff_equ_text.String = {['Difference Equation:'], str};
catch
diff_equ_text.String = {['b coefficients: ' num2str(b)], ['a coefficients:' num2str(a)]};
end
zero_locations.String = num2str(zero_vals);
case 'draw'
set(gcf,'userdata','')
set(gcf,'windowbuttonmotionfcn','set(gcf,''userdata'',''motion'')')
set(gcf,'windowbuttonupfcn','set(gcf,''userdata'',''up'')')
done = 0;
while ~done
waitfor(gcf,'userdata')
switch get(gcf,'userdata')
case 'motion'
pt = get(draw_plt,'currentpoint');
if(abs(pt(1,1))>1)
pt(1,1) = 1*sign(pt(1,1));
end
if(abs(pt(1,2))>1)
pt(1,2) = 1*sign(pt(1,2));
end
points(ptr) = pt(1,1)+pt(1,2)*j;
xyz_pts.XData(ptr) = real(points(ptr));
xyz_pts.YData(ptr) = imag(points(ptr));
%filter the signal (it's not perfect!)
try
op_points(ptr) = 0;
for k= 1:length(b)
op_points(ptr) = op_points(ptr)+ b(k)*points(ptr-k+1);
end
catch
end
op_xyz_pts.XData(ptr) = real(op_points(ptr));
op_xyz_pts.YData(ptr) = imag(op_points(ptr));
if(ptr>buf_len)
disp_plt.ZLim = [ptr-buf_len ptr];
op_disp_plt.ZLim = [ptr-buf_len ptr];
freq_pts.YData = abs(fftshift(fft(points(ptr-buf_len+1:ptr).*win_func', fft_len)));
op_freq_pts.YData = abs(fftshift(fft(op_points(ptr-buf_len+1:ptr).*win_func', fft_len)));
end
for m =1:tail_len
try
xy_pts(m).XData(1) = real(points(ptr-m+1));
xy_pts(m).YData(1) = imag(points(ptr-m+1));
text_samps.String = num2str(fliplr(round(points(ptr-m+1:ptr)*100)/100));
op_text_samps.String = num2str(fliplr(round(op_points(ptr-m+1:ptr)*100)/100));
catch
text_samps.String = num2str(fliplr(round(points(1:ptr)/100)*100));
end
end
ptr=ptr+1;
if(ptr> max_points)
negfgui('init')
end
case 'up'
done = 1;
end
set(gcf,'userdata','')
end
set(gcf,'windowbuttonmotionfcn','')
set(gcf,'windowbuttonupfcn','')
end
Categories: matlab code, Uncategorized, youtube demo code
FFT based FIR filter design
July 19, 2019
Leave a comment
%% 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)
xlabel('Frequency (radians per sample)')
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)
xlabel('Frequency (radians per sample)')
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)
xlabel('Frequency (radians per sample)')
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)
xlabel('Frequency (radians per sample)')
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
%start with the DC prototype filter
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))
Categories: matlab code, youtube demo code
Linear Phase Filters – why they are used
October 1, 2014
Leave a comment
%% 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
%
% Code available from https://dadorran.wordpress.com
%
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
% Add some noise
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);
%add noise to clean signal
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
xlabel('frequency ( x \pi rads/sample)')
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))
xlabel('frequency ( x \pi rads/sample)')
ylabel('Magnitude')
title('Magnitude Response of filter')
subplot(2,1,2)
plot(w/pi,angle(H_iir))
xlabel('frequency ( x \pi rads/sample)')
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
Categories: matlab code, youtube demo code
pitch/period tracking using autocorrelation
September 24, 2014
Leave a comment
%% Using Autocorrelation to track the local period of a signal
% This code is used as part of a youtube video demonstration
% See http://youtube.com/ddorran
%
% Code available at https://dadorran.wordpress.com
%
% The following wav file can be downloaded from
% https://www.dropbox.com/s/3y25abf1xuqpizj/speech_demo.wav
%% speech analysis example
[ip fs] = wavread('speech_demo.wav');
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
%function available from https://dadorran.wordpress.com
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)
Categories: matlab code, youtube demo code
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
Categories: matlab code, youtube demo code
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) > sig(k-1) & sig(k) > sig(k+1) & sig(k) > 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;
Categories: matlab code, Uncategorized, youtube demo code
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');
Categories: matlab code, Uncategorized, youtube demo code
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
Categories: matlab code, Uncategorized, youtube demo code
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)
Categories: matlab code, Uncategorized, youtube demo code
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.
Contents
- Load Example Data
- Quick view of double-sided (two-sided) magnitude spectrum
- Double-sided magnitude spectrum with frequency axis (in bins)
- Single-sided magnitude spectrum with frequency axis in bins
- Single-sided magnitude spectrum with frequency axis in Hertz
- Single-sided magnitude spectrum with frequency axis normalised
- Single-sided magnitude spectrum – frequency in rads per sample
- Double-sided magnitude spectrum showing negative frequencies
- Single-sided magnitiude spectrum in decibels and Hertz
- Single-sided power spectrum in decibels and Hertz
- Single-sided power spectrum in dB and frequency on a log scale
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
Categories: matlab code, Uncategorized, youtube demo code
dadorran 