typedef signed long s32;
#define ABS(x) ((x) < 0 ? -(x) : (x))
#define CLAMP(x,min,max) {if ((x) <= (min)) (x) = (min); else if ((x) >= (max)) (x) = (max);}
#define KLPF_MAX_ADAPT 232L
#define KLPF_MIN_ADAPT 0L
#define KLPF_DEFAULT_ADAPT 160L
#define ACCEL_MIN_ADAPT 5
#define ACCEL_MAX_ADAPT 15
s32 gnKLpf = KLPF_DEFAULT_ADAPT; // IIR filter coefficient, 0 to 255
s32 gnKLpfMax = KLPF_MAX_ADAPT;
s32 gnAccel = 0; // rate of change of climbrate
s32 gnCpsx256 = 0; // climbrate in centimeters per second, with 8 bit fractional precision
s32 gnCps; // climbrate in centimeters per second
// IIR low pass filter using fixed point integer arithmetic. 8bits of fractional precision for IIR coefficient K.
// So the actual floating point coefficient is k = K/256, 0.0 <= k < 1.0. The filter computes
// climbRateFiltered = climbRateFiltered*k + newClimbRate*(1.0 - k). So K values close to 256 will result in
// heavy damping, K values close to 0 will result in almost no damping. The IIR low pass filtered output gnCps
// is the rounded up integer value, the 8bit fraction precision value gnCpsx256 is a global variable, so is
// maintained between calls to the filter. No integer divides required.
void sns_LPFClimbRate(void) {
s32 tmp;
tmp = gnCpsx256*gnKLpf + (gnCps<<8)*(256L-gnKLpf);
gnCpsx256 = (tmp >= 0 ? ((tmp + 128L)>>8) : ((tmp - 128L)>>8));
gnCps = (gnCpsx256>>8);
}
// Adapt the IIR filter coeffient to the rate of change. If samples are not changing much, increase the damping.
// If the samples are changing by a large amount, reduce the damping. This is done to reduce the response delay to a
// a step change, while still being able to heavily damp out jitter in steady state noisy samples.
// This function basically linearly interpolates KLpf between KLPF_MAX and KLPF_MIN based on the acceleration.
// Values of acceleration outside experimentally determined limits are clamped to the limits, depends on the
// application. A variable is used for KLpfMax to allow user to adjust the maximum allowed damping.
s32 sns_AdaptKLpf(void) {
s32 klpf,nAccel;
nAccel = ABS(gnAccel);
CLAMP(nAccel, ACCEL_MIN_ADAPT, ACCEL_MAX_ADAPT);
klpf = gnKLpfMax + ((KLPF_MIN_ADAPT-gnKLpfMax)*(nAccel - ACCEL_MIN_ADAPT))/(ACCEL_MAX_ADAPT-ACCEL_MIN_ADAPT);
CLAMP(klpf, KLPF_MIN_ADAPT, KLPF_MAX_ADAPT);
return klpf;
}
// usage :
// For each new data sample that arrives
// 1. calculate new unfiltered climbrate nCps from sample data buffer
// 2. gnAccel = (nCps - gnCps);
// 3. (optionally low pass filter gnAccel using same type of IIR filter, if required)
// 4. gnCps = nCps;
// 5. gnKLpf = sns_AdaptKLpf();
// 6. sns_LPFClimbRate();
% 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'])
clear all; close all;
test = 1; %see comments below
n = 10000; %number of test samples
switch test
case 1, z = linspace(.02,-.01,n); %gentle zero crossing
case 2, z = linspace(.0005,-.0005,n); %excess dc at zero crossing
case 3, z = randsrc(1,n,[.25 -.25]); %random sweep, creates messy signal
case 4, z = randsrc(1,n,[-.25:.001,.25]); %random sweep, creates messy signal
case 5, z = ones(1,n/4)*.01; z = [z -z z -z]; %discontinuity at zero crossing
end
%%%%%%%%%%%%%%%%%%%%%%%%% Finite NCO model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = 2^15 - 1; %max amplitude
M = 2^12; %NCO accumulator phase resolution
inc = round(M*z); %NCO accumulator phase increment
k = 2^12; %lut phase resolution (lut size),
lsb = log2(M) - log2(k); %LSBs discarded when addressing
lut = round(A*exp(j*2*pi*(0:k-1)/k)); %lut, one cycle cos/sin data
ptr = 0;
addr = 0;
for i = 1:n
y(i) = lut(addr+1); %add 1 for matlab LUT
ptr = mod(ptr + inc(i), M); %phase accumulator(ramp)
addr = round(ptr/2^lsb); %discard LSBs
addr(addr >= k) = addr - k; %check address overflow
end
%display output
plot(real(y),'-');hold
plot(imag(y),'r-');
//Program to generate different window functions
//For FIR Filter design based on windowing techniques
clear all;
close;
clc
M =11 ;
for n = 1:M
h_Rect(n) = 1;
h_hann(n) = 0.5-0.5*cos(2*%pi*(n-1)/(M-1));
h_hamm(n) = 0.54-0.46*cos(2*%pi*(n-1)/(M-1));
h_balckmann(n) = 0.42-0.5*cos(2*%pi*n/(M-1))+0.08*cos(4*%pi*n/(M-1));
end
plot2d(1:M,[h_Rect,h_hann,h_hamm,h_balckmann],[2,5,7,9]);
legend(['Rectangular Window';'Hanning';'Hamming';'Balckmann']);
title('Window Functions for Length M = 11')
//Band Pass FIlter of length M = 16
//Lower Cutoff frequency fp = 0.2 and Upper Cutoff frequency fs = 0.3
// Choose the number of cosine functions and create a dense grid
// in [0,0.1) and [0.2,0.35] and [0.425,0.5]
//magnitude for pass band = 1 & stop band = 0 (i.e) [0 1 0]
//Weighting function =[10 1 10]
clear all;
clc;
close;
hn = 0;
hm = 0;
hn=eqfir(16,[0 .1;.2 .35;.425 .5],[0 1 0],[10 1 10]);
[hm,fr]=frmag(hn,256);
disp(hn,'The Filter Coefficients are:')
figure
plot(.5*(0:255)/256,20*log10(frmag(hn,256)));
a = gca();
xlabel('Normalized Digital Frequency fr');
ylabel('Magnitude in dB');
title('Frequency Response of FIR BPF using REMEZ algorithm M=16')
xgrid(2)
function thr = oracleshrink1(CD,T)
%function used to calculate threshold for oracleshrink method
CD = CD(:)';
n = length(CD);
sx2 = (T-CD).^2;
b = cumsum(sx2);
[risk,best] = min(b);
hr = sqrt(sx2(best));
function [T,sigma1,sigma2] = model3_soft(CD,CH,CH1)
%function used to calculate the threshold and variance using modified
%bivariate model
%CD and CH : obtained after first stage wavelet decomposition
%CH1: obtained after second stage 2D-wavelet decomposition
sigman = (median(median(abs(CD))))/0.6745;
sigmay1 = 0;
sigmay2 = 0;
[m,n] = size(CD);
[m1,n1] = size(CH1);
%To calculate sigma1 value
for i =1:m
for j = 1:n
sigmay1 = sigmay1+((CH(i,j))^2);
end
end
sigmay1 = sqrt(2)*sigmay1/(m*n);
sigma1 = sqrt(max((((sigmay1))-((sigman)^2)),0));
if sigma1~=0
T = sqrt(3)*(sigman^2);
else
T = max(max(abs(CH)));
end
%To calculate sigma2 value
for i =1:m1
for j = 1:n1
sigmay2 = sigmay2+((CH1(i,j))^2);
end
end
sigmay2 = sqrt(2)*sigmay2/(m1*n1);
sigma2 = sqrt(max((((sigmay2))-((sigman)^2)),0));
function T = bishrink2(CD,CY)
sigman = (median(median(abs(CD))))/0.6745;
[m,n] = size(CY);
sigmay = 0;
for i =1:m
for j = 1:n
sigmay = sigmay+((CY(i,j))^2);
end
end
sigmay = sqrt(2)*sigmay/(m*n);
sigma = sqrt(max((((sigmay))-((sigman)^2)),0)); % Variance calculation
if sigma~=0
T = sqrt(3)*(sigman^2)/sigma; %Threshold
else
T = max(max(abs(CY)));
end
function [w1] = bishrink1(y1,y2,T)
% Bivariate Shrinkage Function
% y1 - a noisy coefficient value
% y2 - the corresponding parent value
% T - threshold value
% w1 - the denoised coefficient
R = sqrt(abs(y1).^2 + abs(y2).^2);
R = R - T;
R = R .* (R > 0);
w1 = y1 .* R./(R+T);
function [b,a] = allpass(order, Fc, Fs, Q)
% Returns allpass filter coefficients.
% Currently only 1st and 2nd orders are supported.
%
% Usage: [b,a] = ALLPASS(N, FC, FS, Q);
%
% N is the order of the allpass
% FC is the frequency a the 90deg phase shift point
% FS is the sampling rate
% Q is quality factor describing the slope of phase shift
%
% Author: sparafucile17 01/07/2004
if(order > 2)
error('Only 1st and 2nd orders are supported');
end
%Bilinear transform
g = tan(pi*(Fc/Fs));
if(order == 2)
d = 1/Q;
K = 1/(1 + d*g + g^2);
b0 = (1 - g*d + g^2) * K;
b1 = 2 * (g^2 - 1) * K;
b2 = 1;
a1 = b1;
a2 = b0;
else
b0 = (g-1)/(g+1);
b1 = 1;
b2 = 0;
a1 = b0;
a2 = 0;
end
b = [b0 b1 b2];
a = [1 a1 a2];
%In this example we are going to apply a low-pass filter to all WAV files,
%but it could be a multitude of different "processing" methods. This is
%only used to illustrate the batch processing example.
Fs = 44100; %Hz
Fc = 1000; %Hz
[b,a] = butter(2, Fc/(Fs/2), 'low');
%These are the input and output directories relative to this M-file
input_directory = 'test_input_database\';
output_direcory = 'test_output\';
%What extensions are you looking for? In our case, we want to batch
%process audio files that are store in the uncompressed *.WAV format
extension = '*.wav';
%Get the files
files=dir([input_directory '*.wav']);
N_files=numel(files);
%Loop through one file at a time
for i=1:N_files
%This is a stereo file and for this example,
%I only care about channel number 1
x = wavread([input_directory files(i).name ]);
x = x(:,1);
%Process the data by applying our filter
y = filter(b,a,x);
%Output the data as a unique filename based on input
wavwrite(y, Fs, [output_directory 'processed' files(i).name]);
disp(['Processed File: ' input_directory files(i).name 'as: ' ...
output_directory 'processed' files(i).name]);
end
