% 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
