Code Snippets Submitted by sparafucile17
Zero Crossing Counter
function count = zero_crossings(x)
% Count the number of zero-crossings from the supplied time-domain
% input vector. A simple method is applied here that can be easily
% ported to a real-time system that would minimize the number of
% if-else conditionals.
%
% Usage: COUNT = zero_crossings(X);
%
% X is the input time-domain signal (one dimensional)
% COUNT is the amount of zero-crossings in the input signal
%
% Author: sparafucile17 06/27/04
%
% initial value
count = 0;
% error checks
if(length(x) == 1)
error('ERROR: input signal must have more than one element');
end
if((size(x, 2) ~= 1) && (size(x, 1) ~= 1))
error('ERROR: Input must be one-dimensional');
end
% force signal to be a vector oriented in the same direction
x = x(:);
num_samples = length(x);
for i=2:num_samples
% Any time you multiply to adjacent values that have a sign difference
% the result will always be negative. When the signs are identical,
% the product will always be positive.
if((x(i) * x(i-1)) < 0)
count = count + 1;
end
endPhase changing Allpass IIR Filter
function [b,a] = allpass(order, Fc, Fs, Q)
% Returns allpass filter coefficients.
% Currently only 1st and 2nd orders are supported.
%
% Usage: [b,a] = ALLPASS(N, FC, FS, Q);
%
% N is the order of the allpass
% FC is the frequency a the 90deg phase shift point
% FS is the sampling rate
% Q is quality factor describing the slope of phase shift
%
% Author: sparafucile17 01/07/2004
if(order > 2)
error('Only 1st and 2nd orders are supported');
end
%Bilinear transform
g = tan(pi*(Fc/Fs));
if(order == 2)
d = 1/Q;
K = 1/(1 + d*g + g^2);
b0 = (1 - g*d + g^2) * K;
b1 = 2 * (g^2 - 1) * K;
b2 = 1;
a1 = b1;
a2 = b0;
else
b0 = (g-1)/(g+1);
b1 = 1;
b2 = 0;
a1 = b0;
a2 = 0;
end
b = [b0 b1 b2];
a = [1 a1 a2];Measuring Peak-to-Peak
function y = peak2peak(signal)
% Return the peak-to-peak amplitude of the supplied signal. This is the
% same as max(signal) minus min(signal).
%
% Usage: y = PEAK2PEAK(signal);
%
% SIGNAL is your one-dimensional input array
%
% Author: sparafucile17
% Input must have some length
if(length(signal) == 1)
error('ERROR: input signal must have more than one element');
end
% This function only supports one-dimensional arrays
if((size(signal, 2) ~= 1) && (size(signal, 1) ~= 1))
error('ERROR: Input must be one-dimensional');
end
% Find the peak and return it
min_sig = min(signal);
max_sig = max(signal);
% Answer
y = max_sig - min_sig;A-weighting filter
%Sampling Rate
Fs = 48000;
%Analog A-weighting filter according to IEC/CD 1672.
f1 = 20.598997;
f2 = 107.65265;
f3 = 737.86223;
f4 = 12194.217;
A1000 = 1.9997;
pi = 3.14159265358979;
NUM = [ (2*pi*f4)^2*(10^(A1000/20)) 0 0 0 0 ];
DEN = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]);
DEN = conv(conv(DEN,[1 2*pi*f3]),[1 2*pi*f2]);
%Bilinear transformation of analog design to get the digital filter.
[b,a] = bilinear(NUM,DEN,Fs);Percent Error
function [per_error] = percent_error(measured, actual)
% Creates the percent error between measured value and actual value
%
% Usage: percent_error(MEASURED, ACTUAL);
%
% Measured is the your result
% Actual is the value that your result should be
%
% Author: sparafucile17
per_error = abs(( (measured - actual) ./ actual ) * 100);Finding local minima
function [time,amp] = minima(signal, len)
%
% Finds the local minimum values of the given signal.
%
% Usage: [IND,AMP] = MINIMA(X, N);
%
% N is the number of points that need to be evaluated
% (Normally equal to LENGTH(X))
% X is the data
% IND contains the indices of X that contain the minima data
% AMP contains the minima amplitudes
%
% Author: sparafucile17 10/02/2003
%Initialize data
index = 1;
prev = signal(1);
local_time = 0;
local_amp = 0;
prev_slope = 1; %allow a maxima point at the second sample
%prevent length error
if(len > length(signal))
len = length(signal)
end
%Loop through data and find local minimums
for loop_i=2:len
cur = signal(loop_i);
slope = cur - prev;
if((cur < prev) & (prev_slope > 0)) %Positive slope and amplitude dips
local_time(index) = loop_i-1;
local_amp(index) = prev;
index = index+1;
end
prev = cur;
prev_slope = slope;
end
%After loop assign data to output variables
time = local_time;
amp = local_amp;Putting a Matlab figure to a specific screen location
function set_fig_position(h, top, left, height, width)
% Matlab has a wierd way of positioning figures so this function
% simplifies the poisitioning scheme in a more conventional way.
%
% Usage: SET_FIG_POISITION(h, top, left, height, width);
%
% H is the handle to the figure. This can be obtain in the
% following manner: H = figure(1);
% TOP is the "y" screen coordinate for the top of the figure
% LEFT is the "x" screen coordinate for the left side of the figure
% HEIGHT is how tall you want the figure
% WIDTH is how wide you want the figure
%
% Author: sparafucile17
% PC's active screen size
screen_size = get(0,'ScreenSize');
pc_width = screen_size(3);
pc_height = screen_size(4);
%Matlab also does not consider the height of the figure's toolbar...
%Or the width of the border... they only care about the contents!
toolbar_height = 77;
window_border = 5;
% The Format of Matlab is this:
%[left, bottom, width, height]
m_left = left + window_border;
m_bottom = pc_height - height - top - toolbar_height - 1;
m_height = height;
m_width = width - 1;
%Set the correct position of the figure
set(h, 'Position', [m_left, m_bottom, m_width, m_height]);
%If you want it to print to file correctly, this must also be set
% and you must use the "-r72" scaling to get the proper format
set(h, 'PaperUnits', 'points');
set(h, 'PaperSize', [width, height]);
set(h, 'PaperPosition', [0, 0, width, height]); %[ left, bottom, width, height]Batch WAV file processing
%In this example we are going to apply a low-pass filter to all WAV files,
%but it could be a multitude of different "processing" methods. This is
%only used to illustrate the batch processing example.
Fs = 44100; %Hz
Fc = 1000; %Hz
[b,a] = butter(2, Fc/(Fs/2), 'low');
%These are the input and output directories relative to this M-file
input_directory = 'test_input_database\';
output_direcory = 'test_output\';
%What extensions are you looking for? In our case, we want to batch
%process audio files that are store in the uncompressed *.WAV format
extension = '*.wav';
%Get the files
files=dir([input_directory '*.wav']);
N_files=numel(files);
%Loop through one file at a time
for i=1:N_files
%This is a stereo file and for this example,
%I only care about channel number 1
x = wavread([input_directory files(i).name ]);
x = x(:,1);
%Process the data by applying our filter
y = filter(b,a,x);
%Output the data as a unique filename based on input
wavwrite(y, Fs, [output_directory 'processed' files(i).name]);
disp(['Processed File: ' input_directory files(i).name 'as: ' ...
output_directory 'processed' files(i).name]);
endC-weighting filter
%Sampling Rate
Fs = 48000;
%Analog C-weighting filter according to IEC/CD 1672.
1672.
f1 = 20.598997;
f4 = 12194.217;
C1000 = 0.0619;
pi = 3.14159265358979;
NUM = [ (2*pi*f4)^2*(10^(C1000/20)) 0 0 ];
DEN = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]);
%Bilinear transformation of analog design to get the digital [b,a] = bilinear(NUM,DEN,Fs);Interleave Zeros
function out = interleave_zeros(a)
% Interleave vector with zeros.
% Only 1D matrices are supported.
%
% Usage: b = INTERLEAVE_ZEROS(a);
%
% A is the vector that will have zeros inserted
% B will contain the vector set A and the interleaved
% zeros (twice the length)
%
% Author: sparafucile17 2/13/04
my_zeros = zeros(1, length(a));
out = matintrlv([a my_zeros],2,length(a));Saving a Matlab figure to file using command-line
function save_fig(h, name, format)
% Matlab has a wierd way of saving figures and keeping the plot looking
% exactly the way it does on your screen. Thi function simplifies the
% entire process in a way that is barely documented.
%
% Usage: SAVE_FIG(h, name, format);
%
% H is the handle to the figure. This can be obtain in the
% following manner: H = figure(1);
% NAME is the filename without extension
% FORMAT is the graphic format. Options are:
%
% 'bmpmono' BMP monochrome BMP
% 'bmp16m' BMP 24-bit BMP
% 'bmp256' BMP 8-bit (256-color)
% 'bmp' BMP 24-bit
% 'meta' EMF
% 'eps' EPS black and white
% 'epsc' EPS color
% 'eps2' EPS Level 2 black and white
% 'epsc2' EPS Level 2 color
% 'fig' FIG Matlab figure file
% 'hdf' HDF 24-bit
% 'ill' ILL (Adobe Illustrator)
% 'jpeg' JPEG 24-bit
% 'pbm' PBM (plain format) 1-bit
% 'pbmraw' PBM (raw format) 1-bit
% 'pcxmono' PCX 1-bit
% 'pcx24b' PCX 24-bit color PCX file format
% 'pcx256' PCX 8-bit newer color (256-color)
% 'pcx16' PCX 4-bit older color (16-color)
% 'pdf' PDF Color PDF file format
% 'pgm' PGM Portable Graymap (plain format)
% 'pgmraw' PGM Portable Graymap (raw format)
% 'png' PNG 24-bit
% 'ppm' PPM Portable Pixmap (plain format)
% 'ppmraw' PPM Portable Pixmap (raw format)
% 'svg' SVG Scalable Vector Graphics
% 'tiff' TIFF 24-bit
%
% Author: sparafucile17
if(strcmp(format,'fig'))
saveas(h, name, 'fig');
else
options.Format = format;
hgexport(h, name, options);
endPutting a Matlab figure to a specific screen location
function set_fig_position(h, top, left, height, width)
% Matlab has a wierd way of positioning figures so this function
% simplifies the poisitioning scheme in a more conventional way.
%
% Usage: SET_FIG_POISITION(h, top, left, height, width);
%
% H is the handle to the figure. This can be obtain in the
% following manner: H = figure(1);
% TOP is the "y" screen coordinate for the top of the figure
% LEFT is the "x" screen coordinate for the left side of the figure
% HEIGHT is how tall you want the figure
% WIDTH is how wide you want the figure
%
% Author: sparafucile17
% PC's active screen size
screen_size = get(0,'ScreenSize');
pc_width = screen_size(3);
pc_height = screen_size(4);
%Matlab also does not consider the height of the figure's toolbar...
%Or the width of the border... they only care about the contents!
toolbar_height = 77;
window_border = 5;
% The Format of Matlab is this:
%[left, bottom, width, height]
m_left = left + window_border;
m_bottom = pc_height - height - top - toolbar_height - 1;
m_height = height;
m_width = width - 1;
%Set the correct position of the figure
set(h, 'Position', [m_left, m_bottom, m_width, m_height]);
%If you want it to print to file correctly, this must also be set
% and you must use the "-r72" scaling to get the proper format
set(h, 'PaperUnits', 'points');
set(h, 'PaperSize', [width, height]);
set(h, 'PaperPosition', [0, 0, width, height]); %[ left, bottom, width, height]Least Squares fit to a linear equation
function [m, b] = ls_linear(x_observed, y_observed)
%
% Perform Least Squares curve-fitting techniques to solve the coefficients
% of the linear equation: y = m*x + b. Input and Output equation
% observations must be fed into the algorithm and a best fit equation will
% be calculated.
%
% Usage: [M,B] = ls_linear(observed_input, observed_output);
%
% observed_input is the x-vector observed data
% observed_output is the y-vector observed data
%
% Author: sparafucile17 06/04/03
%
if(length(x_observed) ~= length(y_observed))
error('ERROR: Both X and Y vectors must be the same size!');
end
%Calculate vector length once
vector_length = length(x_observed);
%
% Theta = [ x1 1 ]
% [ x2 1 ]
% [ x3 1 ]
% [ x4 1 ]
% [ ... ... ]
%
Theta = ones(vector_length, 2);
Theta(1:end, 1) = x_observed;
%
% Theta_0 = [ y1 ]
% [ y2 ]
% [ y3 ]
% [ y4 ]
% [ ... ]
%
Theta_0 = y_observed;
%
% Theta_LS = [ M ]
% [ B ]
%
Theta_LS = ((Theta' * Theta)^-1 ) * Theta' * Theta_0;
%Return the answer
m = Theta_LS(1);
b = Theta_LS(2);Measuring Peak-to-Peak
function y = peak2peak(signal)
% Return the peak-to-peak amplitude of the supplied signal. This is the
% same as max(signal) minus min(signal).
%
% Usage: y = PEAK2PEAK(signal);
%
% SIGNAL is your one-dimensional input array
%
% Author: sparafucile17
% Input must have some length
if(length(signal) == 1)
error('ERROR: input signal must have more than one element');
end
% This function only supports one-dimensional arrays
if((size(signal, 2) ~= 1) && (size(signal, 1) ~= 1))
error('ERROR: Input must be one-dimensional');
end
% Find the peak and return it
min_sig = min(signal);
max_sig = max(signal);
% Answer
y = max_sig - min_sig;Percent Error
function [per_error] = percent_error(measured, actual)
% Creates the percent error between measured value and actual value
%
% Usage: percent_error(MEASURED, ACTUAL);
%
% Measured is the your result
% Actual is the value that your result should be
%
% Author: sparafucile17
per_error = abs(( (measured - actual) ./ actual ) * 100);Finding local minima
function [time,amp] = minima(signal, len)
%
% Finds the local minimum values of the given signal.
%
% Usage: [IND,AMP] = MINIMA(X, N);
%
% N is the number of points that need to be evaluated
% (Normally equal to LENGTH(X))
% X is the data
% IND contains the indices of X that contain the minima data
% AMP contains the minima amplitudes
%
% Author: sparafucile17 10/02/2003
%Initialize data
index = 1;
prev = signal(1);
local_time = 0;
local_amp = 0;
prev_slope = 1; %allow a maxima point at the second sample
%prevent length error
if(len > length(signal))
len = length(signal)
end
%Loop through data and find local minimums
for loop_i=2:len
cur = signal(loop_i);
slope = cur - prev;
if((cur < prev) & (prev_slope > 0)) %Positive slope and amplitude dips
local_time(index) = loop_i-1;
local_amp(index) = prev;
index = index+1;
end
prev = cur;
prev_slope = slope;
end
%After loop assign data to output variables
time = local_time;
amp = local_amp;Finding local maxima
function [time,amp] = maxima(signal, thresh, len)
% Finds the local maximum values of the given signal.
%
% Usage: [IND,AMP] = MAXIMA(X, THRESH, N);
%
% X is the data
% THRESH is the threshold that signal must be greater
% than in order to be counted. (Default = 0.0)
% N is the number of points that need to be evaluated
% (Defaults = LENGTH(X))
% IND contains the indices of X that contain the maxima data
% AMP contains the maxima amplitudes
%
% Author: sparafucile17 10/02/2003
if(exist('len') == 0)
len = length(signal);
end
if(exist('thresh') == 0)
thresh = 0.0;
end
%Initialize data
index = 1;
prev = signal(1);
local_time = 0;
local_amp = 0;
prev_slope = 1; %allow a maxima point at the second sample
%prevent length error
if(len > length(signal))
len = length(signal)
end
%Loop through data and find local maximums
for loop_i=2:len
cur = signal(loop_i);
slope = cur - prev;
if((cur < prev) & (prev_slope > 0) & (cur > thresh)) %Positive slope and amplitude dips
local_time(index) = loop_i-1;
local_amp(index) = prev;
index = index+1;
end
prev = cur;
prev_slope = slope;
end
%After loop assign data to output variables
time = local_time;
amp = local_amp;Pink Noise Generator
function out = create_pink_noise(Fs, Sec, Amp)
% Creates a pink noise signal and saves it as a wav file
%
% Usage: create_noise(Fs, Sec, Amp);
%
% Fs is the desired sampling rate
% Sec is the duration of the signal in seconds
% Amp is the amplitude in dB of the signal (0dB to -144dB)
%
% Author: sparafucile17 06/14/02
%error trapping
if((Amp > 0) || (Amp < -144))
error('Amplitude is not within the range of 0dB to -144dB');
end
%Create Whitenoise
white_noise = randn((Fs*Sec)+1,1);
%Apply weighted sum of first order filters to approximate a -10dB/decade
%filter. This is Paul Kellet's "refined" method (a.k.a instrumentation
%grade) It is accurate to within +/-0.05dB above 9.2Hz
b=zeros(7,1);
for i=1:((Fs*Sec)+1)
b(1) = 0.99886 * b(1) + white_noise(i) * 0.0555179;
b(2) = 0.99332 * b(2) + white_noise(i) * 0.0750759;
b(3) = 0.96900 * b(3) + white_noise(i) * 0.1538520;
b(4) = 0.86650 * b(4) + white_noise(i) * 0.3104856;
b(5) = 0.55000 * b(5) + white_noise(i) * 0.5329522;
b(6) = -0.7616 * b(6) - white_noise(i) * 0.0168980;
pink_noise(i) = b(1) + b(2) + b(3) + b(4) + b(5) + b(6) + b(7) + white_noise(i) * 0.5362;
b(7) = white_noise(i) * 0.115926;
end
%Normalize to +/- 1
if(abs(min(pink_noise)) > max(pink_noise))
pink_noise = pink_noise / abs(min(pink_noise));
else
pink_noise = pink_noise / max(pink_noise);
end
%Normalize to prevent positive saturation (We can't represent +1.0)
pink_noise = pink_noise /abs(((2^31)-1)/(2^31));
%Scale signal to match desired level
pink_noise = pink_noise * 10^(Amp/20);
%Output noise signal
out = pink_noise(1:end-1);C-weighting filter
%Sampling Rate
Fs = 48000;
%Analog C-weighting filter according to IEC/CD 1672.
1672.
f1 = 20.598997;
f4 = 12194.217;
C1000 = 0.0619;
pi = 3.14159265358979;
NUM = [ (2*pi*f4)^2*(10^(C1000/20)) 0 0 ];
DEN = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]);
%Bilinear transformation of analog design to get the digital [b,a] = bilinear(NUM,DEN,Fs);A-weighting filter
%Sampling Rate
Fs = 48000;
%Analog A-weighting filter according to IEC/CD 1672.
f1 = 20.598997;
f2 = 107.65265;
f3 = 737.86223;
f4 = 12194.217;
A1000 = 1.9997;
pi = 3.14159265358979;
NUM = [ (2*pi*f4)^2*(10^(A1000/20)) 0 0 0 0 ];
DEN = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]);
DEN = conv(conv(DEN,[1 2*pi*f3]),[1 2*pi*f2]);
%Bilinear transformation of analog design to get the digital filter.
[b,a] = bilinear(NUM,DEN,Fs);
