% Filename: FFT_Twiddles_Find_DSPrelated.m
% Computes 'Decimation in Frequency' or 'Decimation
% in Time' Butterfly twiddle factors, for radix-2 FFTs
% with in-order input indices and scrambled output indices.
%
% To use, do two things: (1) define FFT size 'N'; and
% (2) define the desired 'Structure', near line 17-18,
% as 'Dec_in_Time' or 'Dec_in_Freq'.
%
% Author: Richard Lyons, November, 2011
%******************************************
clear, clc
% Define input parameters
N = 8; % FFT size (Must be an integer power of 2)
Structure = 'Dec_in_Time'; % Choose Dec-in-time butterflies
%Structure = 'Dec_in_Freq'; % Choose Dec-in-frequency butterflies
% Start of processing
Num_Stages = log2(N); % Number of stages
StageStart = 1; % First stage to compute
StageStop = Num_Stages; % Last stage to compute
ButterStart = 1; %First butterfly to compute
ButterStop = N/2; %Last butterfly to compute
Pointer = 0; %Init 'results' row pointer
for Stage_Num = StageStart:StageStop
if Structure == 'Dec_in_Time'
for Butter_Num = ButterStart:ButterStop
Twid = floor((2^Stage_Num*(Butter_Num-1))/N);
% Compute bit reversal of Twid
Twid_Bit_Rev = 0;
for I = Num_Stages-2:-1:0
if Twid>=2^I
Twid_Bit_Rev = Twid_Bit_Rev + 2^(Num_Stages-I-2);
Twid = Twid -2^I;
else, end
end %End bit reversal 'I' loop
A1 = Twid_Bit_Rev; %Angle A1
A2 = Twid_Bit_Rev + N/2; %Angle A2
Pointer = Pointer +1;
Results(Pointer,:) = [Stage_Num,Butter_Num,A1,A2];
end
else
for Twiddle_Num = 1:N/2^Stage_Num
Twid = (2^Stage_Num*(Twiddle_Num-1))/2; %Compute integer
Pointer = Pointer +1;
Results(Pointer,:) = [Stage_Num,Twiddle_Num,Twid];
end
end % End 'if'
end % End Stage_Num loop
Results(:,1:3), disp(' Stage# Twid# A'), disp(' ')
% WCDMA channelization codes
% source:
% 3GPP TS 25.213 V10.0.0 section 4.3.1.1
%
% parameters:
% sf: spreading factor
% k: code number
%
% The returned code is a column vector with length sf
%
% this code has not been tested extensively. Please verify
% independently that it does what it promises.
function code=UTRA_FDD_chanCode(sf, k)
persistent flag;
persistent codes;
% * ********************************************
% * Build codebook
% * ********************************************
if isempty(flag)
codes={};
flag=1;
% start of recursion: Identity code for sf=1
codes{1, 1}=1;
for sfi=1:8
sfg=2 ^ sfi;
for kgDest=0:2:sfg-2
kgSrc=kgDest/2;
prev=codes{sfg/2, kgSrc+1};
% generation method:
% - duplicate a lower-sf code
% - duplicate and change sign
codes{sfg, kgDest+1}=[prev prev];
codes{sfg, kgDest+2}=[prev -prev];
end
end
end
% * ********************************************
% * look up the requested code from codebook
% * ********************************************
code=transpose(codes{sf, k+1});
% ################## CUT HERE #########################
% Example use (put this into a separate file)
sf=128; codenum=3;
chanCode=UTRA_FDD_chanCode(sf, codenum);
sig=[1 0 0 1 0 0 ]; % signal, row vector
len=size(sig, 2),
% convolve:
s=chanCode * sig; % now matrix, one row per repetition
len=len*sf;
s=reshape(s, 1, len);
% plot
stem(s);
% **************************************************************
% Vectorized Farrow resampler
% M. Nentwig, 2011
% Note: Uses cyclic signals (wraps around)
% **************************************************************
close all; clear all;
% inData contains input to the resampling process (instead of function arguments)
inData = struct();
% **************************************************************
% example coefficients.
% Each column [c0; c1; c2; c3] describes a polynomial for one tap coefficent in fractional time ft [0, 1]:
% tapCoeff = c0 + c1 * ft + c2 * ft ^ 2 + c3 * ft ^ 3
% Each column corresponds to one tap.
% the matrix size may be changed arbitrarily.
%
% The example filter is based on a 6th order Chebyshev Laplace-domain prototype.
% **************************************************************
inData.cMatrix =[ -8.57738278e-3 7.82989032e-1 7.19303539e+000 6.90955718e+000 -2.62377450e+000 -6.85327127e-1 1.44681608e+000 -8.79147907e-1 7.82633997e-2 1.91318985e-1 -1.88573400e-1 6.91790782e-2 3.07723786e-3 -6.74800912e-3
2.32448021e-1 2.52624309e+000 7.67543936e+000 -8.83951796e+000 -5.49838636e+000 6.07298348e+000 -2.16053205e+000 -7.59142947e-1 1.41269409e+000 -8.17735712e-1 1.98119464e-1 9.15904145e-2 -9.18092030e-2 2.74136108e-2
-1.14183319e+000 6.86126458e+000 -6.86015957e+000 -6.35135894e+000 1.10745051e+001 -3.34847578e+000 -2.22405694e+000 3.14374725e+000 -1.68249886e+000 2.54083065e-1 3.22275037e-1 -3.04794927e-1 1.29393976e-1 -3.32026332e-2
1.67363115e+000 -2.93090391e+000 -1.13549165e+000 5.65274939e+000 -3.60291782e+000 -6.20715544e-1 2.06619782e+000 -1.42159644e+000 3.75075865e-1 1.88433333e-1 -2.64135123e-1 1.47117661e-1 -4.71871047e-2 1.24921920e-2
] / 12.28;
% **************************************************************
% Create example signal
% **************************************************************
nIn = 50; % used for test signal generation only
if true
% complex signal
inData.signal = cos(2*pi*(0:(nIn-1)) / nIn) + cos(2*2*pi*(0:(nIn-1)) / nIn) + 1.5;
else
% impulse response
inData.signal = zeros(1, nIn); inData.signal(1) = 1; %
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(1, 1) = 1; % enable to show constant c in first tap
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(2, 1) = 1; % enable to show linear c in first tap
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(3, 1) = 1; % enable to show quadratic c in first tap
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(1, 2) = 1; % enable to show constant c in second tap
end
% **************************************************************
% Resample to the following number of output samples
% must be integer, otherwise arbitrary
% **************************************************************
inData.nSamplesOut = floor(nIn * 6.28);
% **************************************************************
% Set up Farrow resampling
% **************************************************************
nSamplesIn = size(inData.signal, 2);
nSamplesOut = inData.nSamplesOut;
order = size(inData.cMatrix, 1) - 1; % polynomial order
nTaps = size(inData.cMatrix, 2); % number of input samples that contribute to one output sample (FIR size)
% pointer to the position in the input stream for each output sample (row vector, real numbers), starting at 0
inputIndex = (0:nSamplesOut-1) / nSamplesOut * nSamplesIn;
% split into integer part (0..nSamplesIn - 1) ...
inputIndexIntegerPart = floor(inputIndex);
% ... and fractional part [0, 1[
inputIndexFractionalPart = inputIndex - inputIndexIntegerPart;
% **************************************************************
% Calculate output stream
% First constant term (conventional FIR), then linear, quadratic, cubic, ...
% **************************************************************
outStream = zeros(1, inData.nSamplesOut);
for ixOrder = 0 : order
x = inputIndexFractionalPart .^ ixOrder;
% **************************************************************
% Add the contribution of each tap one-by-one
% **************************************************************
for ixTap = 0 : nTaps - 1
c = inData.cMatrix(ixOrder+1, ixTap+1);
% index of input sample that contributes to output via the current tap
% higher tap index => longer delay => older input sample => smaller data index
dataIx = inputIndexIntegerPart - ixTap;
% wrap around
dataIx = mod(dataIx, nSamplesIn);
% array indexing starts at 1
dataIx = dataIx + 1;
delayed = inData.signal(dataIx);
% for each individual output sample (index in row vector),
% - evaluate f = c(order, tapindex) * fracPart .^ order
% - scale the delayed input sample with f
% - accumulate the contribution of all taps
% this implementation performs the operation for all output samples in parallel (row vector)
outStream = outStream + c * delayed .* x;
end % for ixTap
end % for ixOrder
% **************************************************************
% plot
% **************************************************************
xIn = linspace(0, 1, nSamplesIn + 1); xIn = xIn(1:end-1);
xOut = linspace(0, 1, nSamplesOut + 1); xOut = xOut(1:end-1);
figure(); grid on; hold on;
stem(xIn, inData.signal, 'k+-');
plot(xOut, outStream, 'b+-');
legend('input', 'output');
title('Farrow resampling. Signals are cyclic');
function [Windowed_Spec] = Wind_Flattop(Spec)
% Given an input spectral sequence 'Spec', that is the
% FFT of some time sequence 'x', Wind_Flattop(Spec)
% returns a spectral sequence that is equivalent
% to the FFT of a flat-top windowed version of time
% sequence 'x'. The peak magnitude values of output
% sequence 'Windowed_Spec' can be used to accurately
% estimate the peak amplitudes of sinusoidal components
% in time sequence 'x'.
% Input: 'Spec' (a sequence of complex FFT sample values)
%
% Output: 'Windowed_Spec' (a sequence of complex flat-top
% windowed FFT sample values)
%
% Based on Lyons': "Reducing FFT Scalloping Loss Errors
% Without Multiplication", IEEE Signal Processing Magazine,
% DSP Tips & Tricks column, March, 2011, pp. 112-116.
%
% Richard Lyons [December, 2011]
N = length(Spec);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform freq-domain convolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g_Coeffs = [1, -0.94247, 0.44247];
% Compute first two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(1) = g_Coeffs(3)*Spec(N-1) ...
+g_Coeffs(2)*Spec(N) + Spec(1) ...
+g_Coeffs(2)*Spec(2) + g_Coeffs(3)*Spec(3);
Windowed_Spec(2) = g_Coeffs(3)*Spec(N) ...
+g_Coeffs(2)*Spec(1) + Spec(2) ...
+g_Coeffs(2)*Spec(3) + g_Coeffs(3)*Spec(4);
% Compute last two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(N-1) = g_Coeffs(3)*Spec(N-3) ...
+g_Coeffs(2)*Spec(N-2) + Spec(N-1) ...
+g_Coeffs(2)*Spec(N) + g_Coeffs(3)*Spec(1);
Windowed_Spec(N) = g_Coeffs(3)*Spec(N-2) ...
+g_Coeffs(2)*Spec(N-1) + Spec(N) ...
+g_Coeffs(2)*Spec(1) + g_Coeffs(3)*Spec(2);
% Compute convolved spec samples for the middle of the spectrum
for K = 3:N-2
Windowed_Spec(K) = g_Coeffs(3)*Spec(K-2) ...
+g_Coeffs(2)*Spec(K-1) + Spec(K) ...
+g_Coeffs(2)*Spec(K+1) + g_Coeffs(3)*Spec(K+2);
end % % End of 'Wind_Flattop(Spec)' function
% Lagrange interpolation for resampling
% References:
% [1] A digital signal processing approach to Interpolation
% Ronald W. Schafer and Lawrence R. Rabiner
% Proc. IEEE vol 61, pp.692-702, June 1973
% [2] https://ccrma.stanford.edu/~jos/Interpolation/Lagrange_Interpolation.html
% [3] Splitting the Unit delay: Tools for fractional delay filter design
% T. I. Laakso, V. Valimaki, M. Karjalainen, and U. K. Laine
% IEEE Signal Processing Magazine, vol. 13, no. 1, pp. 30..60, January 1996
function lagrangeResamplingDemo()
originDefinition = 0; % see comment for LagrangeBasisPolynomial() below
% Regarding order: From [1]
% "Indeed, it is easy to show that whenever Q is odd, none of the
% impulse responses corresponding to Lagrange interpolation can have
% linear phase"
% Here, order = Q+1 => odd orders are preferable
order = 3;
% *******************************************************************
% Set up signals
% *******************************************************************
nIn = order + 1;
nOut = nIn * 5 * 20;
% Reference data: from [1] fig. 8, linear-phase type
ref = [-0.032, -0.056, -0.064, -0.048, 0, ...
0.216, 0.448, 0.672, 0.864, 1, ...
0.864, 0.672, 0.448, 0.216, 0, ...
-0.048, -0.064, -0.056, -0.032, 0];
tRef = (1:size(ref, 2)) / size(ref, 2);
% impulse input to obtain impulse response
inData = zeros(1, nIn);
inData(1) = 1;
outData = zeros(1, nOut);
outTime = 0:(nOut-1);
outTimeAtInput = outTime / nOut * nIn;
outTimeAtInputInteger = floor(outTimeAtInput);
outTimeAtInputFractional = outTimeAtInput - outTimeAtInputInteger;
evalFracTime = outTimeAtInputFractional - 0.5 + originDefinition;
% odd-order modification
if mod(order, 2) == 0
% Continuity of the impulse response is achieved when support points are located at
% the intersections between adjacent segments "at +/- 0.5"
% For an even order polynomial (odd number of points), this is only possible with
% an asymmetric impulse response
offset = 0.5;
%offset = -0.5; % alternatively, its mirror image
else
offset = 0;
end
% *******************************************************************
% Apply Lagrange interpolation to input data
% *******************************************************************
for ixTap = 0:order
% ixInSample is the input sample that appears at FIR tap 'ixTap' to contribute
% to the output sample
% Row vector, for all output samples in parallel
ixInSample = mod(outTimeAtInputInteger + ixTap - order, nIn) + 1;
% the value of said input sample, for all output samples in parallel
d = inData(ixInSample);
% Get Lagrange polynomial coefficients of this tap
c = lagrangeBasisPolynomial(order, ixTap, originDefinition + offset);
% Evaluate the Lagrange polynomial, resulting in the time-varying FIR tap weight
cTap = polyval(c, evalFracTime);
% FIR operation: multiply FIR tap weight with input sample and add to
% output sample (all outputs in parallel)
outData = outData + d .* cTap;
end
% *******************************************************************
% Plot
% *******************************************************************
figure(); hold on;
h = plot((0:nOut-1) / nOut, outData, 'b-'); set(h, 'lineWidth', 3);
stem(tRef, ref, 'r'); set(h, 'lineWidth', 3);
legend('impulse response', 'reference result');
title('Lagrange interpolation, impulse response');
end
% Returns the coefficients of a Lagrange basis polynomial
% 1 <= order: Polynomial order
% 0 <= evalIx <= order: index of basis function.
%
% At the set of support points, the basis polynomials evaluate as follows:
% evalIx = 1: [1 0 0 ...]
% evalIx = 2: [0 1 0 ...]
% evalIx = 3: [0 0 1 ...]
%
% The support point are equally spaced.
% Example, using originDefinition=0:
% order = 1: [-0.5 0.5]
% order = 2: [-1 0 1]
% order = 3: [-1.5 -0.5 0.5 1.5]
%
% The area around the mid-point corresponds to -0.5 <= x <= 0.5.
% If a resampler implementation uses by convention 0 <= x <= 1 instead, set
% originDefinition=0.5 to offset
% the polynomial.
function polyCoeff = lagrangeBasisPolynomial(order, evalIx, originDefinition)
assert(evalIx >= 0);
assert(evalIx <= order);
tapLocations = -0.5*(order) + (0:order) + originDefinition;
polyCoeff = [1];
for loopIx = 0:order
if loopIx ~= evalIx
% numerator: places a zero in the polynomial via (x-xTap(k)), with k != evalIx
% denominator: scales to 1 amplitude at x=xTap(evalIx)
term = [1 -tapLocations(loopIx+1)] / (tapLocations(evalIx+1)-tapLocations(loopIx+1));
% multiply polynomials => convolve coefficients
polyCoeff = conv(polyCoeff, term);
end
end
% TEST:
% The Lagrange polynomial evaluates to 1 at the location of the tap
% corresponding to evalIx
thisTapLocation = tapLocations(evalIx+1);
pEval = polyval(polyCoeff, thisTapLocation);
assert(max(abs(pEval) - 1) < 1e6*eps);
% The Lagrange polynomial evaluates to 0 at all other tap locations
tapLocations(evalIx+1) = [];
pEval = polyval(polyCoeff, tapLocations);
assert(max(abs(pEval)) < 1e6*eps);
end
%For a real-valued, single tone sine or cosine wave, estimate the frequency
%using an fft. Resolution will be limited by the sampling rate of the tone
%and the number of points in the fft.
%Procedure:
%1. Generate the test tone with a given Fs and N
%2. Take the fft and keep the first half
%3. Detect the maximum peak and its index
%4. Equate this index to a frequency
%Error Estimate:
% (Fs/N)= Hz/bin
%Things to keep in mind:
%Adding noise to the input will skew the results
%windowing will help with the fft
%% Clear and close everything
clc;
close all;
clear all; %remove this line if you are trying to use breakpoints
%% Just copy into m file, save and run
Fs = 1e6; %1MHz
fi = 1.333e3; %1kHz
t = 0:1/Fs:.5;
y = sin(2*pi*fi*t);
%Plot the time and frequency domain. Be sure to zoom in to see the waveform
%and spectrum
figure;
subplot(2,1,1);
temp = (2/fi)*Fs;
plot(y);
xlabel('time (s)');
subplot(2,1,2);
sX=fft(y);
N=length(sX);
sXM = abs(sX(1:floor(end/2))).^2; %take the magnitude and only keep 0:Fs/2
plot(0:Fs/N:(Fs/2)-1,sXM);
xlabel('Frequency')
ylabel('Magnitude')
[vv, ii]=max(sXM); %find the index of the largest value in the spectrum
freqEst = (ii-1)*Fs/N; %units are Hz
resMin = Fs/N; %units are Hz
display(freqEst);%display to the command window the results
display(resMin);
%Normalized Least Mean Square Adaptive Filter
%for Echo Cancellation
function [e,h,t] = adapt_filt_nlms(x,B,h,delta,l,l1)
%x = the signal from the speaker
%B = signal through echo path
%h = impulse response of adaptive filter
%l = length of the signal x
%l1 = length of the adaptive filter
%delta = step size
%e = error
%h = adaptive filter impulse response
tic;
for k = 1:150
for n= 1:l
xcap = x((n+l1-1):-1:(n+l1-1)-l1+1);
yout(n) = h*xcap';
e(n) = B(n) - yout(n);
xnorm = 0.000001+(xcap*xcap');
h = h+((delta*e(n))*(xcap/xnorm));
end
eold = 0.0;
for i = 1:l
mesquare = eold+(e(i)^2);
eold = mesquare;
end
if mesquare <=0.0001
break;
end
end
t = toc; %to get the time elapsed
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 [Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB)
%
% SNR_Set(x, Desired_SNR_dB) returns the real-valued
% input 'Signal' contaminated with normally-distributed,
% zero-mean, random noise. The signal-to-noise ratio
% (SNR in dB) of the output 'Noisy_Signal' signal is
% controlled by the input 'Desired_SNR_dB' argument measured
% in dB.
% Example:
% Npts = 128; % Number of time samples
% n = 0:Npts-1; % Time-domain index
% Signal = 3*sin(2*pi*3*n/Npts); % Real-valued signal
% Desired_SNR_dB = 3; % Set SNR of output 'Noisy_Signal' to +3 dB
% [Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB);
%
% Author: Richard Lyons [December 2011]
%******************************************
Npts = length(Signal); % Number of input time samples
Noise = randn(1,Npts); % Generate initial noise; mean zero, variance one
Signal_Power = sum(abs(Signal).*abs(Signal))/Npts;
Noise_Power = sum(abs(Noise).*abs(Noise))/Npts;
%Initial_SNR = 10*(log10(Signal_Power/Noise_Power));
K = (Signal_Power/Noise_Power)*10^(-Desired_SNR_dB/10); % Scale factor
New_Noise = sqrt(K)*Noise; % Change Noise level
%New_Noise_Power = sum(abs(New_Noise).*abs(New_Noise))/Npts
%New_SNR = 10*(log10(Signal_Power/New_Noise_Power))
Noisy_Signal = Signal + New_Noise;
'SNR_Set()' Function Test Code:
% Filename SNR_Set_test.m
%
% Tests the 'SNR_Set()" function. Adds a predefined
% amount of random noise to a noise-free signal such that
% the noisy signal has a desired signal-to-noise ratio (SNR).
%
% Author: Richard Lyons [December 2011]
clear, clc
% Create a noise-free signal
Npts = 128; % Number of time samples
n = 0:Npts-1; % Time-domain index
Cycles = 5; % Integer number of cycles in noise-free sinwave signal
Signal = 3*sin(2*pi*Cycles*n/Npts); % Real-valued noise-free signal
Desired_SNR_dB = 3 % Set desired SNR in dB
[Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB); % Generate noisy signal
% Plot original and 'noisy' signals
figure(1)
subplot(2,1,1)
plot(n, Signal, '-bo', 'markersize', 2)
title('Original Signal')
grid on, zoom on
subplot(2,1,2)
plot(n, Noisy_Signal, '-bo', 'markersize', 2)
title('Noisy Signal')
xlabel('Time-samples')
grid on, zoom on
% Measure SNR in freq domain
Spec = fft(Noisy_Signal);
Spec_Mag = abs(Spec); % Spectral magnitude
figure(2)
plot(Spec_Mag, '-bo', 'markersize', 2)
title('Spec Mag of Noisy Signal')
xlabel('Freq-samples'), ylabel('Linear')
grid on, zoom on
Signal_Power = Spec_Mag(Cycles+1)^2 + Spec_Mag(Npts-Cycles+1)^2;
Noise_Power = sum(Spec_Mag.^2) -Signal_Power;
Measured_SNR = 10*log10(Signal_Power/Noise_Power)
% *************************************************
% alias- and in-channel error analysis for root-raised
% cosine filter with upsampler FIR cascade
% Markus Nentwig, 25.12.2011
%
% * plots the aliases at the output of each stage
% * plots the error spectrum (deviation from ideal RRC-
% response)
% *************************************************
function eval_RRC_resampler()
1;
% variant 1
% conventional RRC filter and FIR resampler
smode = 'evalConventional';
% export resampler frequency response to design equalizing RRC filter
%smode = 'evalIdeal';
% variant 2
% equalizing RRC filter and FIR resampler
% smode = 'evalEqualized';
% *************************************************
% load impulse responses
% *************************************************
switch smode
case 'evalConventional'
% conventionally designed RRC filter
h0 = load('RRC.dat');
case 'evalIdeal'
h0 = 1;
case 'evalEqualized'
% alternative RRC design that equalizes the known frequency
% response of the resampler
h0 = load('RRC_equalized.dat');
otherwise assert(false);
end
h1 = load('ip1.dat');
h2 = load('ip2.dat');
h3 = load('ip3.dat');
h4 = load('ip4.dat');
% *************************************************
% --- signal source ---
% *************************************************
n = 10000; % test signal, number of symbol lengths
rate = 1;
s = zeros(1, n);
s(1) = 1;
p = {};
p = addPlot(p, s, rate, 'k', 5, 'sym stream r=1');
% *************************************************
% --- upsampling RRC filter ---
% *************************************************
rate = rate * 2;
s = upsample(s, 2); % insert one zero after every sample
p = addPlot(p, s, rate, 'k', 2, 'sym stream r=2');
s = filter(h0, [1], s);
p = addPlot(p, s, rate, 'b', 3, 'sym stream r=2, filtered');
p = addErrPlot(p, s, rate, 'error');
figure(1); clf; grid on; hold on;
doplot(p, 'interpolating pulse shaping filter');
ylim([-70, 2]);
p = {};
% *************************************************
% --- first interpolator by 2 ---
% *************************************************
rate = rate * 2;
s = upsample(s, 2); % insert one zero after every sample
p = addPlot(p, s, rate, 'k', 3, 'interpolator 1 input');
s = filter(h1, [1], s);
p = addPlot(p, s, rate, 'b', 3, 'interpolator 1 output');
p = addErrPlot(p, s, rate, 'error');
figure(2); clf; grid on; hold on;
doplot(p, 'first interpolator by 2');
ylim([-70, 2]);
p = {};
% *************************************************
% --- second interpolator by 2 ---
% *************************************************
rate = rate * 2;
s = upsample(s, 2); % insert one zero after every sample
p = addPlot(p, s, rate, 'k', 3, 'interpolator 2 input');
s = filter(h2, [1], s);
p = addPlot(p, s, rate, 'b', 3, 'interpolator 2 output');
p = addErrPlot(p, s, rate, 'error');
figure(3); clf; grid on; hold on;
doplot(p, 'second interpolator by 2');
ylim([-70, 2]);
p = {};
% *************************************************
% --- third interpolator by 2 ---
% *************************************************
rate = rate * 2;
s = upsample(s, 2); % insert one zero after every sample
p = addPlot(p, s, rate, 'k', 3, 'interpolator 3 input');
s = filter(h3, [1], s);
p = addPlot(p, s, rate, 'b', 3, 'interpolator 3 output');
p = addErrPlot(p, s, rate, 'error');
figure(4); clf; grid on; hold on;
doplot(p, 'third interpolator by 2');
ylim([-70, 2]);
p = {};
% *************************************************
% --- fourth interpolator by 4 ---
% *************************************************
rate = rate * 4;
s = upsample(s, 4); % insert three zeros after every sample
p = addPlot(p, s, rate, 'k', 3, 'interpolator 4 input');
s = filter(h4, [1], s);
p = addPlot(p, s, rate, 'b', 3, 'final output');
p = addErrPlot(p, s, rate, 'error at output');
figure(5); clf; grid on; hold on;
doplot(p, 'fourth interpolator by 4');
ylim([-70, 2]);
figure(334);
stem(real(s(1:10000)));
% *************************************************
% export resampler frequency response
% *************************************************
switch smode
case 'evalIdeal'
exportFrequencyResponse(s, rate, 'interpolatorFrequencyResponse.mat');
end
end
% ************************************
% put frequency response plot data into p
% ************************************
function p = addPlot(p, s, rate, plotstyle, linewidth, legtext)
p{end+1} = struct('sig', s, 'rate', rate, 'plotstyle', plotstyle, 'linewidth', linewidth, 'legtext', legtext);
end
% ************************************
% determine the error spectrum, compared to an ideal filter (RRC)
% and add a plot to p
% ************************************
function p = addErrPlot(p, s, rate, legtext)
ref = RRC_impulseResponse(numel(s), rate);
% refB is scaled and shifted (sub-sample resolution) replica of ref
% that minimizes the energy in (s - refB)
[coeff, refB, deltaN] = fitSignal_FFT(s, ref);
err = s - refB;
err = brickwallFilter(err, rate, 1.15); % 1+alpha
% signal is divided into three parts:
% - A) wanted in-channel energy (correlated with ref)
% - B) unwanted in-channel energy (uncorrelated with ref)
% - C) unwanted out-of-channel energy (aliases)
% the error vector magnitude is B) relative to A)
energySig = refB * refB';
energyErr = err * err';
EVM_dB = 10*log10(energyErr / energySig);
legtext = sprintf('%s; EVM=%1.2f dB', legtext, EVM_dB);
p{end+1} = struct('sig', err, 'rate', rate, 'plotstyle', 'r', 'linewidth', 3, 'legtext', legtext);
end
% ************************************
% helper function, plot data in p
% ************************************
function doplot(p, t)
leg = {};
for ix = 1:numel(p)
pp = p{ix};
fb = FFT_frequencyBasis(numel(pp.sig), pp.rate);
fr = 20*log10(abs(fft(pp.sig) + eps));
h = plot(fftshift(fb), fftshift(fr), pp.plotstyle);
set(h, 'lineWidth', pp.linewidth);
xlabel('frequency, normalized to symbol rate');
ylabel('dB');
leg{end+1} = pp.legtext;
end
legend(leg);
title(t);
end
% ************************************
% ideal RRC filter (impulse response is as
% long as test signal)
% ************************************
function ir = RRC_impulseResponse(n, rate)
alpha = 0.15;
fb = FFT_frequencyBasis(n, rate);
% bandwidth is 1
c = abs(fb / 0.5);
c = (c-1)/(alpha); % -1..1 in the transition region
c=min(c, 1);
c=max(c, -1);
RRC_h = sqrt(1/2+cos(pi/2*(c+1))/2);
ir = real(ifft(RRC_h));
end
% ************************************
% remove any energy at frequencies > BW/2
% ************************************
function s = brickwallFilter(s, rate, BW)
n = numel(s);
fb = FFT_frequencyBasis(n, rate);
ix = find(abs(fb) > BW / 2);
s = fft(s);
s(ix) = 0;
s = real(ifft(s));
end
% ************************************
% for an impulse response s at rate, write the
% frequency response to fname
% ************************************
function exportFrequencyResponse(s, rate, fname)
fb = fftshift(FFT_frequencyBasis(numel(s), rate));
fr = fftshift(fft(s));
figure(335); grid on;
plot(fb, 20*log10(abs(fr)));
title('exported frequency response');
xlabel('normalized frequency');
ylabel('dB');
save(fname, 'fb', 'fr');
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
% *******************************************************
% delay-matching between two signals (complex/real-valued)
%
% * matches the continuous-time equivalent waveforms
% of the signal vectors (reconstruction at Nyquist limit =>
% ideal lowpass filter)
% * Signals are considered cyclic. Use arbitrary-length
% zero-padding to turn a one-shot signal into a cyclic one.
%
% * output:
% => coeff: complex scaling factor that scales 'ref' into 'signal'
% => delay 'deltaN' in units of samples (subsample resolution)
% apply both to minimize the least-square residual
% => 'shiftedRef': a shifted and scaled version of 'ref' that
% matches 'signal'
% => (signal - shiftedRef) gives the residual (vector error)
%
% Example application
% - with a full-duplex soundcard, transmit an arbitrary cyclic test signal 'ref'
% - record 'signal' at the same time
% - extract one arbitrary cycle
% - run fitSignal
% - deltaN gives the delay between both with subsample precision
% - 'shiftedRef' is the reference signal fractionally resampled
% and scaled to optimally match 'signal'
% - to resample 'signal' instead, exchange the input arguments
% *******************************************************
function [coeff, shiftedRef, deltaN] = fitSignal_FFT(signal, ref)
n=length(signal);
% xyz_FD: Frequency Domain
% xyz_TD: Time Domain
% all references to 'time' and 'frequency' are for illustration only
forceReal = isreal(signal) && isreal(ref);
% *******************************************************
% Calculate the frequency that corresponds to each FFT bin
% [-0.5..0.5[
% *******************************************************
binFreq=(mod(((0:n-1)+floor(n/2)), n)-floor(n/2))/n;
% *******************************************************
% Delay calculation starts:
% Convert to frequency domain...
% *******************************************************
sig_FD = fft(signal);
ref_FD = fft(ref, n);
% *******************************************************
% ... calculate crosscorrelation between
% signal and reference...
% *******************************************************
u=sig_FD .* conj(ref_FD);
if mod(n, 2) == 0
% for an even sized FFT the center bin represents a signal
% [-1 1 -1 1 ...] (subject to interpretation). It cannot be delayed.
% The frequency component is therefore excluded from the calculation.
u(length(u)/2+1)=0;
end
Xcor=abs(ifft(u));
% figure(); plot(abs(Xcor));
% *******************************************************
% Each bin in Xcor corresponds to a given delay in samples.
% The bin with the highest absolute value corresponds to
% the delay where maximum correlation occurs.
% *******************************************************
integerDelay = find(Xcor==max(Xcor));
% (1): in case there are several bitwise identical peaks, use the first one
% Minus one: Delay 0 appears in bin 1
integerDelay=integerDelay(1)-1;
% Fourier transform of a pulse shifted by one sample
rotN = exp(2i*pi*integerDelay .* binFreq);
uDelayPhase = -2*pi*binFreq;
% *******************************************************
% Since the signal was multiplied with the conjugate of the
% reference, the phase is rotated back to 0 degrees in case
% of no delay. Delay appears as linear increase in phase, but
% it has discontinuities.
% Use the known phase (with +/- 1/2 sample accuracy) to
% rotate back the phase. This removes the discontinuities.
% *******************************************************
% figure(); plot(angle(u)); title('phase before rotation');
u=u .* rotN;
% figure(); plot(angle(u)); title('phase after rotation');
% *******************************************************
% Obtain the delay using linear least mean squares fit
% The phase is weighted according to the amplitude.
% This suppresses the error caused by frequencies with
% little power, that may have radically different phase.
% *******************************************************
weight = abs(u);
constRotPhase = 1 .* weight;
uDelayPhase = uDelayPhase .* weight;
ang = angle(u) .* weight;
r = [constRotPhase; uDelayPhase] .' \ ang.'; %linear mean square
%rotPhase=r(1); % constant phase rotation, not used.
% the same will be obtained via the phase of 'coeff' further down
fractionalDelay=r(2);
% *******************************************************
% Finally, the total delay is the sum of integer part and
% fractional part.
% *******************************************************
deltaN = integerDelay + fractionalDelay;
% *******************************************************
% provide shifted and scaled 'ref' signal
% *******************************************************
% this is effectively time-convolution with a unit pulse shifted by deltaN
rotN = exp(-2i*pi*deltaN .* binFreq);
ref_FD = ref_FD .* rotN;
shiftedRef = ifft(ref_FD);
% *******************************************************
% Again, crosscorrelation with the now time-aligned signal
% *******************************************************
coeff=sum(signal .* conj(shiftedRef)) / sum(shiftedRef .* conj(shiftedRef));
shiftedRef=shiftedRef * coeff;
if forceReal
shiftedRef = real(shiftedRef);
end
end
% ****************************************************************
% baseband equivalent source of local oscillator with phase noise
% Markus Nentwig, 18.12.2011
% ****************************************************************
function pn_generator()
close all;
% ****************************************************************
% PN generator configuration
% ****************************************************************
srcPar = struct();
srcPar.n = 2 ^ 18; % generated number of output samples
srcPar.rate_Hz = 7.68e6; % sampling rate
srcPar.f_Hz = [0, 10e3, 1e6, 9e9]; % phase noise spectrum, frequencies
srcPar.g_dBc1Hz = [-80, -80, -140, -140]; % phase noise spectrum, magnitude
srcPar.spursF_Hz = [300e3, 400e3, 700e3]; % discrete spurs (set [] if not needed)
srcPar.spursG_dBc = [-50, -55, -60]; % discrete spurs, power relative to carrier
% ****************************************************************
% run PN generator
% ****************************************************************
s = PN_src(srcPar);
if false
% ****************************************************************
% export phase noise baseband-equivalent signal for use in simulator etc
% ****************************************************************
tmp = [real(s); imag(s)] .';
save('phaseNoiseSample.dat', 'tmp', '-ascii');
end
if exist('spectrumAnalyzer', 'file')
% ****************************************************************
% spectrum analyzer configuration
% ****************************************************************
SAparams = struct();
SAparams.rate_Hz = srcPar.rate_Hz; % sampling rate of the input signal
SAparams.pRef_W = 1e-3; % unity signal represents 0 dBm (1/1000 W)
SAparams.pNom_dBm = 0; % show 0 dBm as 0 dB;
SAparams.filter = 'brickwall';
SAparams.RBW_window_Hz = 1000; % convolve power spectrum with a 1k filter
SAparams.RBW_power_Hz = 1; % show power density as dBc in 1 Hz
SAparams.noisefloor_dB = -250; % don't add artificial noise
SAparams.logscale = true; % use logarithmic frequency axis
% plot nominal spectrum
figure(1); grid on;
h = semilogx(max(srcPar.f_Hz, 100) / 1e6, srcPar.g_dBc1Hz, 'k+-');
set(h, 'lineWidth', 3);
hold on;
spectrumAnalyzer('signal', s, SAparams, 'fMin_Hz', 100, 'fig', 1);
ylabel('dBc in 1 Hz');
legend('nominal PSD', 'output spectrum');
title('check match with nominal PSD');
spectrumAnalyzer('signal', s, SAparams, 'fMin_Hz', 100, 'RBW_power_Hz', 'sine', 'fig', 2);
title('check carrier level (0 dBc); check spurs level(-50/-55/-60 dBc)');
ylabel('dBc for continuous-wave signal');
spectrumAnalyzer('signal', s, SAparams, 'fig', 3, 'logscale', false);
ylabel('dBc in 1 Hz');
end
end
function pn_td = PN_src(varargin)
def = {'includeCarrier', true, ...
'spursF_Hz', [], ...
'spursG_dBc', [], ...
'fMax_Hz', []};
p = vararginToStruct(def, varargin);
% length of signal in the time domain (after ifft)
len_s = p.n / p.rate_Hz
% FFT bin frequency spacing
deltaF_Hz = 1 / len_s
% construct AWGN signal in the frequency domain
% a frequency domain bin value of n gives a time domain power of 1
% for example ifft([4 0 0 0]) => 1 1 1 1
% each bin covers a frequency interval of deltaF_Hz
mag = p.n;
% scale "unity power in one bin" => "unity power per Hz":
% multiply with sqrt(deltaF_Hz):
mag = mag * sqrt(deltaF_Hz);
% Create noise according to mag in BOTH real- and imaginary value
mag = mag * sqrt(2);
% both real- and imaginary part contribute unity power => divide by sqrt(2)
pn_fd = mag / sqrt(2) * (randn(1, p.n) + 1i * randn(1, p.n));
% frequency vector corresponding to the FFT bins (0, positive, negative)
fb_Hz = FFT_freqbase(p.n, deltaF_Hz);
% interpolate phase noise spectrum on frequency vector
% note: interpolate dB on logarithmic frequency axis
H_dB = interp1(log(p.f_Hz+eps), p.g_dBc1Hz, log(abs(fb_Hz)+eps), 'linear');
% dB => magnitude
H = 10 .^ (H_dB / 20);
% H = 1e-6; % sanity check: enforce flat -120 dBc in 1 Hz
% apply filter to noise spectrum
pn_fd = pn_fd .* H;
% set spurs
for ix = 1:numel(p.spursF_Hz)
fs = p.spursF_Hz(ix);
u = abs(fb_Hz - fs);
ix2 = find(u == min(u), 1);
% random phase
rot = exp(2i*pi*rand());
% bin value of n: unity (carrier) power (0 dBc)
% scale with sqrt(2) because imaginary part will be discarded
% scale with sqrt(2) because the tone appears at both positive and negative frequencies
smag = 2 * p.n * 10 ^ (p.spursG_dBc(ix) / 20);
pn_fd(ix2) = smag * rot;
end
% limit bandwidth (tool to avoid aliasing in an application
% using the generated phase noise signal)
if ~isempty(p.fMax_Hz)
pn_fd(find(abs(fb_Hz) > p.fMax_Hz)) = 0;
end
% convert to time domain
pn_td = ifft(pn_fd);
% discard imaginary part
pn_td = real(pn_td);
% Now pn_td is a real-valued random signal with a power spectral density
% as specified in f_Hz / g_dBc1Hz.
% phase-modulate to carrier
% note: d/dx exp(x) = 1 near |x| = 1
% in other words, the phase modulation has unity gain for small phase error
pn_td = exp(i*pn_td);
if ~p.includeCarrier
% remove carrier
% returns isolated phase noise component
pn_td = pn_td - 1;
end
end
% returns 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
% *************************************************************
% 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
% discrete-time model for Laplace-domain expression
% Markus Nentwig, 30.12.2011
function sn_model()
close all;
run_demo1(10);
run_demo2(11);
end
% ************************************
% Constructs a FIR model for a relatively
% narrow-band continuous-time IIR filter.
% At the edge of the first Nyquist zone
% (+/- fSample/2), the frequency response
% is down by about 70 dB, which makes
% the modeling unproblematic.
% The impact of different windowing options
% is visible both at high frequencies, but
% also as deviation between original frequency
% response and model at the passband edge.
% ************************************
function run_demo1(fig)
[b, a] = getContTimeExampleFilter();
fc_Hz = 0.5e6; % frequency corresponding to omegaNorm == 1
commonParameters = struct(...
's_a', a, ...
's_b', b, ...
'z_rate_Hz', 3e6, ...
's_fNorm_Hz', fc_Hz, ...
'fig', fig);
% sample impulse response without windowing
ir = sampleLaplaceDomainImpulseResponse(...
commonParameters, ...
'plotstyle_s', 'k-', ...
'plotstyle_z', 'b-');
% use mild windowing
ir = sampleLaplaceDomainImpulseResponse(...
commonParameters, ...
'plotstyle_z', 'r-', ...
'winLen_percent', 4);
% use heavy windowing - error shows at passband edge
ir = sampleLaplaceDomainImpulseResponse(...
commonParameters, ...
'plotstyle_z', 'm-', ...
'winLen_percent', 100);
legend('s-domain', 'sampled, no window', 'sampled, 4% RC window', 'sampled, 100% RC window')
figure(33); clf;
h = stem(ir);
set(h, 'markersize', 2);
set(h, 'lineWidth', 2);
title('sampled impulse response');
end
% ************************************
% model for a continuous-time IIR filter
% that is relatively wide-band.
% The frequency response requires some
% manipulation at the edge of the Nyquist zone.
% Otherwise, there would be an abrupt change
% that would result in an excessively long impulse
% response.
% ************************************
function run_demo2(fig)
[b, a] = getContTimeExampleFilter();
fc_Hz = 1.4e6; % frequency corresponding to omegaNorm == 1
commonParameters = struct(...
's_a', a, ...
's_b', b, ...
'z_rate_Hz', 3e6, ...
's_fNorm_Hz', fc_Hz, ...
'fig', fig);
% sample impulse response without any manipulations
ir1 = sampleLaplaceDomainImpulseResponse(...
commonParameters, ...
'aliasZones', 0, ...
'plotstyle_s', 'k-', ...
'plotstyle_z', 'b-');
% use artificial aliasing (introduces some passband error)
ir2 = sampleLaplaceDomainImpulseResponse(...
commonParameters, ...
'aliasZones', 5, ...
'plotstyle_z', 'r-');
% use artificial low-pass filter (freq. response
% becomes invalid beyond +/- BW_Hz / 2)
ir3 = sampleLaplaceDomainImpulseResponse(...
commonParameters, ...
'aliasZones', 0, ...
'BW_Hz', 2.7e6, ...
'plotstyle_z', 'm-');
line([0, 2.7e6/2, 2.7e6/2], [-10, -10, -50]);
legend('s-domain', ...
sprintf('unmodified (%i taps)', numel(ir1)), ...
sprintf('artificial aliasing (%i taps)', numel(ir2)), ...
sprintf('artificial LP filter (%i taps)', numel(ir3)));
title('2nd example: wide-band filter model frequency response');
figure(350); grid on; hold on;
subplot(3, 1, 1);
stem(ir1, 'b'); xlim([1, numel(ir1)])
title('wide-band filter model: unmodified');
subplot(3, 1, 2);
stem(ir2, 'r');xlim([1, numel(ir1)]);
title('wide-band filter model: art. aliasing');
subplot(3, 1, 3);
stem(ir3, 'm');xlim([1, numel(ir1)]);
title('wide-band filter model: art. LP filter');
end
% Build example Laplace-domain model
function [b, a] = getContTimeExampleFilter()
if true
order = 6;
ripple_dB = 1.2;
omegaNorm = 1;
[b, a] = cheby1(order, ripple_dB, omegaNorm, 's');
else
% same as above, if cheby1 is not available
b = 0.055394;
a = [1.000000 0.868142 1.876836 1.109439 0.889395 0.279242 0.063601];
end
end
% ************************************
% * Samples the impulse response of a Laplace-domain
% component
% * Adjusts group delay and impulse response length so that
% the discarded part of the impulse response is
% below a threshold.
% * Applies windowing
%
% Mandatory arguments:
% s_a, s_b: Laplace-domain coefficients
% s_fNorm_Hz: normalization frequency for a, b coefficients
% z_rate_Hz: Sampling rate of impulse response
%
% optional arguments:
% truncErr_dB: Maximum truncation error of impulse response
% nPts: Computed impulse response length before truncation
% winLen_percent: Part of the IR where windowing is applied
% BW_Hz: Apply artificical lowpass filter for |f| > +/- BW_Hz / 2
%
% plotting:
% fig: Figure number
% plotstyle_s: set to plot Laplace-domain frequency response
% plotstyle_z: set to plot z-domain model freq. response
% ************************************
function ir = sampleLaplaceDomainImpulseResponse(varargin)
def = {'nPts', 2^18, ...
'truncErr_dB', -60, ...
'winLen_percent', -1, ...
'fig', 99, ...
'plotstyle_s', [], ...
'plotstyle_z', [], ...
'aliasZones', 1, ...
'BW_Hz', []};
p = vararginToStruct(def, varargin);
% FFT bin frequencies
fbase_Hz = FFT_frequencyBasis(p.nPts, p.z_rate_Hz);
% instead of truncating the frequency response at +/- z_rate_Hz,
% fold the aliases back into the fundamental Nyquist zone.
% Otherwise, we'd try to build a near-ideal lowpass filter,
% which is essentially non-causal and requires a long impulse
% response with artificially introduced group delay.
H = 0;
for alias = -p.aliasZones:p.aliasZones
% Laplace-domain "s" in normalized frequency
s = 1i * (fbase_Hz + alias * p.z_rate_Hz) / p.s_fNorm_Hz;
% evaluate polynomial in s
H_num = polyval(p.s_b, s);
H_denom = polyval(p.s_a, s);
Ha = H_num ./ H_denom;
H = H + Ha;
end
% plot |H(f)| if requested
if ~isempty(p.plotstyle_s)
figure(p.fig); hold on; grid on;
fr = fftshift(20*log10(abs(H) + eps));
h = plot(fftshift(fbase_Hz), fr, p.plotstyle_s);
set(h, 'lineWidth', 3);
xlabel('f/Hz'); ylabel('|H(f)| / dB');
xlim([0, p.z_rate_Hz/2]);
end
% apply artificial lowpass filter
if ~isempty(p.BW_Hz)
% calculate RC transition bandwidth
BW_RC_trans_Hz = p.z_rate_Hz - p.BW_Hz;
% alpha (RC filter parameter) implements the
% RC transition bandwidth:
alpha_RC = BW_RC_trans_Hz / p.z_rate_Hz / 2;
% With a cutoff frequency of BW_RC, the RC filter
% reaches |H(f)| = 0 at f = z_rate_Hz / 2
% BW * (1+alpha) = z_rate_Hz / 2
BW_RC_Hz = p.z_rate_Hz / ((1+alpha_RC));
HRC = raisedCosine(fbase_Hz, BW_RC_Hz, alpha_RC);
H = H .* HRC;
end
% frequency response => impulse response
ir = ifft(H);
% assume s_a, s_b describe a real-valued impulse response
ir = real(ir);
% the impulse response peak is near the first bin.
% there is some earlier ringing, because evaluating the s-domain
% expression only for s < +/- z_rate_Hz / 2 implies an ideal,
% non-causal low-pass filter.
ir = fftshift(ir);
% now the peak is near the center
% threshold for discarding samples
peak = max(abs(ir));
thr = peak * 10 ^ (p.truncErr_dB / 20);
% first sample that is above threshold
% determines the group delay of the model
ixFirst = find(abs(ir) > thr, 1, 'first');
% last sample above threshold
% determines the length of the impulse response
ixLast = find(abs(ir) > thr, 1, 'last');
% truncate
ir = ir(ixFirst:ixLast);
% apply windowing
if p.winLen_percent > 0
% note: The window drops to zero for the first sample that
% is NOT included in the impulse response.
v = linspace(-1, 0, numel(ir)+1);
v = v(1:end-1);
v = v * 100 / p.winLen_percent;
v = v + 1;
v = max(v, 0);
win = (cos(v*pi) + 1) / 2;
ir = ir .* win;
end
% plot frequency response, if requested
if ~isempty(p.plotstyle_z)
irPlot = zeros(1, p.nPts);
irPlot(1:numel(ir)) = ir;
figure(p.fig); hold on;
fr = fftshift(20*log10(abs(fft(irPlot)) + eps));
h = plot(fftshift(fbase_Hz), fr, p.plotstyle_z);
set(h, 'lineWidth', 3);
xlabel('f/Hz'); ylabel('|H(f)| / dB');
xlim([0, p.z_rate_Hz/2]);
end
end
% ************************************
% raised cosine frequency response
% ************************************
function H = raisedCosine(f, BW, alpha)
c=abs(f * 2 / BW);
% c=-1 for lower end of transition region
% c=1 for upper end of transition region
c=(c-1) / alpha;
% clip to -1..1 range
c=min(c, 1);
c=max(c, -1);
H = 1/2+cos(pi/2*(c+1))/2;
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
% *************************************************************
% 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
% Code for testing the 'Wind_Flattop(Spec)' function in
% reducing 'scalloping loss' errors in time signal amplitude
% estimation.
%
% Generates a time-domain sinusoid, computes its FFT, and passes
% that FFT sequence to the 'Wind_Flattop(Spec)' function.
%
% The maximum output sample of the 'Wind_Flattop(Spec)' function
% is used to estimate the peak amplitude of the original
% sinusoidal time-domain test signal.
%
% The user controls the values for the test sinusoid's
% 'Test_Freq' and peak amplitude 'Peak_Amp' near lines 23 & 24.
%
% Richard Lyons [December, 2011]
clear all, clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define test parameters
Test_Freq = 7.22; % Test tone's freq. Must be less that N/2
Peak_Amp = 5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 64; % Number of time samples
Index = (0:N-1);
X = Peak_Amp*cos(2*pi*(Test_Freq)*Index/N + pi/3);
figure(1), clf
subplot(2,1,1)
plot(X,':ko', 'markersize', 4)
title('Original time signal'), grid on, zoom on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT the input 'X' sequence and call 'Wind_Flattop()' function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Spec = fft(X);
[Windowed_Spec] = Wind_Flattop(Spec);
subplot(2,1,2), plot(abs(Spec),'ko', 'markersize', 3)
title('SpecMag of unwindowed time signal'), grid on, zoom on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display results accuracy (error)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
disp(['Test Freq = ',num2str(Test_Freq),...
', True Peak Amplitude = ',num2str(Peak_Amp)])
Mag_peak_unwindowed = max(abs(Spec));
Unwindowed_Amp_Estimate = 2*Mag_peak_unwindowed/N;
Unwindowed_Amp_Estimate_Error_in_dB = ...
20*log10(Unwindowed_Amp_Estimate/Peak_Amp);
disp(' ')
disp(['Unwindowed Peak Amplitude Estimate = ',...
num2str(Unwindowed_Amp_Estimate)])
disp(['Unwindowed Estimate Error in dB = ',...
num2str(Unwindowed_Amp_Estimate_Error_in_dB),' dB'])
M_peak_windowed = max(abs(Windowed_Spec));
Windowed_Amp_Estimated = 2*M_peak_windowed/N;
Windowed_Amp_Estimation_Error_in_dB = ...
20*log10(Windowed_Amp_Estimated/Peak_Amp);
disp(' ')
disp(['Windowed Peak Amplitude Estimate = ',...
num2str(Windowed_Amp_Estimated)])
disp(['Windowed Estimate Error in dB = ',...
num2str(Windowed_Amp_Estimation_Error_in_dB),' dB'])
function [Windowed_Spec] = Wind_Flattop(Spec)
% Given an input spectral sequence 'Spec', that is the
% FFT of some time sequence 'x', Wind_Flattop(Spec)
% returns a spectral sequence that is equivalent
% to the FFT of a flat-top windowed version of time
% sequence 'x'. The peak magnitude values of output
% sequence 'Windowed_Spec' can be used to accurately
% estimate the peak amplitudes of sinusoidal components
% in time sequence 'x'.
% Input: 'Spec' (a sequence of complex FFT sample values)
%
% Output: 'Windowed_Spec' (a sequence of complex flat-top
% windowed FFT sample values)
%
% Based on Lyons': "Reducing FFT Scalloping Loss Errors
% Without Multiplication", IEEE Signal Processing Magazine,
% DSP Tips & Tricks column, March, 2011, pp. 112-116.
%
% Richard Lyons [December, 2011]
N = length(Spec);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform freq-domain convolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g_Coeffs = [1, -0.94247, 0.44247];
% Compute first two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(1) = g_Coeffs(3)*Spec(N-1) ...
+g_Coeffs(2)*Spec(N) + Spec(1) ...
+g_Coeffs(2)*Spec(2) + g_Coeffs(3)*Spec(3);
Windowed_Spec(2) = g_Coeffs(3)*Spec(N) ...
+g_Coeffs(2)*Spec(1) + Spec(2) ...
+g_Coeffs(2)*Spec(3) + g_Coeffs(3)*Spec(4);
% Compute last two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(N-1) = g_Coeffs(3)*Spec(N-3) ...
+g_Coeffs(2)*Spec(N-2) + Spec(N-1) ...
+g_Coeffs(2)*Spec(N) + g_Coeffs(3)*Spec(1);
Windowed_Spec(N) = g_Coeffs(3)*Spec(N-2) ...
+g_Coeffs(2)*Spec(N-1) + Spec(N) ...
+g_Coeffs(2)*Spec(1) + g_Coeffs(3)*Spec(2);
% Compute convolved spec samples for the middle of the spectrum
for K = 3:N-2
Windowed_Spec(K) = g_Coeffs(3)*Spec(K-2) ...
+g_Coeffs(2)*Spec(K-1) + Spec(K) ...
+g_Coeffs(2)*Spec(K+1) + g_Coeffs(3)*Spec(K+2);
end % % End of 'Wind_Flattop(Spec)' function
% ************************************
% 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
