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

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

Gabriel Rivas 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()) {
            /*When there's new sample at your ADC or CODEC input*/
            /*Read the sample*/
            xin = read_sample();
            /*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

Markus Nentwig 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.
% For large frequency offsets, one would need more information on
% 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];
h = plot(maskX_MHz, maskY_dB, 'k'); set(h, 'linewidth', 3);

legend('|H(f)|', 'GSM mask example');
xlabel('f/MHz');
ylabel('PSD/dB');
title('GSM spectrum');

FPGA IIR Lowpass Direct Form I Filter Generator

Christopher Felton 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()