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
Computing the Nth Roots of a Number
function [nth_Roots] = roots_nth(Num, n, Plot)
% ROOTS_NTH computes the nth roots of Num.
%
% Inputs:
% Num is the number(s) for which the nth roots will be computed.
% Num can be real-valued or complex-valued.
% If Num has more than one element, then Num must be a row vector.
% n must be a positive integer.
% Plot is a string equal to 'Y' if a complex plane plot of
% the n complex roots of each element in Num is desired.
% Principle root is plotted as a blue square. Remaining roots
% are red dots. If input Plot string does not exist, then no
% plot is generated.
%
% Output:
% nth_Roots is a column vector (with n rows) giving the nth roots of Num.
%
% Example:
% clear, clc % clear all variables
% Num = [2+j*11, 3*exp(j*pi)]; % Two numbers
% n = 3; % find the 3rd roots
% [nth_Roots] = roots_nth(Num, n, 'Y')
%
% Results are:
% nth_Roots =
%
% 2.0000 + 1.0000i 0.7211 + 1.2490i
% -1.8660 + 1.2321i -1.4422 + 0.0000i
% -0.1340 - 2.2321i 0.7211 - 1.2490i
%
% [R. Lyons, Jan. 2013]
% Input argument check
if (nargin<2),
disp('Not enough input arguments!!')
return
end;
if (nargin==2)
Plot == 'N'
end
Num_of_Columns = max(size(Num)); % How many elements are in Num?
for Column_Loop = 1:Num_of_Columns;
Mag(Column_Loop) = abs(Num(Column_Loop));
Angle_Degrees(Column_Loop) = angle(Num(Column_Loop))*180/pi;
for k = 1:n
Root_Mag = Mag(Column_Loop)^(1/n);
Root_Angle_Deg(k) = Angle_Degrees(Column_Loop)/n + (k-1)*360/n;
nth_Roots(k,Column_Loop) = Root_Mag*exp(j*Root_Angle_Deg(k)*pi/180);
end
% Plot roots, on complex plane(s), if necessary
if Plot == 'Y'
figure
% Plot unit circle
Num_Pts = 200; % Number of circle points
Index = 0 : Num_Pts-1;
Alpha = 1/Num_Pts;
Points = Root_Mag*exp(j*2*pi*Alpha*Index);
plot(Points, ':k','markersize',5)
grid on
hold on
% Plot principle root (square blue dot) and list results in plot
axis([-1.1*Root_Mag, 1.1*Root_Mag, -1.1*Root_Mag, 1.1*Root_Mag]);
plot(Root_Mag*cos(Root_Angle_Deg(1)*pi/180), ...
Root_Mag*sin(Root_Angle_Deg(1)*pi/180), 'bs');
text(-Root_Mag/2, 3*Root_Mag/4, ['Magnitudes = ',num2str(Root_Mag)]);
text(-Root_Mag/2, 3*Root_Mag/4-Root_Mag/10, ...
['Principle Angle (deg.) = ',num2str(Root_Angle_Deg(1))]);
% Plot remaining roots as red dots
for k = 2:n
plot(Root_Mag*cos(Root_Angle_Deg(k)*pi/180), ...
Root_Mag*sin(Root_Angle_Deg(k)*pi/180), 'ro');
text(-Root_Mag/2, 3*Root_Mag/4-k*Root_Mag/10, ...
[' Angle ',num2str(k),' (deg.) = ', ...
num2str(Root_Angle_Deg(k))]);
end % End k loop
xlabel('Real'); ylabel('Imag.');
hold off
end % End 'if Plot == 'Y''
end % End Column_Loop
2D Matrix Downsample
function [ Y ] = downsample2d(M)
x = 2;
y = 2;
N = downsample(M,x);
N = N';
P = downsample(N,y);
P = P'
end
Spectrum analyzer
% ************************************
% spectrum analyzer
% Markus Nentwig, 10.12.2011
%
% Note: If the signal is not cyclic, apply conventional windowing before
% calling spectrumAnalyzer
% ************************************
function satest()
close all;
rate_Hz = 30.72e6; % sampling rate
n = 100000; % number of samples
noise = randn(1, n);
pulse = zeros(1, n); pulse(1) = 1;
% ************************************
% example (linear frequency scale, RBW filter)
% ************************************
tmp1 = normalize(RRC_filter('signal', noise, 'rate_Hz', rate_Hz, 'alpha', 0.22, 'BW_Hz', 3.84e6, 'fCenter_Hz', 5e6));
% pRef = 0.001: a normalized signal will show as 0 dBm = 0.001 W
saPar = struct('rate_Hz', rate_Hz, 'fig', 1, 'pRef_W', 1e-3, 'RBW_power_Hz', 3.84e6, 'filter', 'brickwall', 'plotStyle', 'b-');
spectrumAnalyzer(saPar, 'signal', tmp1, 'RBW_window_Hz', 1e3);
spectrumAnalyzer(saPar, 'signal', tmp1, 'RBW_window_Hz', 30e3, 'plotStyle', 'k-');
spectrumAnalyzer(saPar, 'signal', tmp1, 'RBW_window_Hz', 300e3, 'plotStyle', 'r-');
spectrumAnalyzer(saPar, 'signal', tmp1, 'RBW_window_Hz', 30e3, 'filter', 'gaussian', 'plotStyle', 'g-');
legend('1k RBW filter', '30k RBW filter', '300k RBW filter', '30k Gaussian filter (meas. instrument model)');
title('WCDMA-like signal');
% ************************************
% example (logarithmic frequency scale,
% no RBW filter)
% ************************************
fPar = struct('rate_Hz', rate_Hz, ...
'chebyOrder', 6, ...
'chebyRipple_dB', 1, ...
'fCorner_Hz', 1e5);
saPar = struct('rate_Hz', rate_Hz, ...
'pRef_W', 1e-3, ...
'fMin_Hz', 1e3, ...
'RBW_power_Hz', 1e5, ...
'RBW_window_Hz', 1e3, ...
'filter', 'none', ...
'plotStyle', 'b-', ...
'logscale', true, ...
'fig', 2, ...
'noisefloor_dB', -300);
tmp1 = normalize(IIR_filterExample('signal', noise, fPar));
tmp2 = normalize(IIR_filterExample('signal', pulse, fPar));
spectrumAnalyzer('signal', tmp1, saPar);
spectrumAnalyzer('signal', tmp2, saPar, 'plotStyle', 'k-');
end
function sig = normalize(sig)
sig = sig / sqrt(sig * sig' / numel(sig));
end
% ************************************
% calculates the frequency that corresponds to each
% FFT bin (negative, zero, positive)
% ************************************
function fb_Hz = FFT_frequencyBasis(n, rate_Hz)
fb = 0:(n - 1);
fb = fb + floor(n / 2);
fb = mod(fb, n);
fb = fb - floor(n / 2);
fb = fb / n; % now [0..0.5[, [-0.5..0[
fb_Hz = fb * rate_Hz;
end
% ************************************
% root-raised cosine filter (to generate
% example signal)
% ************************************
function sig = RRC_filter(varargin)
def = struct('fCenter_Hz', 0);
p = vararginToStruct(def, varargin);
n = numel(p.signal);
fb_Hz = FFT_frequencyBasis(n, p.rate_Hz);
% frequency relative to cutoff frequency
c = abs((fb_Hz - p.fCenter_Hz) / (p.BW_Hz / 2));
% c=-1 for lower end of transition region
% c=1 for upper end of transition region
c=(c-1) / p.alpha;
% clip to -1..1 range
c=min(c, 1);
c=max(c, -1);
% raised-cosine filter
H = 1/2+cos(pi/2*(c+1))/2;
% root-raised cosine filter
H = sqrt(H);
% apply filter
sig = ifft(fft(p.signal) .* H);
end
% ************************************
% continuous-time filter example
% ************************************
function sig = IIR_filterExample(varargin)
p = vararginToStruct(varargin);
% design continuous-time IIR filter
[IIR_b, IIR_a] = cheby1(p.chebyOrder, p.chebyRipple_dB, 1, 's');
% evaluate on the FFT frequency grid
fb_Hz = FFT_frequencyBasis(numel(p.signal), p.rate_Hz);
% Laplace domain operator, normalized to filter cutoff frequency
% (the above cheby1 call designs for unity cutoff)
s = 1i * fb_Hz / p.fCorner_Hz;
% polynomial in s
H_num = polyval(IIR_b, s);
H_denom = polyval(IIR_a, s);
H = H_num ./ H_denom;
% apply filter
sig = ifft(fft(p.signal) .* H);
end
% ----------------------------------------------------------------------
% optionally: Move the code below into spectrumAnalyzer.m
function [fr, fb, handle] = spectrumAnalyzer(varargin)
def = struct(...
'noisefloor_dB', -150, ...
'filter', 'none', ...
'logscale', false, ...
'NMax', 10000, ...
'freqUnit', 'MHz', ...
'fig', -1, ...
'plotStyle', 'b-', ...
'originMapsTo_Hz', 0 ...
);
p = vararginToStruct(def, varargin);
signal = p.signal;
handle=nan; % avoid warning
% A resolution bandwidth value of 'sine' sets the RBW to the FFT bin spacing.
% The power of a pure sine wave now shows correctly from the peak in the spectrum (unity => 0 dB)
singleBinMode=strcmp(p.RBW_power_Hz, 'sine');
nSamples = numel(p.signal);
binspacing_Hz = p.rate_Hz / nSamples;
windowBW=p.RBW_window_Hz;
noisefloor=10^(p.noisefloor_dB/10);
% factor in the scaling factor that will be applied later for conversion to
% power in RBW
if ~singleBinMode
noisefloor = noisefloor * binspacing_Hz / p.RBW_power_Hz;
end
% fr: "f"requency "r"esponse (plot y data)
% fb: "f"requency "b"ase (plot x data)
fr = fft(p.signal);
scale_to_dBm=sqrt(p.pRef_W/0.001);
% Normalize amplitude to the number of samples that contribute
% to the spectrum
fr=fr*scale_to_dBm/nSamples;
% magnitude
fr = fr .* conj(fr);
[winLeft, winRight] = spectrumAnalyzerGetWindow(p.filter, singleBinMode, p.RBW_window_Hz, binspacing_Hz);
% winLeft is the right half of the window, but it appears on the
% left side of the FFT space
winLen=0;
if ~isempty(winLeft)
% zero-pad the power spectrum in the middle with a length
% equivalent to the window size.
% this guarantees that there is enough bandwidth for the filter!
% one window length is enough, the spillover from both sides overlaps
% Scale accordingly.
winLen=size(winLeft, 2)+size(winRight, 2);
% note: fr is the power shown in the plot, NOT a frequency
% domain representation of a signal.
% There is no need to renormalize because of the length change
center=floor(nSamples/2)+1;
rem=nSamples-center;
fr=[fr(1:center-1), zeros(1, winLen-1), fr(center:end)];
% construct window with same length as fr
win=zeros(size(fr));
win(1:1+size(winLeft, 2)-1)=winLeft;
win(end-size(winRight, 2)+1:end)=winRight;
assert(isreal(win));
assert(isreal(fr));
assert(size(win, 2)==size(fr, 2));
% convolve using FFT
win=fft(win);
fr=fft(fr);
fr=fr .* win;
fr=ifft(fr);
fr=real(fr); % chop off roundoff error imaginary part
fr=max(fr, 0); % chop off small negative numbers
% remove padding
fr=[fr(1:center-1), fr(end-rem+1:end)];
end
% ************************************
% build frequency basis and rotate 0 Hz to center
% ************************************
fb = FFT_frequencyBasis(numel(fr), p.rate_Hz);
fr = fftshift(fr);
fb = fftshift(fb);
if false
% use in special cases (very long signals)
% ************************************
% data reduction:
% If a filter is used, details smaller than windowBW are
% suppressed already by the filter, and using more samples
% gives no additional information.
% ************************************
if numel(fr) > p.NMax
switch(p.filter)
case 'gaussian'
% 0.2 offset from the peak causes about 0.12 dB error
incr=floor(windowBW/binspacing_Hz/5);
case 'brickwall'
% there is no error at all for a peak
incr=floor(windowBW/binspacing_Hz/3);
case 'none'
% there is no filter, we cannot discard data at this stage
incr=-1;
end
if incr > 1
fr=fr(1:incr:end);
fb=fb(1:incr:end);
end
end
end
% ************************************
% data reduction:
% discard beyond fMin / fMax, if given
% ************************************
indexMin = 1;
indexMax = numel(fb);
flag=0;
if isfield(p, 'fMin_Hz')
indexMin=min(find(fb >= p.fMin_Hz));
flag=1;
end
if isfield(p, 'fMax_Hz')
indexMax=max(find(fb <= p.fMax_Hz));
flag=1;
end
if flag
fb=fb(indexMin:indexMax);
fr=fr(indexMin:indexMax);
end
if p.NMax > 0
if p.logscale
% ************************************
% Sample number reduction for logarithmic
% frequency scale
% ************************************
assert(isfield(p, 'fMin_Hz'), 'need fMin_Hz in logscale mode');
assert(p.fMin_Hz > 0, 'need fMin_Hz > 0 in logscale mode');
if ~isfield(p, 'fMax_Hz')
p.fMax_Hz = p.rate_Hz / 2;
end
% averaged output arrays
fbOut=zeros(1, p.NMax-1);
frOut=zeros(1, p.NMax-1);
% candidate output frequencies (it's not certain yet
% that there is actually data)
s=logspace(log10(p.fMin_Hz), log10(p.fMax_Hz), p.NMax);
f1=s(1);
nextStartIndex=max(find(fb < f1));
if isempty(nextStartIndex)
nextStartIndex=1;
end
% iterate through all frequency regions
% collect data
% average
for index=2:size(s, 2)
f2=s(index);
endIndex=max(find(fb < f2));
% number of data points in bin
n=endIndex-nextStartIndex+1;
if n > 0
% average
ix=nextStartIndex:endIndex;
fbOut(index-1)=sum(fb(ix))/n;
frOut(index-1)=sum(fr(ix))/n;
nextStartIndex=endIndex+1;
else
% mark point as invalid (no data)
fbOut(index-1)=nan;
end
end
% remove bins where no data point contributed
ix=find(~isnan(fbOut));
fbOut=fbOut(ix);
frOut=frOut(ix);
fb=fbOut;
fr=frOut;
else
% ************************************
% Sample number reduction for linear
% frequency scale
% ************************************
len=size(fb, 2);
decim=ceil(len/p.NMax); % one sample overlength => decim=2
if decim > 1
% truncate to integer multiple
len=floor(len / decim)*decim;
fr=fr(1:len);
fb=fb(1:len);
fr=reshape(fr, [decim, len/decim]);
fb=reshape(fb, [decim, len/decim]);
if singleBinMode
% apply peak hold over each segment (column)
fr=max(fr);
else
% apply averaging over each segment (column)
fr = sum(fr) / decim;
end
fb=sum(fb)/decim; % for each column the average
end % if sample reduction necessary
end % if linear scale
end % if sample number reduction
% ************************************
% convert magnitude to dB.
% truncate very small values
% using an artificial noise floor
% ************************************
fr=(10*log10(fr+noisefloor));
if singleBinMode
% ************************************
% The power reading shows the content of each
% FFT bin. This is accurate for single-frequency
% tones that fall exactly on the frequency grid
% (an integer number of cycles over the signal length)
% ************************************
else
% ************************************
% convert sensed power density from FFT bin spacing
% to resolution bandwidth
% ************************************
fr=fr+10*log10(p.RBW_power_Hz/binspacing_Hz);
end
% ************************************
% Post-processing:
% Translate frequency axis to RF
% ************************************
fb = fb + p.originMapsTo_Hz;
% ************************************
% convert to requested units
% ************************************
switch(p.freqUnit)
case 'Hz'
case 'kHz'
fb = fb / 1e3;
case 'MHz'
fb = fb / 1e6;
case 'GHz'
fb = fb / 1e9;
otherwise
error('unsupported frequency unit');
end
% ************************************
% Plot (if requested)
% ************************************
if p.fig > 0
% *************************************************************
% title
% *************************************************************
if isfield(p, 'title')
t=['"', p.title, '";'];
else
t='';
end
switch(p.filter)
case 'gaussian'
t=[t, ' filter: Gaussian '];
case 'brickwall'
t=[t, ' filter: ideal bandpass '];
case 'none'
t=[t, ' filter: none '];
otherwise
assert(0)
end
if ~strcmp(p.filter, 'none')
t=[t, '(', mat2str(p.RBW_window_Hz), ' Hz)'];
end
if isfield(p, 'pNom_dBm')
% *************************************************************
% show dB relative to a given absolute power in dBm
% *************************************************************
fr=fr-p.pNom_dBm;
yUnit='dB';
else
yUnit='dBm';
end
% *************************************************************
% plot
% *************************************************************
figure(p.fig);
if p.logscale
handle = semilogx(fb, fr, p.plotStyle);
else
handle = plot(fb, fr, p.plotStyle);
end
hold on; % after plot, otherwise prevents logscale
if isfield(p, 'lineWidth')
set(handle, 'lineWidth', p.lineWidth);
end
% *************************************************************
% axis labels
% *************************************************************
xlabel(p.freqUnit);
if singleBinMode
ylabel([yUnit, ' (continuous wave)']);
else
ylabel([yUnit, ' in ', mat2str(p.RBW_power_Hz), ' Hz integration BW']);
end
title(t);
% *************************************************************
% adapt y range, if requested
% *************************************************************
y=ylim();
y1=y(1); y2=y(2);
rescale=false;
if isfield(p, 'yMin')
y(1)=p.yMin;
rescale=true;
end
if isfield(p, 'yMax')
y(2)=p.yMax;
rescale=true;
end
if rescale
ylim(y);
end
end
end
function [winLeft, winRight] = spectrumAnalyzerGetWindow(filter, singleBinMode, RBW_window_Hz, binspacing_Hz)
switch(filter)
case 'gaussian'
% construct Gaussian filter
% -60 / -3 dB shape factor 4.472
nRBW=6;
nOneSide=ceil(RBW_window_Hz/binspacing_Hz*nRBW);
filterBase=linspace(0, nRBW, nOneSide);
winLeft=exp(-(filterBase*0.831) .^2);
winRight=fliplr(winLeft(2:end)); % omit 0 Hz bin
case 'brickwall'
nRBW=1;
n=ceil(RBW_window_Hz/binspacing_Hz*nRBW);
n1 = floor(n/2);
n2 = n - n1;
winLeft=ones(1, n1);
winRight=ones(1, n2);
case 'none'
winLeft=[];
winRight=[];
otherwise
error('unknown RBW filter type');
end
% the window is not supposed to shift the spectrum.
% Therefore, the first bin is always in winLeft (0 Hz):
if size(winLeft, 2)==1 && isempty(winRight)
% there is no use to convolve with one-sample window
% it's always unity
winLeft=[];
winRight=[];
tmpwin=[];
end
if ~isempty(winLeft)
% (note: it is not possible that winRight is empty, while winLeft is not)
if singleBinMode
% normalize to unity at 0 Hz
s=winLeft(1);
else
% normalize to unity area under the filter
s=sum(winLeft)+sum(winRight);
end
winLeft=winLeft / s;
winRight=winRight / s;
end
end
% *************************************************************
% helper function: Parse varargin argument list
% allows calling myFunc(A, A, A, ...)
% where A is
% - key (string), value (arbitrary) => result.key = value
% - a struct => fields of A are copied to result
% - a cell array => recursive handling using above rules
% *************************************************************
function r = vararginToStruct(varargin)
% note: use of varargin implicitly packs the caller's arguments into a cell array
% that is, calling vararginToStruct('hello') results in
% varargin = {'hello'}
r = flattenCellArray(varargin, struct());
end
function r = flattenCellArray(arr, r)
ix=1;
ixMax = numel(arr);
while ix <= ixMax
e = arr{ix};
if iscell(e)
% cell array at 'key' position gets recursively flattened
% becomes struct
r = flattenCellArray(e, r);
elseif ischar(e)
% string => key.
% The following entry is a value
ix = ix + 1;
v = arr{ix};
% store key-value pair
r.(e) = v;
elseif isstruct(e)
names = fieldnames(e);
for ix2 = 1:numel(names)
k = names{ix2};
r.(k) = e.(k);
end
else
e
assert(false)
end
ix=ix+1;
end % while
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;
FFT for arbitrary time and frequency grids (Chirp-z Transform)
function Demo
while 1
%Frequency grid with arbitrary start value, spacing, and number of elements. ...
%Mathematically, the grid is
%
% fo+Df*m, m=0, 1,...,M-1.
%
fo = 10*rand(1)-0.5; %Start value.
Df = rand(1)+0.01; %Frequency spacing.
M = round(rand(1)*1000); %Number of frequencies.
%Time grid with arbitrary start value, spacing, and number of elements. Mathematically, ...
%the grid is
%
% to+Dt*n, n=0, 1,...,N-1.
%
to = 10*rand(1)-0.5; %Start value; (instant of first sample).
Dt = rand(1)+0.01; %Time spacing.
N = round(rand(1)*1000); %Number of samples.
%Random vector of samples.
x = randn(1,N)+1i*randn(1,N);
%x(n) is the sample at instant to+Dt*(n-1), n=1,2,...,N.
%We proceed to compute the DFT of x at frequencies fo+m*Df, m=0, 1,..., M-1.
%Compute DFT using the direct and chirp-z methods.
tic, XD = DirectGridDFT(x,to,Dt,fo,Df,M); tmD = toc;
tic, XF = chirpzDFT(x,to,Dt,fo,Df,M); tmF = toc;
disp(['Timing direct method: ', num2str(tmD), ' sec.'])
disp(['Timing chirp z: ',num2str(tmF),' sec.'])
disp(['Relative difference: ', num2str(max(abs(XD-XF))/max(abs(XD)))])
disp(' ')
input('(Press a key for another trial or Ctrl-C) ')
disp(' ')
end
function X = chirpzDFT(x,to,Dt,fo,Df,M)
%Author: J. Selva.
%Date: November 2011.
%
%This function computes the DFT for two regular time and frequency grids with arbitrary
%starting points and spacings, using the Chirp-z Transform. Specifically, it computes
%
% N
% X(k) = sum x(n)*exp(-j*2*pi*(fo+Df*(k-1))*(to+Dt*(n-1))), 1 <= k <= M.
% n=1
%
%Input parameters:
%
% x: Signal samples.
% to: Instant of first sample in x.
% Dt: Time spacing between consecutive samples of x.
% fo: First frequency in which the Fourier sum is computed.
% Df: Spacing of the desired frequency grid.
% M: Desired number of output samples.
%
%The algorithm is explained in Sec. 9.6.2, page 656 of
%
% A.V. Oppenheim, R. W. Schafer, J. R. Buck, "Discrete-time signal processing," second
% edition, 1998.
%
x = x(:).';
N = numel(x);
P = 2^ceil(log2(M+N-1));
%The next three lines do not depend on the vector x, and so they can be pre-computed if
%the time and frequency grids do not change, (i.e, for computing the transform of
%additional vectors). Thus, this algorithm just involves two FFTs for fixed grids and
%three if the grids change.
a = exp(-1i*2*pi*((1:N)*Dt*(fo-Df)+Dt*Df*(1:N).^2/2));
b = exp(-1i*2*pi*((to-Dt)*(fo+(0:M-1)*Df)+Dt*Df*(1:M).^2/2));
phi = fft(exp(1i*2*pi*Dt*Df*(1-N:M-1).^2/2),P); %FFT of chirp pulse.
%Weigh x using a and perform FFT convolution with phi.
X = ifft( fft(x.*a,P) .* phi );
%Truncate the convolution tails and weigh using b.
X = X(N:M+N-1) .* b;
function X = DirectGridDFT(x,to,Dt,fo,Df,M)
%Direct (and slow) version of the Fourier sum with arbitrary and regular time and
%frequency grids. Used for validation of chirpzDFT.
x = x(:).';
N = length(x);
X = zeros(1,M);
for k = 1:M
for n = 1:N
X(k) = X(k) + x(n)*exp(-1i*2*pi*(to+Dt*(n-1))*(fo+Df*(k-1)));
end
end
Waterfall Spectrograph
function waterfallspect(s, fs, spectsize, spectnum)
%WATERFALLSPECT Display spectrum of signal s as a 3-D plot.
% s - input signal, use wavload to load a WAV file
% fs - Sampling frequency (Hz).
% spectsize - FFT window size
% spectnum - number of windows to analyze
frequencies = [0:fs/spectsize:fs/2];
offset = floor((length(s)-spectsize)/spectnum);
for i=0:(spectnum-1)
start = i*offset;
A = abs(fft(s((1+start):(start+spectsize))));
magnitude(:,(i+1)) = A(1:spectnum/2+1);
end
waterfall(frequencies, 0:(spectnum-1), magnitude');
xlabel('frequency');
ylabel('time');
zlabel('magnitude');
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);
Beat Notes
% Filename: Beat_Frequency.m
%
% [Richard Lyons, Feb. 2013]
clear, clc
Fs = 8192; % Sample rate of dig. samples
N = 8192; % Number of time samples
n = 0:N-1;
Wave_1 = sin(2*pi*210*n/Fs); % First tone, 210 Hz
Wave_2 = sin(2*pi*200*n/Fs); % Second tone, 200 Hz
% Plot the two tones
figure(1)
plot(n/Fs, Wave_1, '-b')
ylabel('200 Hz'); xlabel('Time (sec.)');
hold on
plot(n/Fs, Wave_2, '-r')
axis([0, 0.05, -1.2, 1.5]);
ylabel('Input tones'); xlabel('Time (sec.)');
title('red = 200 Hz tone, blue = 210 Hz tone');
grid on, zoom on
hold off
Product = Wave_1.*Wave_2;
Sum = Wave_1 + Wave_2;
% Plot the tones' product and sum
figure(2)
subplot(2,1,1)
plot(n/Fs, Product, '-b'),
ylabel('Product'); xlabel('Time (sec.)');
grid on, zoom on
hold on
Red_Curve = 0.5*cos(2*pi*10*n/Fs) + 0.5; % Used for plotting only
plot(n/Fs, Red_Curve, '-r')
axis([0, 0.3, -1.25, 1.5]);
hold off
grid on, zoom on
subplot(2,1,2)
plot(n/Fs, Sum, '-b')
hold on
Red_Curve = 2*cos(2*pi*5*n/Fs); % Used for plotting only
plot(n/Fs, Red_Curve, '-r')
axis([0, 0.3, -2.4, 3]);
hold off
ylabel('Sum'); xlabel('Time (sec.)');
grid on, zoom on
% Play all the signals
sound(Wave_1, Fs)
pause(1.2)
sound(Wave_2, Fs)
pause(1.2)
sound(Product, Fs)
pause(1.2)
sound(Sum, Fs)
% Spec analysis of the "Sum" signal
Spec = fft(Sum);
Spec_Mag = abs(Spec);
Freq = n*Fs/N; % Freq vector in Hz
figure (3) % Plot positive-freq spec amg
plot(Freq(1:N/16), Spec_Mag(1:N/16))
title('Spec Mag of "Sum" signal')
ylabel('Mag.'), xlabel('Hz')
grid on, zoom on
Computing the Nth Roots of a Number
function [nth_Roots] = roots_nth(Num, n, Plot)
% ROOTS_NTH computes the nth roots of Num.
%
% Inputs:
% Num is the number(s) for which the nth roots will be computed.
% Num can be real-valued or complex-valued.
% If Num has more than one element, then Num must be a row vector.
% n must be a positive integer.
% Plot is a string equal to 'Y' if a complex plane plot of
% the n complex roots of each element in Num is desired.
% Principle root is plotted as a blue square. Remaining roots
% are red dots. If input Plot string does not exist, then no
% plot is generated.
%
% Output:
% nth_Roots is a column vector (with n rows) giving the nth roots of Num.
%
% Example:
% clear, clc % clear all variables
% Num = [2+j*11, 3*exp(j*pi)]; % Two numbers
% n = 3; % find the 3rd roots
% [nth_Roots] = roots_nth(Num, n, 'Y')
%
% Results are:
% nth_Roots =
%
% 2.0000 + 1.0000i 0.7211 + 1.2490i
% -1.8660 + 1.2321i -1.4422 + 0.0000i
% -0.1340 - 2.2321i 0.7211 - 1.2490i
%
% [R. Lyons, Jan. 2013]
% Input argument check
if (nargin<2),
disp('Not enough input arguments!!')
return
end;
if (nargin==2)
Plot == 'N'
end
Num_of_Columns = max(size(Num)); % How many elements are in Num?
for Column_Loop = 1:Num_of_Columns;
Mag(Column_Loop) = abs(Num(Column_Loop));
Angle_Degrees(Column_Loop) = angle(Num(Column_Loop))*180/pi;
for k = 1:n
Root_Mag = Mag(Column_Loop)^(1/n);
Root_Angle_Deg(k) = Angle_Degrees(Column_Loop)/n + (k-1)*360/n;
nth_Roots(k,Column_Loop) = Root_Mag*exp(j*Root_Angle_Deg(k)*pi/180);
end
% Plot roots, on complex plane(s), if necessary
if Plot == 'Y'
figure
% Plot unit circle
Num_Pts = 200; % Number of circle points
Index = 0 : Num_Pts-1;
Alpha = 1/Num_Pts;
Points = Root_Mag*exp(j*2*pi*Alpha*Index);
plot(Points, ':k','markersize',5)
grid on
hold on
% Plot principle root (square blue dot) and list results in plot
axis([-1.1*Root_Mag, 1.1*Root_Mag, -1.1*Root_Mag, 1.1*Root_Mag]);
plot(Root_Mag*cos(Root_Angle_Deg(1)*pi/180), ...
Root_Mag*sin(Root_Angle_Deg(1)*pi/180), 'bs');
text(-Root_Mag/2, 3*Root_Mag/4, ['Magnitudes = ',num2str(Root_Mag)]);
text(-Root_Mag/2, 3*Root_Mag/4-Root_Mag/10, ...
['Principle Angle (deg.) = ',num2str(Root_Angle_Deg(1))]);
% Plot remaining roots as red dots
for k = 2:n
plot(Root_Mag*cos(Root_Angle_Deg(k)*pi/180), ...
Root_Mag*sin(Root_Angle_Deg(k)*pi/180), 'ro');
text(-Root_Mag/2, 3*Root_Mag/4-k*Root_Mag/10, ...
[' Angle ',num2str(k),' (deg.) = ', ...
num2str(Root_Angle_Deg(k))]);
end % End k loop
xlabel('Real'); ylabel('Imag.');
hold off
end % End 'if Plot == 'Y''
end % End Column_Loop
Delay estimation revisited
% ****************************************************************
% find sub-sample delay and scaling factor between two cyclic signals
% to maximize crosscorrelation
% Markus Nentwig, 120609_v1.1
% ****************************************************************
function iterDelayEstDemo();
close all;
n = 1024;
% n = 1024 * 256; disp('*** test: long signal enabled ***');
% ****************************************************************
% random signal
% ****************************************************************
fd = randn(1, n) + 1i * randn(1, n);
% ****************************************************************
% lowpass filter
% ****************************************************************
f = binFreq(n);
fd(abs(f) > 0.045) = 0;
s1 = real(ifft(fd)) * sqrt(n);
% ****************************************************************
% create delayed 2nd signal
% ****************************************************************
dTest_samples = 12.3456;
cTest = 1.23456;
% cTest = cTest + 1i; disp('*** test: complex coeff enabled ***');
% cTest = -cTest; disp('*** test: negative coeff enabled ***');
s2 = cTest * cs_delay(s1, 1, dTest_samples);
% s2 = s2 + 0.5*randn(size(s2)); disp('*** test: noise enabled ***');
% ****************************************************************
% estimate delay
% ****************************************************************
[delay_samples, coeff] = iterDelayEst(s1, s2);
% ****************************************************************
% correct it
% ****************************************************************
s2a = cs_delay(s1, 1, delay_samples);
s2b = s2a * coeff;
figure(); hold on;
h = plot(real(s1), 'k'); set(h, 'lineWidth', 3);
h = plot(real(s2), 'b'); set(h, 'lineWidth', 3);
h = plot(real(s2a), 'r'); set(h, 'lineWidth', 1);
h = plot(real(s2b), 'm'); set(h, 'lineWidth', 1);
xlim([1, numel(s1)]);
xlabel('samples');
legend('s1', 's2', 's2 un-delayed', 's2 un-delayed and scaled');
title('test signals');
format long;
disp('nominal delay of s2 relative to s1')';
dTest_samples
disp('iterDelayEst() returned:');
delay_samples
disp('original scaling factor:');
cTest
disp('estimated scaling factor:');
coeff
end
% ****************************************************************
% estimates delay and scaling factor
% ****************************************************************
function [delay_samples, coeff] = iterDelayEst(s1, s2)
s1 = s1(:) .'; % force row vectors
s2 = s2(:) .';
rflag = isreal(s1) && isreal(s2);
n = numel(s1);
halfN = floor(n/2);
assert(numel(s2) == n, 'signals must have same length');
% ****************************************************************
% constants
% ****************************************************************
% exit if uncertainty below threshold
thr_samples = 1e-7;
% exit after fixed number of iterations
nIter = 25;
% frequency domain representation of signals
fd1 = fft(s1);
fd2 = fft(s2);
% first round: No delay was applied
tau = [];
fd2Tau = fd2; % delayed s2 in freq. domain
% frequency corresponding to each FFT bin -0.5..0.5
f = binFreq(n);
% uncertainty plot data
e = [];
% normalization factor
nf = sqrt((fd1 * fd1') * (fd2 * fd2')) / n; % normalizes to 1
% search window:
% known maximum and two surrounding points
x1 = -1;
x2 = -1;
x3 = -1;
y1 = -1;
y2 = -1;
y3 = -1;
% ****************************************************************
% iteration loop
% ****************************************************************
for count = 1:nIter
% ****************************************************************
% crosscorrelation with time-shifted signal
% ****************************************************************
xcorr = abs(ifft(fd2Tau .* conj(fd1)))/ nf;
% ****************************************************************
% detect peak
% ****************************************************************
if isempty(tau)
% ****************************************************************
% startup
% initialize with three adjacent bins around peak
% ****************************************************************
ix = find(xcorr == max(xcorr));
ix = ix(1); % use any, if multiple bitwise identical peaks
% indices of three bins around peak
ixLow = mod(ix-1-1, n) + 1; % one below
ixMid = ix;
ixHigh = mod(ix-1+1, n) + 1; % one above
% delay corresponding to the three bins
tauLow = mod(ixLow -1 + halfN, n) - halfN;
tauMid = mod(ixMid -1 + halfN, n) - halfN;
tauHigh = mod(ixHigh -1 + halfN, n) - halfN;
% crosscorrelation value for the three bins
xcLow = xcorr(ixLow);
xcMid = xcorr(ixMid);
xcHigh = xcorr(ixHigh);
x1 = tauLow;
x2 = tauMid;
x3 = tauHigh;
y1 = xcLow;
y2 = xcMid;
y3 = xcHigh;
else
% ****************************************************************
% only main peak at first bin is of interest
% ****************************************************************
tauMid = tau;
xcMid = xcorr(1);
if xcMid > y2
% ****************************************************************
% improve midpoint
% ****************************************************************
if tauMid > x2
% midpoint becomes lower point
x1 = x2;
y1 = y2;
else
% midpoint becomes upper point
x3 = x2;
y3 = y2;
end
x2 = tauMid;
y2 = xcMid;
elseif tauMid < x2
% ****************************************************************
% improve low point
% ****************************************************************
assert(tauMid >= x1); % bitwise identical is OK
assert(tauMid > x1 || xcMid > y1); % expect improvement
x1 = tauMid;
y1 = xcMid;
elseif tauMid > x2
% ****************************************************************
% improve high point
% ****************************************************************
assert(tauMid <= x3); % bitwise identical is OK
assert((tauMid < x3) || (xcMid > y3)); % expect improvement
x3 = tauMid;
y3 = xcMid;
else
assert(false, '?? evaluated for existing tau ??');
end
end
% ****************************************************************
% calculate uncertainty (window width)
% ****************************************************************
eIter = abs(x3 - x1);
e = [e, eIter];
if eIter < thr_samples
% disp('threshold reached, exiting');
break;
end
if y1 == y2 || y2 == y3
% reached limit of numerical accuracy on one side
usePoly = 0;
else
% ****************************************************************
% fit 2nd order polynomial and find maximum
% ****************************************************************
num = (x2^2-x1^2)*y3+(x1^2-x3^2)*y2+(x3^2-x2^2)*y1;
denom = (2*x2-2*x1)*y3+(2*x1-2*x3)*y2+(2*x3-2*x2)*y1;
if denom ~= 0
tau = num / denom;
% is the point within [x1, x3]?
usePoly = ((tau > x1) && (tau < x3));
else
usePoly = 0;
end
end
if ~usePoly
% revert to linear interpolation on the side with the
% less-accurate outer sample
% Note: There is no guarantee that the side with the more accurate
% outer sample is the right one, as the samples aren't
% placed on a regular grid!
% Therefore, iterate to improve the "worse" side, which will
% eventually become the "better side", and iteration converges.
tauLow = (x1 + x2) / 2;
tauHigh = (x2 + x3) / 2;
if y1 < y3
o = [tauLow, tauHigh];
else
o = [tauHigh, tauLow];
end
% don't choose point that is identical to one that is already known
tau = o(1);
if tau == x1 || tau == x2 || tau == x3
tau = o(2);
if tau == x1 || tau == x2 || tau == x3
break;
end
end
end
% ****************************************************************
% advance 2nd signal according to location of maximum
% phase shift in frequency domain - delay in time domain
% ****************************************************************
fd2Tau = fd2 .* exp(2i * pi * f * tau);
end % for
% ****************************************************************
% plot the uncertainty (window size) over the number of iterations
% ****************************************************************
if true
figure(); semilogy(e, '+-'); grid on;
xlabel('iteration');
title('uncertainty in delay');
end
% ****************************************************************
% the delay estimate is the final location of the delay that
% maximized crosscorrelation (center of window).
% ****************************************************************
delay_samples = x2;
% ****************************************************************
% Coefficient: Turn signal 1 into signal 2
% ****************************************************************
coeff = fd2Tau * fd1' ./ (fd1 * fd1');
% ****************************************************************
% chop roundoff error, if input signals are known to be
% real-valued.
% ****************************************************************
if rflag
coeff = real(coeff);
end
end
% ****************************************************************
% frequency corresponding to FFT bin
% ****************************************************************
function f = binFreq(n)
f = (mod(((0:n-1)+floor(n/2)), n)-floor(n/2))/n;
end
% ****************************************************************
% delay by phase shift
% needed only for demo code
% ****************************************************************
function waveform = cs_delay(waveform, rate_Hz, delay_s)
rflag = isreal(waveform);
n = numel(waveform);
cycLen_s = n / rate_Hz;
nCyc = delay_s / cycLen_s();
f = 0:(n - 1);
f = f + floor(n / 2);
f = mod(f, n);
f = f - floor(n / 2);
phase = -2 * pi * f * nCyc;
rot = exp(1i*phase);
waveform = ifft(fft(waveform) .* rot);
if rflag
waveform = real(waveform);
end
end
Circular Convolution
%circular convolution
%for testing you may use:
h = fir1(20,.3);
x = randn(1,1024);
%function y = conv_circ(h,x)
y = conv(h,x);
L1 = length(h);
L2 = length(x);
%add end to start, add start to end
temp = y(1:L1-1);
y(1:L1-1) = y(1:L1-1) + y(L2+(1:L1-1));
y(L2+(1:L1-1)) = y(L2+(1:L1-1)) + temp;
%compare to direct convolution
y2 = conv(h,x);
plot(y,'o-');hold
plot(y2,'r.--')
legend('circular','direct')
Power Computation of a Digital Stream
clear all;
%example vector having ac & dc power
x = complex(randn(1,2048)+.2,randn(1,2048)-.1);
%total power, three equivalent methods
pwr1_total = mean(abs(x).^2); %mean of squared values
pwr2_total = mean(real(x).^2 + imag(x).^2);
pwr3_total = mean(x .* conj(x));
%total expressed as rms
rms_total = sqrt(pwr1_total);
%dc power
pwr1_dc = mean(real(x))^2 + mean(imag(x))^2; %square of mean of values
%ac power
pwr1_ac = pwr1_total - pwr1_dc; %mean of squares - square of mean
%relation of ac power to statistical variance
pwr2_ac = var(x); %approximately
%ac expressed as rms
rms1_ac = sqrt(pwr1_ac);
%ac relation to standard of deviation, std = sqrt(var)
rms2_ac = std(x); %approximately
%dc relation to variance
pwr2_dc = pwr1_total - var(x); %approximately
fprintf('----------------------------------------------------\r');
fprintf('power(total), : %1.5f, %1.5f, %1.5f\r',pwr1_total,pwr2_total,pwr3_total);
fprintf('rms(total) : %1.5f\r',rms_total);
fprintf('power(ac), variance : %1.5f, %1.5f\r',pwr1_ac,pwr2_ac);
fprintf('rms(ac), std : %1.5f, %1.5f\r',rms1_ac,rms2_ac);
fprintf('power(dc), (total-var) : %1.5f, %1.5f\r',pwr1_dc,pwr2_dc);
Computing CIC Filter Register Pruning Using Matlab
% Filename: CIC_Word_Truncation.m
%
% Computes CIC decimation filter accumulator register
% truncation in each filter stage based on Hogenauer's
% 'accumulator register pruning' technique.
%
% Inputs:
% N = number of decimation CIC filter stages (filter order).
% R = CIC filter rate change factor (decimation factor).
% M = differential delay.
% Bin = number of bits in an input data word.
% Bout = number of bits in the filter's final output data word.
% Outputs:
% Stage number (ranges from 1 -to- 2*N+1).
% Bj = number of least significant bits that can be truncated
% at the input of a filter stage.
% Accumulator widths = number of a stage's necessary accumulator
% bits accounting for truncation.
% Richard Lyons Feb., 2012
clear, clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define CIC filter parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%N = 4; R = 25; M = 1; Bin = 16; Bout = 16; % Hogenauer paper, pp. 159
%N = 3; R = 32; M = 2; Bin = 8; Bout = 10; % Meyer Baese book, pp. 268
%N = 3; R = 16; M = 1; Bin = 16; Bout = 16; % Thorwartl's PDF file
N = 3; R = 8; M = 1; Bin = 12; Bout = 12; % Lyons' blog Figure 2 example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find h_sub_j and "F_sub_j" values for (N-1) cascaded integrators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
for j = N-1:-1:1
h_sub_j = [];
h_sub_j((R*M-1)*N + j -1 + 1) = 0;
for k = 0:(R*M-1)*N + j -1
for L = 0:floor(k/(R*M)) % Use uppercase "L" for loop variable
Change_to_Result = (-1)^L*nchoosek(N, L)...
*nchoosek(N-j+k-R*M*L,k-R*M*L);
h_sub_j(k+1) = h_sub_j(k+1) + Change_to_Result;
end % End "L" loop
end % End "k" loop
F_sub_j(j) = sqrt(sum(h_sub_j.^2));
end % End "j" loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define "F_sub_j" values for up to seven cascaded combs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j_for_many_combs = sqrt([2, 6, 20, 70, 252, 924, 3432]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute F_sub_j for last integrator stage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j(N) = F_sub_j_for_many_combs(N-1)*sqrt(R*M); % Last integrator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute F_sub_j for N cascaded filter's comb stages
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2*N:-1:N+1
F_sub_j(j) = F_sub_j_for_many_combs(2*N -j + 1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define "F_sub_j" values for the final output register truncation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j(2*N+1) = 1; % Final output register truncation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute column vector of minus log base 2 of "F_sub_j" values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Minus_log2_of_F_sub_j = -log2(F_sub_j)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute total "Output_Truncation_Noise_Variance" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CIC_Filter_Gain = (R*M)^N;
Num_of_Bits_Growth = ceil(log2(CIC_Filter_Gain));
% The following is from Hogenauer's Eq. (11)
%Num_Output_Bits_With_No_Truncation = Num_of_Bits_Growth +Bin -1;
Num_Output_Bits_With_No_Truncation = Num_of_Bits_Growth +Bin;
Num_of_Output_Bits_Truncated = Num_Output_Bits_With_No_Truncation -Bout;
Output_Truncation_Noise_Variance = (2^Num_of_Output_Bits_Truncated)^2/12;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute log base 2 of "Output_Truncation_Noise_Standard_Deviation" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Output_Truncation_Noise_Standard_Deviation = ...
sqrt(Output_Truncation_Noise_Variance);
Log_base2_of_Output_Truncation_Noise_Standard_Deviation = ...
log2(Output_Truncation_Noise_Standard_Deviation);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute column vector of "half log base 2 of 6/N" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Half_Log_Base2_of_6_over_N = 0.5*log2(6/N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute desired "B_sub_j" vector
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B_sub_j = floor(Minus_log2_of_F_sub_j ...
+ Log_base2_of_Output_Truncation_Noise_Standard_Deviation ...
+ Half_Log_Base2_of_6_over_N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' '), disp(' ')
disp(['N = ',num2str(N),', R = ',num2str(R),', M = ',num2str(M),...
', Bin = ', num2str(Bin),', Bout = ',num2str(Bout)])
disp(['Num of Bits Growth Due To CIC Filter Gain = ', ...
num2str(Num_of_Bits_Growth)])
disp(['Num of Accumulator Bits With No Truncation = ', ...
num2str(Num_Output_Bits_With_No_Truncation)])
% disp(['Output Truncation Noise Variance = ', ...
% num2str(Output_Truncation_Noise_Variance)])
% disp(['Log Base2 of Output Truncation Noise Standard Deviation = ',...
% num2str(Log_base2_of_Output_Truncation_Noise_Standard_Deviation)])
% disp(['Half Log Base2 of 6/N = ', num2str(Half_Log_Base2_of_6_over_N)])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create and display "Results" matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for Stage = 1:2*N
Results(Stage,1) = Stage;
Results(Stage,2) = F_sub_j(Stage);
Results(Stage,3) = Minus_log2_of_F_sub_j(Stage);
Results(Stage,4) = B_sub_j(Stage);
Results(Stage,5) = Num_Output_Bits_With_No_Truncation -B_sub_j(Stage);
end
% Include final output stage truncation in "Results" matrix
Results(2*N+1,1) = 2*N+1; % Output stage number
Results(2*N+1,2) = 1;
Results(2*N+1,4) = Num_of_Output_Bits_Truncated;
Results(2*N+1,5) = Bout;
%Results % Display "Results" matrix in raw float-pt.form
% % Display "F_sub_j" values if you wish
% disp(' ')
% disp(' Stage Fj -log2(Fj) Bj Accum width')
% for Stage = 1:2*N+1
% disp([' ',sprintf('%2.2g',Results(Stage,1)),sprintf('\t'),sprintf('%12.3g',Results(Stage,2)),...
% sprintf('\t'),sprintf('%7.5g',Results(Stage,3)),sprintf('\t'),...
% sprintf('%7.5g',Results(Stage,4)),sprintf('\t'),sprintf('%7.5g',Results(Stage,5))])
% end
% Display Stage number, # of truncated input bits, & Accumulator word widths
disp(' ')
disp(' Stage(j) Bj Accum (adder) width')
for Stage = 1:2*N
disp([' ',sprintf('%2.0g',Results(Stage,1)),...
sprintf('\t'),...
sprintf('%5.5g',Results(Stage,4)),sprintf('\t'),...
sprintf('%7.5g',Results(Stage,5))])
end
disp([' ',sprintf('%2.0g',Results(2*N+1,1)),...
sprintf('\t'),...
sprintf('%5.5g',Results(2*N+1,4)),sprintf('\t'),...
sprintf('%7.5g',Results(2*N+1,5)),' (final truncation)'])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Reading text files in Matlab
clear all;
% create test file
data = [1 5 2 4 3 3 4 2 5 1];
filename = 'test_file.txt';
fid = fopen(filename, 'w');
fprintf(fid, '%d %d\n', data);
fclose(fid);
% (1) read file using fscanf
fid = fopen(filename, 'r');
y1 = fscanf(fid, '%d\n'); %interleaves columns
fclose(fid);
% (2) read file using textread (or textscan)
[ya,yb] = textread(filename,'%d%d');
y2 = [ya yb];
% (3) read file using importdata
y3 = importdata(filename);
% (4) read file using load
y4 = load(filename);
disp('-------------------------------------------------------------')
disp(' original vector data')
disp(data)
disp(' file content using fprintf')
disp(y2)
disp(' vector created by fscanf')
disp(y1)
disp(' matrix created by:')
disp(' textread importdata load')
disp([y2 y3 y4])
Ideal interpolation filter
clear all; close all;
%example bandlimited random input & parameters
x = filter(fir1(70,.1),1,randn(1,1024));
up = 3; %Interpolation factor
cutoff = .3;
intorder = 6;
h1 = intfilt(up, intorder, cutoff); %ideal filter
h1 = up*h1/sum(h1);
h2 = fir1(2*up*intorder-2,cutoff); %ordinary LPF
h2 = up*h2/sum(h2);
%upsample
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;
x_f1 = filter(h1,1,x_up);
x_f2 = filter(h2,1,x_up);
figure;
subplot(3,1,1);hold
plot(x_f1,'o--');
plot(x_f2,'r.-');
legend('ideal output','fir1 output');
subplot(3,1,2);
plot(x_f1-x_f2);
legend('error');
subplot(3,1,3);hold
plot(h1,'.-');
plot(h2,'r.-');
legend('ideal filter','fir1 filter');
checking resampling in time domain
clear all; close all;
%example bandlimited random input & parameters
x = filter(fir1(70,.1),1,randn(1,512));
h = fir1(30,.3); %filter used for upsampling
up = 3; %Interpolation factor
dn = 2; %Decimation factor
%%%%%%%%%%%%%%%%%%%%% up/dn model %%%%%%%%%%%%%%%%%%%%%%
%upsample input
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;
x_f = filter(up*h,1,x_up);
%downsample signal by decimation
x_dn = x_f(1:dn:end);
delay = 30/2+1;
figure;
subplot(2,1,1);hold;
plot(x_up,'o--');
plot(x_f(delay:end),'r.-');
legend('input with zeros','upsampled stage');
subplot(2,1,2);hold;
plot(x_up(1:dn:end),'o--');
plot(x_dn(ceil(delay/dn):end),'r.-');
legend('input with zeros','final signal');
Resampling filter performance
%%%%%%%%%%%%%%%%%% inputs for model %%%%%%%%%%%%%%%%
clear all; close all;
%example bandlimited impulse input & parameters
x = filter(fir1(70,.1),1,[1 zeros(1,2^15-1)]);
Fs = 120; %MHz original sample rate
h = fir1(30,.3); %filter used for upsampling
up = 3; %Interpolation factor
dn = 2; %Decimation factor
Fc = 12; %MHz band centre (-Fs/2 ~ +Fs/2)
Fch = 0; %MHz filter centre (-Fs*up/2 ~ +Fs*up/2)
%move signal to its centre
x = x.*exp(j*2*pi*(0:length(x)-1)*Fc/Fs);
%shift filter
h = h.*exp(j*2*pi*(0:length(h)-1)*Fch/(Fs*up));
%%%%%%%%%%%%%%%%%%%%% model %%%%%%%%%%%%%%%%%%%%%%
%check signal in upsampled domain
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;
[P, F] = pwelch(x_up, [], 0, 2^16, Fs*up,'twosided');
F = F - max(F)/2;
P = fftshift(P);
y(find(P == 0)) = -100; %avoid log of zero
y(find(P ~= 0)) = 10*log10(P(find(P ~= 0)));
P_dB = y - 10*log10(max(P)); %normalise
%check filter response in upsampled domain
H = fftshift(20*log10(abs(fft(h,2^16))));
subplot(2,1,1);
hold;grid;
plot(F, P_dB,'.-');
plot(F,H,'m--');
axis([min(F)-1 max(F)+1 -80 1]);
legend('upsampled signal','upsampling filter');
%check signal in downsampled domain
x_f = filter(h,1,x_up);
x_dn = x_f(1:dn:end);
[P, F] = pwelch(x_dn, [], 0, 2^16, Fs*up/dn,'twosided');
F = F - max(F)/2;
P = fftshift(P);
y(find(P == 0)) = -100; %avoid log of zero
y(find(P ~= 0)) = 10*log10(P(find(P ~= 0)));
P_dB = y - 10*log10(max(P)); %normalise
subplot(2,1,2)
plot(F,P_dB,'r.-');
grid;
axis([min(F)-1 max(F)+1 -80 1]);
legend('downsampled signal');