Archive

Archive for the ‘youtube demo code’ Category

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

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

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)


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

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.

Contents

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

plotting_frequency_spectrum_demo_01

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

plotting_frequency_spectrum_demo_02

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

plotting_frequency_spectrum_demo_03

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

plotting_frequency_spectrum_demo_04

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

plotting_frequency_spectrum_demo_05

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

plotting_frequency_spectrum_demo_06

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

plotting_frequency_spectrum_demo_07

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

plotting_frequency_spectrum_demo_08

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

plotting_frequency_spectrum_demo_09

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_frequency_spectrum_demo_10

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