Function to perform Adaptive filtering

Senthilkumar December 28, 2011 Coded in Scilab
function[y,M] = adapt_filt(xlms,B,h,delta,l,l1)
    //x = the signal from the speaker directly
    //B = the signal thorugh the echo path
    //h = impulse response of adaptive filter
    //l = length of signal 'x'
    //l1 = length of adaptive filter order
    for k = 1:150
        for n = 1:l
            xcap = xlms((n+l1-1):-1:(n+l1-1)-l1+1);
            yout(n) = h*xcap;
            e(n) = B(n)- yout(n);
            xnorm = 0.001+(xcap*xcap');
            h = h+((delta*e(n))*(xcap'));
        end
        eold = 0.0;
        for i = 1:l
            MSE = eold+(e(i)^2);
            eold = MSE;
        end
        if MSE <= 0.0001 then
            break;
        end
    end
    y =  zeros(1,length(e));
    M = zeros(1,length(h));
    y  = e;
    M = h;
endfunction

Full band Echo canceller

Senthilkumar December 28, 2011 Coded in Scilab
//Caption : FULL Band Echo Cnacellation using NLMS Adaptive Filter

clc;
//Reading a speech signal
[x,Fs,bits]=wavread("E:\4.wav");
order = 40;  // Adaptive filter order
x = x';
N = length(x);  //length of speech signal
//Delay introduced in echo path
delay = 100;   
//Echo at same speaker
xdelayed = zeros(1,N+delay);
for i = delay+1:N+delay
    xdelayed(i) = x(i-delay);
end
//Initialize the adaptive filter coefficients to zero
hcap = zeros(1,order);
//To avoid negative values generated during convolution
for i = 1:order-1
    xlms(i) = 0.0;
end
for i = order:N+order-1
    xlms(i) = x(i-order+1);
end
Power_X = pow_1(x,N);  //Average power of speech signal
//Calculation of step size of adaptive filter
delta = 1/(10*order*Power_X);
[x_out,Adapt_Filter_IR] = adapt_filt(xlms,xdelayed,hcap,delta,N,order)
figure(1)
subplot(3,1,1)
plot([1:N],x)
title('Speech Signal Generated by some Speaker A')
sound(x,Fs,16)
subplot(3,1,2)
plot([1:length(xdelayed)],xdelayed)
title('Echo signal of speaker A received by speaker A')
sound(xdelayed,Fs,16)
subplot(3,1,3)
plot([1:length(x_out)],x_out)
title('Echo signal at speaker A after using NLMS Adaptive filter Echo Canceller')
figure(2)
plot2d3('gnn',[1:length(Adapt_Filter_IR)],Adapt_Filter_IR)
title('Adaptive Filter (Echo Canceller) Impulse response')

Chaning the bit depth in speech samples

Senthilkumar December 28, 2011 Coded in Scilab
//Caption: Reading a Speech Signal & 
//[1]. Write it in another file
//[2]. Changing the bit depth from 16 bits to 8 bits
clear;
clc;
[y,Fs,bits_y]=wavread("E:\4.wav");
wavwrite(y,Fs,8,"E:\4_8bit.wav");
[x,Fs,bits_x]=wavread("E:\4_8bit.wav");
Ny = length(y); //Number of samples in y (4.wav)
Nx = length(x); //Number of samples in x (4_8bit.wav)
Memory_y = Ny*bits_y;  //memory requirement for 4.wav in bits
Memory_x = Nx*bits_x;  //memory requirement for 4_8bit.wav in bits
disp(Memory_y,'memory requirement for 4.wav in bits =')
disp(Memory_x,'memory requirement for 4_8bit.wav in bits =')
//Result
//memory requirement for 4.wav in bits =   
//     133760.  
// 
//memory requirement for 4_8bit.wav in bits =   
//     66880.

Alias/error simulation of interpolating RRC filter

Markus Nentwig December 25, 2011 Coded in Matlab
% *************************************************
% 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 phase noise model

Markus Nentwig December 18, 20114 comments Coded in Matlab
% ****************************************************************
% 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

Digital model for analog filters

Markus Nentwig December 17, 2011 Coded in Matlab
% 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

Python zplane function

Christopher Felton December 17, 20113 comments Coded in Python
#
# Copyright (c) 2011 Christopher Felton
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

# The following is derived from the slides presented by
# Alexander Kain for CS506/606 "Special Topics: Speech Signal Processing"
# CSLU / OHSU, Spring Term 2011.

import numpy as np
import matplotlib.pyplot as plt
from  matplotlib import patches
from matplotlib.figure import Figure
from matplotlib import rcParams
    
def zplane(b,a,filename=None):
    """Plot the complex z-plane given a transfer function.
    """

    # get a figure/plot
    ax = plt.subplot(111)

    # create the unit circle
    uc = patches.Circle((0,0), radius=1, fill=False,
                        color='black', ls='dashed')
    ax.add_patch(uc)

    # The coefficients are less than 1, normalize the coeficients
    if np.max(b) > 1:
        kn = np.max(b)
        b = b/float(kn)
    else:
        kn = 1

    if np.max(a) > 1:
        kd = np.max(a)
        a = a/float(kd)
    else:
        kd = 1
        
    # Get the poles and zeros
    p = np.roots(a)
    z = np.roots(b)
    k = kn/float(kd)
    
    # Plot the zeros and set marker properties    
    t1 = plt.plot(z.real, z.imag, 'go', ms=10)
    plt.setp( t1, markersize=10.0, markeredgewidth=1.0,
              markeredgecolor='k', markerfacecolor='g')

    # Plot the poles and set marker properties
    t2 = plt.plot(p.real, p.imag, 'rx', ms=10)
    plt.setp( t2, markersize=12.0, markeredgewidth=3.0,
              markeredgecolor='r', markerfacecolor='r')

    ax.spines['left'].set_position('center')
    ax.spines['bottom'].set_position('center')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)

    # set the ticks
    r = 1.5; plt.axis('scaled'); plt.axis([-r, r, -r, r])
    ticks = [-1, -.5, .5, 1]; plt.xticks(ticks); plt.yticks(ticks)

    if filename is None:
        plt.show()
    else:
        plt.savefig(filename)
    

    return z, p, k

Testing the Flat-Top Windowing Function

Rick Lyons December 14, 2011 Coded in Matlab
%    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'])

Flat-Top Windowing Function for the Accurate Measurement of a Sinusoid's Peak Amplitude Based on FFT Data

Rick Lyons December 14, 2011 Coded in Matlab
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

16QAM Modem model (Basic)

Kadhiem Ayob December 12, 2011 Coded in Matlab
clear all; close all;

n = 2^16;   %number of symbols
Fsym = 12.5; %Msps
Fs = 100;   %MHz, IF sampling frequency
Fc = 20;    %MHz, upconverter frequency

%generate 16QAM raw symbols
alphabet = [-1 -1/3 +1/3 +1];
tx = complex(randsrc(1,n,alphabet),randsrc(1,n,alphabet));

%pulse shaping, root raised cosine
h = firrcos(50,Fsym/4,.15,Fsym,'rolloff','sqrt');
tx_up1 = zeros(1,2*n);
tx_up1(1:2:end-1) = tx;
tx_shaped = filter(2*h,1,tx_up1);

%further upsampling by 4
tx_up2 = resample(tx_shaped,4,1);

%upconvert to 10MHz centre
f1 = exp(1i*2*pi*(0:8*n-1)*Fc/Fs);
tx_upconverted = tx_up2 .* f1;

%channel signal
rx_real = real(tx_upconverted);

%Rx shifts signal back to zero 
f2 = exp(1i*2*pi*(0:8*n-1)*-Fc/Fs);
rx_downconverted = rx_real.*f2;

%Rx downsamples back by 4
rx_dn1 = resample(rx_downconverted,1,4);

%rx matched filter
rx_dn2 = filter(2*h,1,rx_dn1);
rx_dn2 = rx_dn2(50:end); %remove initial zeros

%one phase should be correct(odd/even)
rx_e = rx_dn2(1:2:end); %odd
rx_o = rx_dn2(2:2:end); %even

subplot(3,1,1);plot(real(tx),imag(tx),'o');
axis([-2 +2 -2 +2]);grid;
title('Tx constellations')
subplot(3,1,2);plot(real(rx_e),imag(rx_e),'o');
axis([-2 +2 -2 +2]);grid;
title('Rx constellations, first phase')
subplot(3,1,3);plot(real(rx_o),imag(rx_o),'o');
axis([-2 +2 -2 +2]);grid;
title('Rx constellations, second phase')