Code Snippets

Power Computation of a Digital Stream

kaz - ● February 26, 2012 ● Coded in Matlab 
clear all;

%example vector having ac & dc power
x = complex(randn(1,2048)+.2,randn(1,2048)-.1);

%total power, three equivalent methods
pwr1_total = mean(abs(x).^2);                %mean of squared values
pwr2_total = mean(real(x).^2 + imag(x).^2);
pwr3_total = mean(x .* conj(x));

%total expressed as rms
rms_total = sqrt(pwr1_total);

%dc power
pwr1_dc = mean(real(x))^2 + mean(imag(x))^2; %square of mean of values

%ac power
pwr1_ac = pwr1_total - pwr1_dc;              %mean of squares - square of mean

%relation of ac power to statistical variance 
pwr2_ac = var(x);              %approximately

%ac expressed as rms
rms1_ac = sqrt(pwr1_ac);

%ac relation to standard of deviation, std = sqrt(var)
rms2_ac = std(x);              %approximately

%dc relation to variance
pwr2_dc = pwr1_total - var(x); %approximately

fprintf('----------------------------------------------------\r');
fprintf('power(total),          : %1.5f, %1.5f, %1.5f\r',pwr1_total,pwr2_total,pwr3_total);
fprintf('rms(total)             : %1.5f\r',rms_total);
fprintf('power(ac),  variance   : %1.5f, %1.5f\r',pwr1_ac,pwr2_ac);
fprintf('rms(ac),    std        : %1.5f, %1.5f\r',rms1_ac,rms2_ac);
fprintf('power(dc), (total-var) : %1.5f, %1.5f\r',pwr1_dc,pwr2_dc);

Halfband Filter Design with Python/Scipy

Christopher Felton ● February 24, 2012●5 comments ● Coded in Python 
# The following is a Python/scipy snippet to generate the 
# coefficients for a halfband filter.  A halfband filter
# is a filter where the cutoff frequency is Fs/4 and every
# other coeffecient is zero except the cetner tap.
# Note: every other (even except 0) is 0, most of the coefficients
#       will be close to zero, force to zero actual

import numpy
from numpy import log10, abs, pi
import scipy
from scipy import signal
import matplotlib
import matplotlib.pyplot
import matplotlib as mpl

# ~~[Filter Design with Parks-McClellan Remez]~~
N = 32  # Filter order
# Filter symetric around 0.25 (where .5 is pi or Fs/2)
bands = numpy.array([0., .22, .28, .5])
h = signal.remez(N+1, bands, [1,0], [1,1])
h[abs(h) <= 1e-4] = 0.
(w,H) = signal.freqz(h)

# ~~[Filter Design with Windowed freq]~~
b = signal.firwin(N+1, 0.5)
b[abs(h) <= 1e-4] = 0.
(wb, Hb) = signal.freqz(b)

# Dump the coefficients for comparison and verification
print('          remez       firwin')
print('------------------------------------')
for ii in range(N+1):
    print(' tap %2d   %-3.6f    %-3.6f' % (ii, h[ii], b[ii]))

## ~~[Plotting]~~
# Note: the pylab functions can be used to create plots,
#       and these might be easier for beginners or more familiar
#       for Matlab users.  pylab is a wrapper around lower-level
#       MPL artist (pyplot) functions.
fig = mpl.pyplot.figure()
ax0 = fig.add_subplot(211)
ax0.stem(numpy.arange(len(h)), h)
ax0.grid(True)
ax0.set_title('Parks-McClellan (remez) Impulse Response')
ax1 = fig.add_subplot(212)
ax1.stem(numpy.arange(len(b)), b)
ax1.set_title('Windowed Frequency Sampling (firwin) Impulse Response')
ax1.grid(True)
fig.savefig('hb_imp.png')

fig = mpl.pyplot.figure()
ax1 = fig.add_subplot(111)
ax1.plot(w, 20*log10(abs(H)))
ax1.plot(w, 20*log10(abs(Hb)))
ax1.legend(['remez', 'firwin'])
bx = bands*2*pi
ax1.axvspan(bx[1], bx[2], facecolor='0.5', alpha='0.33')
ax1.plot(pi/2, -6, 'go')
ax1.axvline(pi/2, color='g', linestyle='--')
ax1.axis([0,pi,-64,3])
ax1.grid('on')
ax1.set_ylabel('Magnitude (dB)')
ax1.set_xlabel('Normalized Frequency (radians)')
ax1.set_title('Half Band Filter Frequency Response')
fig.savefig('hb_rsp.png')

Computing CIC Filter Register Pruning Using Matlab [Updated]

Rick Lyons ● February 15, 2012●8 comments ● Coded in Matlab 
%  Filename: CIC_Word_Truncation.m
%
%  Computes CIC decimation filter accumulator register  
%  truncation in each filter stage based on Hogenauer's  
%  'accumulator register pruning' technique.
%
%   Inputs:
%     N = number of decimation CIC filter stages (filter order).
%     R = CIC filter rate change factor (decimation factor).
%     M = differential delay.
%     Bin = number of bits in an input data word.
%     Bout = number of bits in the filter's final output data word.
%   Outputs:
%     Stage number (ranges from 1 -to- 2*N+1).
%     Bj = number of least significant bits that can be truncated
%       at the input of a filter stage.
%     Accumulator widths = number of a stage's necessary accumulator
%       bits accounting for truncation.
 
%  Richard Lyons Feb., 2012
 
clear, clc
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Define CIC filter parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%N = 4;  R = 25;  M = 1;  Bin = 16;  Bout = 16; % Hogenauer paper, pp. 159
%N = 3;  R = 32;  M = 2;  Bin = 8;  Bout = 10; % Meyer Baese book, pp. 268
%N = 3;  R = 16;  M = 1;  Bin = 16;  Bout = 16; % Thorwartl's PDF file

N = 3; R = 8; M = 1; Bin = 12; Bout = 12; % Lyons' blog Figure 2 example
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find h_sub_j and "F_sub_j" values for (N-1) cascaded integrators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    disp(' ')
    for j = N-1:-1:1
        h_sub_j = [];
        h_sub_j((R*M-1)*N + j -1 + 1) = 0;
        for k = 0:(R*M-1)*N + j -1
            for L = 0:floor(k/(R*M)) % Use uppercase "L" for loop variable
                Change_to_Result = (-1)^L*nchoosek(N, L)...
                                  *nchoosek(N-j+k-R*M*L,k-R*M*L);
                h_sub_j(k+1) =  h_sub_j(k+1) + Change_to_Result;
            end % End "L" loop
        end % End "k" loop
        F_sub_j(j) = sqrt(sum(h_sub_j.^2));
    end % End "j" loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define "F_sub_j" values for up to seven cascaded combs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j_for_many_combs = sqrt([2, 6, 20, 70, 252, 924, 3432]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Compute F_sub_j for last integrator stage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j(N) = F_sub_j_for_many_combs(N-1)*sqrt(R*M);  % Last integrator   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Compute F_sub_j for N cascaded filter's comb stages
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2*N:-1:N+1
    F_sub_j(j) = F_sub_j_for_many_combs(2*N -j + 1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define "F_sub_j" values for the final output register truncation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j(2*N+1) = 1; % Final output register truncation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute column vector of minus log base 2 of "F_sub_j" values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Minus_log2_of_F_sub_j = -log2(F_sub_j)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute total "Output_Truncation_Noise_Variance" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CIC_Filter_Gain = (R*M)^N;
Num_of_Bits_Growth = ceil(log2(CIC_Filter_Gain));
% The following is from Hogenauer's Eq. (11)
    %Num_Output_Bits_With_No_Truncation = Num_of_Bits_Growth +Bin -1;
    Num_Output_Bits_With_No_Truncation = Num_of_Bits_Growth +Bin;
Num_of_Output_Bits_Truncated = Num_Output_Bits_With_No_Truncation -Bout;
Output_Truncation_Noise_Variance = (2^Num_of_Output_Bits_Truncated)^2/12;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute log base 2 of "Output_Truncation_Noise_Standard_Deviation" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Output_Truncation_Noise_Standard_Deviation = ...
    sqrt(Output_Truncation_Noise_Variance);
Log_base2_of_Output_Truncation_Noise_Standard_Deviation = ...
    log2(Output_Truncation_Noise_Standard_Deviation);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute column vector of "half log base 2 of 6/N" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Half_Log_Base2_of_6_over_N = 0.5*log2(6/N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute desired "B_sub_j" vector
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B_sub_j = floor(Minus_log2_of_F_sub_j ...
          + Log_base2_of_Output_Truncation_Noise_Standard_Deviation ...
          + Half_Log_Base2_of_6_over_N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' '), disp(' ')
disp(['N = ',num2str(N),',   R = ',num2str(R),',   M = ',num2str(M),...
        ',   Bin = ', num2str(Bin),',   Bout = ',num2str(Bout)])
disp(['Num of Bits Growth Due To CIC Filter Gain = ', ...
    num2str(Num_of_Bits_Growth)])
disp(['Num of Accumulator Bits With No Truncation = ', ...
    num2str(Num_Output_Bits_With_No_Truncation)])
% disp(['Output Truncation Noise Variance = ', ...
%     num2str(Output_Truncation_Noise_Variance)])
% disp(['Log Base2 of Output Truncation Noise Standard Deviation = ',...
%         num2str(Log_base2_of_Output_Truncation_Noise_Standard_Deviation)])
% disp(['Half Log Base2 of 6/N = ', num2str(Half_Log_Base2_of_6_over_N)])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create and display "Results" matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for Stage = 1:2*N
    Results(Stage,1) = Stage;
    Results(Stage,2) = F_sub_j(Stage);
    Results(Stage,3) = Minus_log2_of_F_sub_j(Stage);
    Results(Stage,4) = B_sub_j(Stage);
    Results(Stage,5) = Num_Output_Bits_With_No_Truncation -B_sub_j(Stage);
end
% Include final output stage truncation in "Results" matrix
    Results(2*N+1,1) = 2*N+1;  % Output stage number
    Results(2*N+1,2) = 1;
    Results(2*N+1,4) = Num_of_Output_Bits_Truncated;
    Results(2*N+1,5) = Bout;
    %Results % Display "Results" matrix in raw float-pt.form
 
% % Display "F_sub_j" values if you wish
% disp(' ')
% disp(' Stage        Fj        -log2(Fj)    Bj   Accum width')
% for Stage = 1:2*N+1
%   disp(['  ',sprintf('%2.2g',Results(Stage,1)),sprintf('\t'),sprintf('%12.3g',Results(Stage,2)),...
%         sprintf('\t'),sprintf('%7.5g',Results(Stage,3)),sprintf('\t'),...
%         sprintf('%7.5g',Results(Stage,4)),sprintf('\t'),sprintf('%7.5g',Results(Stage,5))])
% end
 
% Display Stage number, # of truncated input bits, & Accumulator word widths
disp(' ')
disp(' Stage(j)   Bj   Accum (adder) width')
for Stage = 1:2*N
    disp(['  ',sprintf('%2.0f',Results(Stage,1)),...
        sprintf('\t'),...
        sprintf('%5.5g',Results(Stage,4)),sprintf('\t'),...
        sprintf('%7.5g',Results(Stage,5))])
end
    disp(['  ',sprintf('%2.0f',Results(2*N+1,1)),...
        sprintf('\t'),...
        sprintf('%5.5g',Results(2*N+1,4)),sprintf('\t'),...
        sprintf('%7.5g',Results(2*N+1,5)),' (final truncation)'])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Reading text files in Matlab

kaz - ● February 13, 2012 ● Coded in Matlab 
clear all;

% create test file
data = [1 5 2 4 3 3 4 2 5 1];
filename = 'test_file.txt';
fid = fopen(filename, 'w');
fprintf(fid, '%d %d\n', data);
fclose(fid);

% (1) read file using fscanf 
fid = fopen(filename, 'r'); 
y1 = fscanf(fid, '%d\n'); %interleaves columns
fclose(fid);

% (2) read file using textread (or textscan)
[ya,yb] = textread(filename,'%d%d');
y2 = [ya yb];

% (3) read file using importdata 
y3 = importdata(filename);

% (4) read file using load 
y4 = load(filename);

disp('-------------------------------------------------------------')
disp('    original vector data')
disp(data)
disp('    file content using fprintf')
disp(y2)
disp('    vector created by fscanf')
disp(y1)
disp('     matrix created by:')
disp('     textread   importdata    load')
disp([y2 y3 y4])

Ideal interpolation filter

kaz - ● January 21, 2012 ● Coded in Matlab 
clear all; close all;

%example bandlimited random input & parameters
x = filter(fir1(70,.1),1,randn(1,1024));

up = 3;             %Interpolation factor
cutoff = .3;

intorder = 6;
h1 = intfilt(up, intorder, cutoff); %ideal filter
h1 = up*h1/sum(h1);

h2 = fir1(2*up*intorder-2,cutoff);  %ordinary LPF 
h2 = up*h2/sum(h2);

%upsample
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;

x_f1 = filter(h1,1,x_up);
x_f2 = filter(h2,1,x_up);

figure;
subplot(3,1,1);hold
plot(x_f1,'o--');
plot(x_f2,'r.-');
legend('ideal output','fir1 output');

subplot(3,1,2);
plot(x_f1-x_f2);
legend('error');

subplot(3,1,3);hold
plot(h1,'.-');
plot(h2,'r.-');
legend('ideal filter','fir1 filter');

checking resampling in time domain

kaz - ● January 21, 2012 ● Coded in Matlab 
clear all; close all;

%example bandlimited random input & parameters
x = filter(fir1(70,.1),1,randn(1,512));
h = fir1(30,.3);    %filter used for upsampling
up = 3;             %Interpolation factor
dn = 2;             %Decimation factor

%%%%%%%%%%%%%%%%%%%%% up/dn model %%%%%%%%%%%%%%%%%%%%%%
%upsample input
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;
x_f = filter(up*h,1,x_up);

%downsample signal by decimation
x_dn = x_f(1:dn:end);

delay = 30/2+1;  
figure;
subplot(2,1,1);hold;
plot(x_up,'o--');
plot(x_f(delay:end),'r.-');
legend('input with zeros','upsampled stage');

subplot(2,1,2);hold;
plot(x_up(1:dn:end),'o--');
plot(x_dn(ceil(delay/dn):end),'r.-');
legend('input with zeros','final signal');

Resampling filter performance

kaz - ● January 14, 2012●3 comments ● Coded in Matlab 
%%%%%%%%%%%%%%%%%% inputs for model %%%%%%%%%%%%%%%%
clear all; close all;

%example bandlimited impulse input & parameters
x = filter(fir1(70,.1),1,[1 zeros(1,2^15-1)]);
Fs = 120;           %MHz original sample rate             
h = fir1(30,.3);    %filter used for upsampling
up = 3;             %Interpolation factor
dn = 2;             %Decimation factor
Fc = 12;             %MHz band centre (-Fs/2 ~ +Fs/2)
Fch = 0;             %MHz filter centre (-Fs*up/2 ~ +Fs*up/2)

%move signal to its centre
x = x.*exp(j*2*pi*(0:length(x)-1)*Fc/Fs);

%shift filter 
h = h.*exp(j*2*pi*(0:length(h)-1)*Fch/(Fs*up));

%%%%%%%%%%%%%%%%%%%%% model %%%%%%%%%%%%%%%%%%%%%%
%check signal in upsampled domain
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;
[P, F] = pwelch(x_up, [], 0, 2^16, Fs*up,'twosided'); 
F = F - max(F)/2; 
P = fftshift(P);
y(find(P == 0)) = -100;  %avoid log of zero 
y(find(P ~= 0)) = 10*log10(P(find(P ~= 0)));
P_dB = y - 10*log10(max(P));  %normalise

%check filter response in upsampled domain
H = fftshift(20*log10(abs(fft(h,2^16))));

subplot(2,1,1);
hold;grid;
plot(F, P_dB,'.-'); 
plot(F,H,'m--');
axis([min(F)-1 max(F)+1 -80 1]);
legend('upsampled signal','upsampling filter');

%check signal in downsampled domain
x_f = filter(h,1,x_up);
x_dn = x_f(1:dn:end);
[P, F] = pwelch(x_dn, [], 0, 2^16, Fs*up/dn,'twosided'); 
F = F - max(F)/2;
P = fftshift(P);
y(find(P == 0)) = -100;   %avoid log of zero 
y(find(P ~= 0)) = 10*log10(P(find(P ~= 0)));
P_dB = y - 10*log10(max(P)); %normalise

subplot(2,1,2)
plot(F,P_dB,'r.-');
grid;
axis([min(F)-1 max(F)+1 -80 1]);
legend('downsampled signal');

Adding a Controlled Amount of Noise to a Noise-Free Signal

Rick Lyons ● January 3, 2012●21 comments ● Coded in Matlab 
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)

Narrow-band moving-average decimator, one addition/sample

Markus Nentwig ● December 31, 2011 ● Coded in Matlab 
% *************************************************
% Moving average decimator 
%
% A decimator for narrow-band signals (~5 % or less bandwidth occupation) 
% can be implemented as follows:
%
% #define DECIM (100)
% double acc = 0.0;
% while (1){
%    int ix;
%    for(ix = DECIM; ix > 0; --ix){
%       acc += getInputSample();
%    } /* for */
%    writeOutputSample(acc / (double)DECIM);
%    acc = 0.0;
% } /* while */
%
% It is conceptually identical to a moving average filter
% http://www.dspguide.com/ch15/4.htm combined with a decimator
%
% Note that the "moving" average jumps ahead in steps of the decimation
% factor. Intermediate output is decimated away, allowing for a very efficient 
% implementation.
% This program calculates the frequency response and alias response,
% based on the decimation factor and bandwidth of the processed signal.
% *************************************************    
function eval_design()
    decimationFactor = 100;
    rateIn_Hz = 48000;
    noDecim = false;
    
    % create illustration with sinc response
    %decimationFactor = 4; noDecim = true; 
    
    % *************************************************
    % signal source: Bandlimited test pulse
    % Does not contain energy in frequency ranges that 
    % cause aliasing
    % *************************************************
    s = zeros(1, 10000 * decimationFactor);    
    fb = FFT_frequencyBasis(numel(s), rateIn_Hz);
    
    % assign energy to frequency bins
    if noDecim
        sPass = ones(size(s));
    else
        sPass = s; sPass(find(abs(fb) < rateIn_Hz / decimationFactor / 2)) = 1;
    end
    sAlias = s; sAlias(find(abs(fb) >= rateIn_Hz / decimationFactor / 2)) = 1;
    
    % convert to time domain
    sPass = fftshift(real(ifft(sPass)));
    sAlias = fftshift(real(ifft(sAlias)));

    % *************************************************    
    % plot spectrum at input
    % *************************************************    
    pPass = {};
    pPass = addPlot(pPass, sPass, rateIn_Hz, 'k', 5, ...
                    'input (passband response)');
    pAlias = {};
    pAlias = addPlot(pAlias, sAlias, rateIn_Hz, 'k', 5, ...
                     'input (alias response)');    
        
    % *************************************************    
    % impulse response
    % *************************************************    
    h = zeros(size(s));
    h (1:decimationFactor) = 1; 
    if noDecim
        h = h / decimationFactor;
        decimationFactor = 1;
    end
    
    % cyclic convolution between signal and impulse response
    sPass = real(ifft(fft(sPass) .* fft(h)));
    sAlias = real(ifft(fft(sAlias) .* fft(h)));
    
    % decimation
    sPass = sPass(decimationFactor:decimationFactor:end);
    sAlias = sAlias(decimationFactor:decimationFactor:end);
    rateOut_Hz = rateIn_Hz / decimationFactor;
    
    % *************************************************    
    % plot spectrum
    % *************************************************    
    pPass = addPlot(pPass, sPass, rateOut_Hz, 'b', 3, ...
                    'decimated (passband response)');
    figure(1); clf; grid on; hold on;
    doplot(pPass, sprintf('passband frequency response over input rate, decim=%i', decimationFactor));
    
    pAlias = addPlot(pAlias, sAlias, rateOut_Hz, 'b', 3, ...
                     'decimated (alias response)');
    figure(2); clf; grid on; hold on;
    doplot(pAlias, sprintf('alias frequency response over input rate, decim=%i', decimationFactor));

    % *************************************************    
    % plot passband ripple
    % *************************************************    
    fb = FFT_frequencyBasis(numel(sPass), 1);
    fr = 20*log10(abs(fft(sPass) + eps));
    ix = find(fb > 0);

    figure(3); clf;
    h = semilogx(fb(ix), fr(ix), 'k');
    set(h, 'lineWidth', 3);
    ylim([-3, 0]);
    title(sprintf('passband gain over output rate, decim=%i', decimationFactor));
    xlabel('frequency relative to output rate');
    ylabel('dB'); grid on;
    
    % *************************************************    
    % plot alias response
    % *************************************************    
    fb = FFT_frequencyBasis(numel(sAlias), 1);
    fr = 20*log10(abs(fft(sAlias) + eps));
    ix = find(fb > 0);
    
    figure(4); clf;
    h = semilogx(fb(ix), fr(ix), 'k');
    set(h, 'lineWidth', 3);
    %    ylim([-80, -20]);
    title(sprintf('alias response over output rate, decim=%i', decimationFactor));
    xlabel('frequency relative to output rate');
    ylabel('dB'); grid on;    
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

% ************************************
% 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('f / Hz');
        ylabel('dB');
        leg{end+1} = pp.legtext;
    end
    legend(leg);
    title(t);
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

PSD estimation window based-scilab code

Senthilkumar ● December 28, 2011 ● Coded in Scilab 
//Caption:Determination of spectrum of a signal
//With maximum normalized frequency f = 0.1
//using Rectangular window and Blackmann window
clear all;
close;
clc;
N = 61;
cfreq = [0.1 0];
[wft,wfm,fr]=wfir('lp',N,cfreq,'re',0);
disp(wft,'Time domain filter coefficients hd(n)=');
disp(wfm,'Frequency domain filter values Hd(w)=');
WFM_dB = 20*log10(wfm);//Frequency response in dB
for n = 1:N
 h_balckmann(n)=0.42-0.5*cos(2*%pi*n/(N-1))+0.08*cos(4*%pi*n/(N-1));
end
wft_blmn = wft'.*h_balckmann;
disp(wft_blmn,'Blackmann window based Filter output h(n)=')
wfm_blmn = frmag(wft_blmn,length(fr));
WFM_blmn_dB =20*log10(wfm_blmn);
subplot(2,1,1)
plot2d(fr,WFM_dB)
xgrid(1)
xtitle('Power Spectrum with Rectangular window Filtered M = 61','Frequency in cycles per samples  f','Energy density in dB')
subplot(2,1,2)
plot2d(fr,WFM_blmn_dB)
xgrid(1)
xtitle('Power Spectrum with Blackmann window Filtered  M = 61','Frequency in cycles per samples  f','Energy density in dB')