## Measuring Peak-to-Peak

September 30, 2011 Coded in Matlab
``````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);

y = max_sig - min_sig;``````

## Percent Error

September 30, 2011 Coded in Matlab
``````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

September 26, 20111 comment Coded in Matlab
``````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

September 26, 20112 comments Coded in Matlab
``````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

September 26, 2011 Coded in Matlab
``````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

September 25, 2011 Coded in Matlab
``````%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

September 25, 20111 comment Coded in Matlab
``````%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);``````

## Auto Wah audio effect (LFO control)

September 24, 2011 Coded in C
``````/*************** AutoWah.c ********************************/

#include "bp_iir.h"
#include "AutoWah.h"

static short center_freq;
static short samp_freq;
static short counter;
static short counter_limit;
static short control;
static short max_freq;
static short min_freq;
static short f_step;
static struct bp_filter H;

/*
This is the auto wah effect initialization function.
This initializes the band pass filter and the effect control variables
*/
void AutoWah_init(short effect_rate,short sampling,short maxf,short minf,short Q,double gainfactor,short freq_step) {
double C;

/*Process variables*/
center_freq = 0;
samp_freq = sampling;
counter = effect_rate;
control = 0;

/*User Parametters*/
counter_limit = effect_rate;

/*Convert frequencies to index ranges*/
min_freq = 0;
max_freq = (maxf - minf)/freq_step;

bp_iir_init(sampling,gainfactor,Q,freq_step,minf);
f_step = freq_step;
}

/*
This function generates the current output value
Note that if the input and output signal are integer
unsigned types, we need to add a half scale offset
*/
double AutoWah_process(int xin) {
double yout;

yout = bp_iir_filter(xin,&H);
#ifdef INPUT_UNSIGNED
yout += 32767;
#endif

return yout;
}

/*
This function will emulate a LFO that will vary according
to the effect_rate parameter set in the AutoWah_init function.
*/
void AutoWah_sweep(void) {
double yout;

if (!--counter) {
if (!control) {
bp_iir_setup(&H,(center_freq+=f_step));
if (center_freq > max_freq) {
control = 1;
}
}
else if (control) {
bp_iir_setup(&H,(center_freq-=f_step));
if (center_freq == min_freq) {
control = 0;
}
}

counter = counter_limit;
}
}

/*************** AutoWah.h ****************************/

#ifndef __AUTOWAH_H__
#define __AUTOWAH_H__

extern void AutoWah_init(short effect_rate,short sampling,short maxf,short minf,short Q,double gainfactor,short freq_step);
extern double AutoWah_process(int xin);
extern void AutoWah_sweep(void);

#endif

/************** Usage Example **************************/

#include "AutoWah.h"

void main(void) {
short xin;
short yout;
AutoWah_init(2000,  /*Effect rate 2000*/
16000, /*Sampling Frequency*/
1000,  /*Maximum frequency*/
500,   /*Minimum frequency*/
4,     /*Q*/
0.707, /*Gain factor*/
10     /*Frequency increment*/
);

while(1) {
if (new_sample_flag()) {
/*Apply AutoWah_process function to the sample*/
yout =AutoWah_process(xin);

/*Send the output value to your DAC or codec output*/
write_output(yout);
/*Makes the LFO vary*/
AutoWah_sweep();
}
}
}``````

## GSM (GMSK) power spectrum equation

September 18, 2011 Coded in Matlab
``````% ************************************************************
% Spectrum model for GSM signal
% Markus Nentwig, 2011
% based on 3GPP TS 45.004 section 2 "Modulation format for GMSK",
% assuming no additional filtering and an infinite-length
% symbol stream (no slot structure)
% ************************************************************

% ************************************************************
% Parameters
% The "baseline" serves as a gentle reminder that the model
% is only valid near the center frequency.
% the particular transmitter used.
% If not available, assume that spectral emission mask requirements
% are met.
% ************************************************************
BW = 285e3;
BW2 = 186e3;
baseline_dB = -76;
% baseline_dB = -999 % disable the constant term

% ************************************************************
% empirical GSM (GMSK narrow-bandwidth pulse) model equation
% ************************************************************
f = (-500e3:1e3:500e3)+1e-3;
gaussPart = exp(-(2*f/BW) .^2);
sincPart = sin(pi*f/BW2) ./ (pi*f/BW2);
flatPart = 10^(baseline_dB/20);
H_dB = 20*log10(abs(gaussPart .* sincPart) + flatPart);

% ************************************************************
% plot the spectrum
% ************************************************************
figure(); grid on; hold on;
h = plot(f/1e6, H_dB, 'b'); set(h, 'linewidth', 2);

% ************************************************************
% plot 'a' GSM spectral emission mask
% note, this is only "an" example
% See 3GPP TS 45.005 section 4.2.1
%   "Spectrum due to the modulation and wide band noise"
% section 4.2.1.1
%   "General requirements for all types of Base stations and MS"
% note the warning regarding measuring the nominal signal level!
% ************************************************************
maskX_MHz = [0, 100e3, 200e3, 250e3, 400e3, 600e3]/1e6;
maskY_dB = [0.5, 0.5, -30, -33, -60, -60];

xlabel('f/MHz');
ylabel('PSD/dB');
title('GSM spectrum');``````

## FPGA IIR Lowpass Direct Form I Filter Generator

August 24, 2011 Coded in Python
``````# Instantiate the SIIR object.  Pass the cutoff frequency
# Fc and the sample rate Fs in Hz.  Also define the input
# and output fixed-point type.  W=(wl, iwl) where
# wl = word-length and iwl = integer word-length.  This
# example uses 23 fraction bits and 1 sign bit.
>>> from siir import SIIR
>>> flt = SIIR(Fstop=1333, Fs=48000, W=(24,0))

# Plot the response of the fixed-point coefficients
>>> plot(flt.hz, 20*log10(flt.h)

# Create a testbench and run a simulation
# (get the simulated response)
>>> from myhdl import Simulation
>>> tb = flt.TestFreqResponse(Nloops=128, Nfft=1024)
>>> Simulation(tb).run()
>>> flt.PlotResponse()

# Happy with results generate the Verilog and VHDL
>>> flt.Convert()``````