Archive

Archive for February, 2014

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

% 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

Categories: matlab code, Uncategorized, youtube demo code

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
Categories: matlab code, Uncategorized, youtube demo code

create_tf_equation

February 3, 2014 Leave a comment

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

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

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

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


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

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