% ************************************************************
% 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');
% ************************************************************
% * Construct a continuous-time impulse response based on a discrete-time filter
% ************************************************************
close all; clear all;
% ************************************************************
% * example filter
% * f = [0 0.7 0.8 1]; a = [1 1 0 0];
% * b = remez(45, f, a);
% * h = b .';
% Illustrates rather clearly the difficulty of "interpolating" (in a geometric sense,
% via polynomials, splines etc) between impulse response samples
% ************************************************************
h = [2.64186706e-003 2.50796828e-003 -4.25450455e-003 4.82982358e-003 -2.51252769e-003 -2.52706568e-003 7.73569146e-003 -9.30398382e-003 4.65100223e-003 5.17152459e-003 -1.49409856e-002 1.76001904e-002 -8.65545521e-003 -9.78478603e-003 2.82103205e-002 -3.36173643e-002 1.68597129e-002 2.01914744e-002 -6.17486493e-002 8.13362871e-002 -4.80981494e-002 -8.05143565e-002 5.87677665e-001 5.87677665e-001 -8.05143565e-002 -4.80981494e-002 8.13362871e-002 -6.17486493e-002 2.01914744e-002 1.68597129e-002 -3.36173643e-002 2.82103205e-002 -9.78478603e-003 -8.65545521e-003 1.76001904e-002 -1.49409856e-002 5.17152459e-003 4.65100223e-003 -9.30398382e-003 7.73569146e-003 -2.52706568e-003 -2.51252769e-003 4.82982358e-003 -4.25450455e-003 2.50796828e-003 2.64186706e-003];
% ************************************************************
% zero-pad the impulse response of the FIR filter to 5x its original length
% (time domain)
% The filter remains functionally identical, since appending zero taps changes nothing
% ************************************************************
timePaddingFactor = 5;
n1 = timePaddingFactor * size(h, 2);
nh = size(h, 2);
nPad = n1 - nh;
nPad1 = floor(nPad / 2);
nPad2 = nPad - nPad1;
h = [zeros(1, nPad1), h, zeros(1, nPad2)];
hOrig = h;
% ************************************************************
% Determine equivalent Fourier series coefficients (frequency domain)
% assume that the impulse response is bandlimited (time-discrete signal by definition)
% and periodic (period length see above, can be extended arbitrarily)
% ************************************************************
h = fft(h); % Fourier transform time domain to frequency domain
% ************************************************************
% Evaluate the Fourier series on an arbitrarily dense grid
% this allows to resample the impulse response on an arbitrarily dense grid
% ************************************************************
% zero-pad Fourier transform
% ideal band-limited oversampling of the impulse response to n2 samples
n2 = 10 * n1;
h = [h(1:ceil(n1/2)), zeros(1, n2-n1), h(ceil(n1/2)+1:end)];
h = ifft(h); % back to time domain
% Note: One may write out the definition of ifft() and evaluate the exp(...) term at an
% arbitrary time to acquire a true continuous-time equation.
% numerical inaccuracy in (i)fft causes some imaginary part ~1e-15
assert(max(abs(imag(h))) / max(abs(real(h))) < 1e-14);
h = real(h);
% scale to maintain amplitude level
h = h * n2 / n1;
% construct x-axis [0, 1[
xOrig = linspace(-timePaddingFactor/2, timePaddingFactor/2, size(hOrig, 2) + 1); xOrig = xOrig(1:end-1);
x = linspace(-timePaddingFactor/2, timePaddingFactor/2, size(h, 2) + 1); x = x(1:end-1);
% ************************************************************
% Plot
% ************************************************************
% ... on a linear scale
% find points where original impulse response is defined (for illustration)
ixOrig = find(abs(hOrig) > 0);
figure(); grid on; hold on;
stem(xOrig(ixOrig), hOrig(ixOrig), 'k');
plot(x, h, 'b');
legend('discrete-time FIR filter', 'continuous-time equivalent');
title('equivalent discrete- and continuous-time filters');
xlabel('time (relative to discrete time IR duration)'); ylabel('h(t)');
% ... and again on a logarithmic scale
myEps = 1e-15;
figure(); grid on; hold on;
u = plot(xOrig(ixOrig), 20*log10(abs(hOrig(ixOrig)) + myEps), 'k+'); set(u, 'lineWidth', 3);
plot(x, 20*log10(abs(h) + myEps), 'b');
legend('discrete-time FIR filter', 'continuous-time equivalent');
title('equivalent discrete- and continuous-time filters');
xlabel('time (relative to discrete time IR duration)'); ylabel('|h(t)| / dB');
% *************************************************************
% Raised cosine window for OFDM intersymbol transitions
% *************************************************************
% Use the window for the guard period between adjacent symbols
% to smoothen the waveform transition, limit sidelobes and prevent
% unwanted emissions into adjacent frequency ranges
%
% reference:
% Data transmission by Frequency Division Multiplexing
% using the Discrete Fourier Transform
% S.B. Weinstein and Paul M. Ebert
% IEEE transactions on Communications technology, No. 5,
% October 1971
% *************************************************************
% Length of the guard period in samples
n=64;
index=linspace(0, 1, n+2);
index=index(2:end-1);
windowNextSymbol=(1-cos(pi*index))/2;
windowPreviousSymbol=fliplr(windowNextSymbol);
% plot
figure(); title('raised cosine window for OFDM symbol transitions');
xlabel('sample index into the guard period between symbols');
ylabel('window weight');
hold on;
plot(windowPreviousSymbol, 'r+');
plot(windowNextSymbol, 'b+');
legend({'cyclic extension of symbol i-1', 'cyclic extension of symbol i'});
% ----------------------------------------------------
% Title: Discrete Wavelet Transform
% by means of a linear convolution matrix
%
% Author: David Valencia
% UPIITA IPN
%
% Posted at: http://www.dsprelated.com/showcode/11.php
%
% Description:
% Computes the Discrete Wavelet Transform
% of n levels (and its branches) from a n sample input.
% Base filters can be of any lenght.
%
% Generalized DWT filter bank computed via the
% convolution matrix method (linear convolution)
%
% Dependencies:
%
% formfilter.m
% http://www.dsprelated.com/showcode/12.php
% upsample2.m
% http://www.dsprelated.com/showcode/10.php
%
% Revisions:
% v1.0a Commented and translated in English
%
% ----------------------------------------------------
close all; clear all; clc;
disp('== GENERALIZED MATRICIAL DWT ==')
%% STEP 1 - Base Filters
% Select type of base filters
typeofbasis = 'o';
typbior = 'bior2.2';
typor = 'db3';
% Obtain base filters
if(typeofbasis == 'b')
[Rf,Df] = biorwavf(typbior);
[h0,h1,g0,g1] = biorfilt(Df,Rf);
elseif (typeofbasis == 'o')
[h0,h1,g0,g1] = wfilters(typor);
end;
%% STEP 2 - Parameter configuration
% One can declare or recover an input vector from another
% program, but here an example vector is built for testing.
N = 16; %Number of input samples
x = (1:N)';
L = length(h0); %Base filter lenght
n_stages = 2; %Of the low-pass stage
n_branches = 2^n_stages; %Number of branches
dec_fact = n_branches; %Decimating factor
% L = Basic filters lenght (h0 ó h1)
% N = Input vector length
% n_stages = number of stages, it generates (2^n_stages) branches
hx = formfilter(n_stages, 1, h0, h1);
Lhx = length(hx);
fprintf('\nLinear convolution matrix processing\n');
%% STEP 3: Build analysis matrices
% -------------------------------------------------------------
% High-pass analysis matrix (single scale) [linear convolution]
% As the DWT has a simple high-pass branch, the matrix is easy
% to build with the base filter.
% Append zeros to each side of the vector to use it as a sliding
% part to build the matrix
movil_1 = [zeros(1,N-1),fliplr(h1),zeros(1,N-1)];
lm = length(movil_1);
n_rows = ceil((L+N-1)/2); %Compute the number of rows
W1 = zeros(n_rows,N); % Allocate memory for the high pass matrix
dec_factor = 2; %Decimate factor for the high-pass matrix
% Build matrix (with sequential shifting)
for i = 0:(n_rows-1)
W1(i+1,:) = movil_1(lm-(i*dec_factor)-N+1 : lm-(i*dec_factor)); %Rama pasaaltas
end
disp('High-pass analysis matrix: ');
W1
disp('High-pass branch results');
y1 = W1 * x
% -------------------------------------------------------------
% Low-pass several stage analysis matrix [linear convolution]
% Get an equivalent filter to get its length
hx = formfilter(n_stages, 1, h0, h1);
Lhx = length(hx);
% Compute the number of rows per branch
rowsperfilter = ceil((N+Lhx)/dec_fact);
% Allocate memory for the low-pass matrix
W0 = zeros((n_branches/2)*rowsperfilter,N);
tic; %Start cronometer
% Build low pass filter
for i = 0:(n_branches/2)-1
% Get the equivalent filter depending on the number of stages
if n_stages < 3
hx = formfilter(n_stages, (n_branches/2) - i, h0, h1);
else
hx = formfilter(n_stages, (n_branches/2) - i - 1, h0, h1);
end
% Append zeros to the vector
movil = [zeros(1,N-1),fliplr(hx),zeros(1,N-1)];
lm = length(movil);
% Shift and add the vector to the matrix
for j = 0:rowsperfilter-1
W0(i*rowsperfilter + j + 1,:) = movil(lm-(dec_fact*j)-N + 1 : lm-(dec_fact*j));
end
end
disp('Low-pass many stage analysis matrix: ');
W0
disp('Low-pass filter bank results: ');
y0 = W0*x
%% STEP 4: Recover signal
disp('Recovered signal: ');
xn1 = (W1')*y1;
xn0 = (W0')*y0;
xn = xn1+xn0
%% STEP 5: Compute error
err = (x - xn).^2;
disp('Error :');
err = (sum(err)/N);
disp(err);
fprintf('\nProcessing time = %d seg \n',toc);
fprintf('\nProgram end >>> \n');
%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
%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);
% ****************************************************************
% RMS phase error from phase noise spectrum
% Markus Nentwig, 2011
%
% - integrates RMS phase error from a phase noise power spectrum
% - generates an example signal and determines the RMS-average
% phase error (alternative method)
% ****************************************************************
function pn_snippet()
close all;
% ****************************************************************
% Phase noise spectrum model
% --------------------------
% Notes:
% - Generally, a phase noise spectrum tends to follow
% |H(f)| = c0 + c1 f^n1 + c2 f^n2 + c3 f^n3 + ...
% Therefore, linear interpolation should be performed in dB
% on a LOGARITHMIC frequency axis.
% - Single-/double sideband definition:
% The phase noise model is defined for -inf <= f <= inf
% in other words, it contributes the given noise density at
% both positive AND negative frequencies.
% Assuming a symmetrical PN spectrum, the model is evaluated
% for |f| => no need to explicitly write out the negative frequency
% side.
% ****************************************************************
% PN model
% first column:
% frequency offset
% second column:
% spectral phase noise density, dB relative to carrier in a 1 Hz
% observation bandwidth
d = [0 -80
10e3 -80
1e6 -140
9e9 -140];
f_Hz = d(:, 1) .';
g_dBc1Hz = d(:, 2) .' -3;
% get RMS phase error by integrating the power spectrum
% (alternative 1)
e1_degRMS = integratePN_degRMS(f_Hz, g_dBc1Hz)
% get RMS phase error based on a test signal
% (alternative 2)
n = 2 ^ 20;
deltaF_Hz = 2;
e2_degRMS = simulatePN_degRMS(f_Hz, g_dBc1Hz, n, deltaF_Hz)
end
% ****************************************************************
% Integrate the RMS phase error from the power spectrum
% ****************************************************************
function r = integratePN_degRMS(f1_Hz, g1_dBc1Hz)
1;
% integration step size
deltaF_Hz = 100;
% list of integration frequencies
f_Hz = deltaF_Hz:deltaF_Hz:5e6;
% interpolate spectrum on logarithmic frequency, dB scale
% unit is dBc in 1 Hz, relative to a unity carrier
fr_dB = interp1(log(f1_Hz+eps), g1_dBc1Hz, log(f_Hz+eps), 'linear');
% convert to power in 1 Hz, relative to a unity carrier
fr_pwr = 10 .^ (fr_dB/10);
% scale to integration step size
% unit: power in deltaF_Hz, relative to unity carrier
fr_pwr = fr_pwr * deltaF_Hz;
% evaluation frequencies are positive only
% phase noise is two-sided
% (note the IEEE definition: one-half the double sideband PSD)
fr_pwr = fr_pwr * 2;
% sum up relative power over all frequencies
pow_relToCarrier = sum(fr_pwr);
% convert the phase noise power to an RMS magnitude, relative to the carrier
pnMagnitude = sqrt(pow_relToCarrier);
% convert from radians to degrees
r = pnMagnitude * 180 / pi;
end
% ****************************************************************
% Generate a PN signal with the given power spectrum and determine
% the RMS phase error
% ****************************************************************
function r = simulatePN_degRMS(f_Hz, g_dBc1Hz, n, deltaF_Hz);
A = 1; % unity amplitude, arbitrary
% indices of positive-frequency FFT bins.
% Does not include DC
ixPos = 2:(n/2);
% indices of negative-frequency FFT bins.
% Does not include the bin on the Nyquist limit
assert(mod(n, 2) == 0, 'n must be even');
ixNeg = n - ixPos + 2;
% Generate signal in the frequency domain
sig = zeros(1, n);
% positive frequencies
sig(ixPos) = randn(size(ixPos)) + 1i * randn(size(ixPos));
% for purely real-valued signal: conj()
% for purely imaginary-valued signal such as phase noise: -conj()
% Note:
% Small-signals are assumed. For higher "modulation indices",
% phase noise would become a circle in the complex plane
sig(ixNeg) = -conj(sig(ixPos));
% normalize energy to unity amplitude A per bin
sig = sig * A / RMS(sig);
% normalize energy to 0 dBc in 1 Hz
sig = sig * sqrt(deltaF_Hz);
% frequency vector corresponding to the FFT bins (0, positive, negative)
fb_Hz = FFT_freqbase(n, deltaF_Hz);
% interpolate phase noise spectrum on frequency vector
% interpolate dB on logarithmic frequency axis
H_dB = interp1(log(f_Hz+eps), g_dBc1Hz, log(abs(fb_Hz)+eps), 'linear');
% dB => magnitude
H = 10 .^ (H_dB / 20);
% plot on linear f axis
figure(); hold on;
plot(fftshift(fb_Hz), fftshift(mag2dB(H)), 'k');
xlabel('f/Hz'); ylabel('dBc in 1 Hz');
title('phase noise spectrum (linear frequency axis)');
% plot on logarithmic f axis
figure(); hold on;
ix = find(fb_Hz > 0);
semilogx(fb_Hz(ix), mag2dB(H(ix)), 'k');
xlabel('f/Hz'); ylabel('dBc in 1 Hz');
title('phase noise spectrum (logarithmic frequency axis)');
% apply frequency response to signal
sig = sig .* H;
% add carrier
sig (1) = A;
% convert to time domain
td = ifft(sig);
% determine phase
% for an ideal carrier, it would be zero
% any deviation from zero is phase error
ph_deg = angle(td) * 180 / pi;
figure();
ix = 1:1000;
plot(ix / n / deltaF_Hz, ph_deg(ix));
xlabel('time / s');
ylabel('phase error / degrees');
title('phase of example signal (first 1000 samples)');
figure();
hist(ph_deg, 300);
title('histogram of example signal phase');
xlabel('phase error / degrees');
% RMS average
% note: exp(1i*x) ~ x
% as with magnitude, power/energy is proportional to phase error squared
r = RMS(ph_deg);
% knowing that the signal does not have an amplitude component, we can alternatively
% determine the RMS phase error from the power of the phase noise component
% exclude carrier:
sig(1) = 0;
% 3rd alternative to calculate RMS phase error
rAlt_deg = sqrt(sig * sig') * 180 / pi
end
% gets a vector of frequencies corresponding to n FFT bins, when the
% frequency spacing between two adjacent bins is deltaF_Hz
function fb_Hz = FFT_freqbase(n, deltaF_Hz)
fb_Hz = 0:(n - 1);
fb_Hz = fb_Hz + floor(n / 2);
fb_Hz = mod(fb_Hz, n);
fb_Hz = fb_Hz - floor(n / 2);
fb_Hz = fb_Hz * deltaF_Hz;
end
% Root-mean-square average
function r = RMS(vec)
r = sqrt(vec * vec' / numel(vec));
end
% convert a magnitude to dB
function r = mag2dB(vec)
r = 20*log10(abs(vec) + eps);
end
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
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]
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);
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;
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);
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;
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;