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
end
```

## Phase 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];
```

## Shelving Filter Design

```
function [b, a] = shelving(G, fc, fs, Q, type)
%
% Derive coefficients for a shelving filter with a given amplitude and
% cutoff frequency. All coefficients are calculated as described in
% Zolzer's DAFX book (p. 50 -55).
%
% Usage: [B,A] = shelving(G, Fc, Fs, Q, type);
%
% G is the logrithmic gain (in dB)
% FC is the center frequency
% Fs is the sampling rate
% Q adjusts the slope be replacing the sqrt(2) term
% type is a character string defining filter type
% Choices are: 'Base_Shelf' or 'Treble_Shelf'
%
% Author: sparafucile17 08/22/05
%
%Error Check
if((strcmp(type,'Base_Shelf') ~= 1) && (strcmp(type,'Treble_Shelf') ~= 1))
error(['Unsupported Filter Type: ' type]);
end
K = tan((pi * fc)/fs);
V0 = 10^(G/20);
root2 = 1/Q; %sqrt(2)
%Invert gain if a cut
if(V0 < 1)
V0 = 1/V0;
end
%%%%%%%%%%%%%%%%%%%%
% BASE BOOST
%%%%%%%%%%%%%%%%%%%%
if(( G > 0 ) & (strcmp(type,'Base_Shelf')))
b0 = (1 + sqrt(V0)*root2*K + V0*K^2) / (1 + root2*K + K^2);
b1 = (2 * (V0*K^2 - 1) ) / (1 + root2*K + K^2);
b2 = (1 - sqrt(V0)*root2*K + V0*K^2) / (1 + root2*K + K^2);
a1 = (2 * (K^2 - 1) ) / (1 + root2*K + K^2);
a2 = (1 - root2*K + K^2) / (1 + root2*K + K^2);
%%%%%%%%%%%%%%%%%%%%
% BASE CUT
%%%%%%%%%%%%%%%%%%%%
elseif (( G < 0 ) & (strcmp(type,'Base_Shelf')))
b0 = (1 + root2*K + K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);
b1 = (2 * (K^2 - 1) ) / (1 + root2*sqrt(V0)*K + V0*K^2);
b2 = (1 - root2*K + K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);
a1 = (2 * (V0*K^2 - 1) ) / (1 + root2*sqrt(V0)*K + V0*K^2);
a2 = (1 - root2*sqrt(V0)*K + V0*K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);
%%%%%%%%%%%%%%%%%%%%
% TREBLE BOOST
%%%%%%%%%%%%%%%%%%%%
elseif (( G > 0 ) & (strcmp(type,'Treble_Shelf')))
b0 = (V0 + root2*sqrt(V0)*K + K^2) / (1 + root2*K + K^2);
b1 = (2 * (K^2 - V0) ) / (1 + root2*K + K^2);
b2 = (V0 - root2*sqrt(V0)*K + K^2) / (1 + root2*K + K^2);
a1 = (2 * (K^2 - 1) ) / (1 + root2*K + K^2);
a2 = (1 - root2*K + K^2) / (1 + root2*K + K^2);
%%%%%%%%%%%%%%%%%%%%
% TREBLE CUT
%%%%%%%%%%%%%%%%%%%%
elseif (( G < 0 ) & (strcmp(type,'Treble_Shelf')))
b0 = (1 + root2*K + K^2) / (V0 + root2*sqrt(V0)*K + K^2);
b1 = (2 * (K^2 - 1) ) / (V0 + root2*sqrt(V0)*K + K^2);
b2 = (1 - root2*K + K^2) / (V0 + root2*sqrt(V0)*K + K^2);
a1 = (2 * ((K^2)/V0 - 1) ) / (1 + root2/sqrt(V0)*K + (K^2)/V0);
a2 = (1 - root2/sqrt(V0)*K + (K^2)/V0) / (1 + root2/sqrt(V0)*K + (K^2)/V0);
%%%%%%%%%%%%%%%%%%%%
% All-Pass
%%%%%%%%%%%%%%%%%%%%
else
b0 = V0;
b1 = 0;
b2 = 0;
a1 = 0;
a2 = 0;
end
%return values
a = [ 1, a1, a2];
b = [ b0, b1, b2];
```

## 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);
```

## First order RC Filter in the Digital Domain

```
function [b,a] = rc_filter(R, C, Fs, filter_type)
% Returns equivalent IIR coefficients for an analog RC filter
%
% Usage: [B,A] = RC_FILTER(r, c, fs, type);
%
% R is the resistance value (in ohms)
% C is the capacitance value (in farrads)
% FS is the digital sample rate (in Hz)
% type is a character string defining filter type
% Choices are: 'high' or 'low'
%
% Highpass filters have the analog structure:
%
%
% | |
% Vi o--------| |----------+---------o Vo
% | | |
% C |
% ---
% | | R
% | |
% ---
% |
% |
% o---------------------+---------o
% GND
%
%
% Lowpass filters have the analog structure:
%
%
% |-----|
% Vi o--------| |------+---------o Vo
% |-----| |
% R |
% ----- C
% -----
% |
% |
% o---------------------+---------o
% GND
%
% This function uses a pre-calculated equation for both of these circuits
% that only requires the resistance and capacitance value to get a true
% digital filter equivalent to a basic analog filter.
%
% The math behind these equations is based off the basic bilinear transform
% technique that can be found in many DSP textbooks. The reference paper
% for this function was "Conversion of Analog to Digital Transfer
% Functions" by C. Sidney Burrus, page 6.
%
% This is also the equivalent of a 1st order butterworth with a cuttoff
% frequency of Fc = 1/(2*pi*R*C);
%
% Author: sparafucile17 07/01/02
%
% Verify that cutoff of the analog filter is below Nyquist
if( (1/(2*pi*R*C)) > (Fs/2) )
error('This analog filter cannot be realized with this sample rate');
end
% Default to allpass if invalid type is selected
b = [ 1 0 ];
a = [ 1 0 ];
% Constants
RC = R * C;
T = 1 / Fs;
% Analog Cutoff Fc
w = 1 / (RC);
% Prewarped coefficient for Bilinear transform
A = 1 / (tan((w*T) / 2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following equations were derived from
%
% s
% T(s) = -------
% s + 1
%
%
% using Bilinear transform of
%
% 1 ( 1 - z^-1 )
% s --> ----------- * ------------
% tan(w*T/2) ( 1 + z^-1 )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(strcmp(filter_type,'high'))
b(1) = (A) / (1 + A);
b(2) = -b(1);
a(2) = (1 - A) / (1 + A);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following equations were derived from
%
% 1
% T(s) = -------
% s + 1
%
%
% using Bilinear transform of
%
% 1 ( 1 - z^-1 )
% s --> ----------- * ------------
% tan(w*T/2) ( 1 + z^-1 )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(strcmp(filter_type,'low'))
b(1) = (1) / (1 + A);
b(2) = b(1);
a(2) = (1 - A) / (1 + A);
end
```

## 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]);
end
```

## 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);
end
```

## 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]
```

## 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);
```