Reed Solomon Codes
function[n,p,r] = ReedSolomon_Codes(m)
// Reed-Solomon Codes
//Single-error-correcting RS code with a 2-bit byte
//m = m-bit symbol
k = 1^m; //number of message bits
t =1; //single bit error correction
n = 2^m-1; //code word length in 2-bit byte
p = n-k; //parity bits length in 2-bit byte
r = k/n; //code rate
disp(n,'code word length in 2-bit byte n =')
disp(p,'parity bits length in 2-bit byte n-k=')
disp(r,'Code rate:r = k/n =')
disp(2*t,'It can correct any error upto =')
endfunction
//Result
//-->exec('C:\Users\SENTHILKUMAR\Desktop\Communication_Toolbox\Digital_Communication\ReedSolomon_Codes.sci', -1)
//
//-->m=2
// m =
//
// 2.
//
//-->[n,p,r] = ReedSolomon_Codes(m)
//
// code word length in 2-bit byte n =
//
// 3.
//
// parity bits length in 2-bit byte n-k=
//
// 2.
//
// Code rate:r = k/n =
//
// 0.3333333
//
// It can correct any error upto =
//
// 2.
// r =
//
// 0.3333333
// p =
//
// 2.
// n =
//
// 3.
//Resampling filter performance
%%%%%%%%%%%%%%%%%% 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');Inplace matrix transposition
/* ------------------------------------------------------------------------- */
//
// This C++ code intents to do inplace transposition for matrix.
// It is more efficient for huge matrix, and quite generic as it uses classes
// and templates.
//
// Feel free to use it and to modify it.
// Allocation can be improved according to your system.
// Compiles with g++
/* ------------------------------------------------------------------------- */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <alloca.h>
#include <iostream>
using std::swap;
using std::cout;
template <typename T>
class matrix_t
{
public:
unsigned int row;
unsigned int col;
T *addr;
/* ------------------------------------------------------------------------- */
// FUNCTION: matrix_t
/* ------------------------------------------------------------------------- */
matrix_t()
{
row = 0;
col = 0;
addr = NULL;
}
/* ------------------------------------------------------------------------- */
// FUNCTION: init
//
// PARAM: _row
// PARAM: _col
//
// Init with allocation
//
// RETURN:
/* ------------------------------------------------------------------------- */
int init(unsigned int _row, unsigned int _col)
{
row = _row;
col = _col;
addr = (T *) calloc(sizeof(T), row * col);
if (addr == NULL)
{
return -1;
}
return 0;
}
/* ------------------------------------------------------------------------- */
// FUNCTION: term
//
// Terminate
//
// RETURN:
/* ------------------------------------------------------------------------- */
int term()
{
if (addr)
{
free(addr);
}
addr = NULL;
return 0;
}
/* ------------------------------------------------------------------------- */
// FUNCTION: fill
//
// Fill with values
/* ------------------------------------------------------------------------- */
void fill()
{
unsigned int x = 0;
int c = 0;
for (x = 0; x < row; x++)
{
unsigned int y = 0;
for (y = 0; y < col; y++)
{
addr[x * col + y] = (T) c++;
}
}
}
/* ------------------------------------------------------------------------- */
// FUNCTION: print
//
// Print to standard output
/* ------------------------------------------------------------------------- */
void print()
{
unsigned int x = 0;
cout << "row: " << row << "\n";
cout << "col: " << col << "\n";
for (x = 0; x < row; x++)
{
unsigned int y = 0;
for (y = 0; y < col; y++)
{
cout.width(3);
cout << addr[x * col + y] << " ";
}
cout << "\n";
}
cout << "\n";
cout << "\n";
}
/* ------------------------------------------------------------------------- */
// FUNCTION: operator==
//
// PARAM: mat2
//
// Comparaison of two matrices
//
// RETURN: 1 if equals; 0 if differents
/* ------------------------------------------------------------------------- */
int operator==(matrix_t mat2)
{
unsigned int x = 0;
if (row != mat2.row)
{
fprintf(stderr, "Error : different number of rows (%d != %d)\n", row, mat2.row);
return 0;
}
if (col != mat2.col)
{
fprintf(stderr, "Error : different number of columns (%d != %d)\n", col, mat2.col);
return 0;
}
for (x = 0; x < row; x++)
{
unsigned int y = 0;
for (y = 0; y < col; y++)
{
unsigned int pos = x * col + y;
if (addr[pos] != mat2.addr[pos])
{
return 0;
}
}
}
return 1;
}
/* ------------------------------------------------------------------------- */
// FUNCTION: transpose
//
// Transposition
/* ------------------------------------------------------------------------- */
int transpose()
{
unsigned int x = 0;
T *addr2 = (T *) calloc(sizeof(T), row * col);
if (addr2 == NULL)
{
return -1;
}
for (x = 0; x < row; x++)
{
unsigned int y = 0;
for (y = 0; y < col; y++)
{
int s = x * col + y;
int d = y * row + x;
addr2[d] = addr[s];
}
}
free(addr);
addr = addr2;
swap<unsigned int>(row, col);
return 0;
}
#define BITSET_POS(ptr,x) ((int)(((unsigned char*)(ptr))[(x)/8]))
#define BITSET_IS(ptr,x) ((int)(((unsigned char*)(ptr))[(x)/8]) & (int)(1<<((x)&7)))
#define BITSET_SET(ptr,x) ((unsigned char*)(ptr))[(x)/8] |= (unsigned char)(1<<((x)&7))
/* ------------------------------------------------------------------------- */
// FUNCTION: transpose_inplace
//
// Inplace transposition
/* ------------------------------------------------------------------------- */
int transpose_inplace()
{
unsigned int len = row * col;
unsigned int size = (len/8) + 1;
unsigned int i = 0;
unsigned char *data = (unsigned char*) calloc(1, size);
if (addr == NULL)
{
return -1;
}
while (i < len)
{
if (BITSET_POS(data, i) == 0xff)
{
i += 8;
continue;
}
if (BITSET_IS(data, i))
{
i++;
continue;
}
unsigned int begin = i;
BITSET_SET(data, i);
T tmp = (T)0;
do
{
swap<T>(addr[i], tmp);
unsigned int _row = i / col;
unsigned int _col = i % col;
i = _col * row + _row;
BITSET_SET(data, i);
}
while (i != begin);
swap<T>(addr[i], tmp);
i = begin + 1;
}
free(data);
swap<unsigned int>(row, col);
return 0;
}
};
/* ------------------------------------------------------------------------- */
// FUNCTION: main
//
// Test case :
// transpose a 9x23 matrix
// transpose inplace a 9x23 matrix
// and compare
//
/* ------------------------------------------------------------------------- */
int main()
{
matrix_t<double> m1, m1_ref;
m1.init(9, 23);
m1_ref.init(9, 23);
m1.fill();
m1_ref.fill();
m1.print();
m1.transpose_inplace();
m1_ref.transpose();
m1.print();
if (m1 == m1_ref)
{
printf("==\n");
}
else
{
printf("!=\n");
}
m1.term();
m1_ref.term();
return 0;
}Full band Echo canceller
//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')Circular Convolution
%circular convolution
%for testing you may use:
h = fir1(20,.3);
x = randn(1,1024);
%function y = conv_circ(h,x)
y = conv(h,x);
L1 = length(h);
L2 = length(x);
%add end to start, add start to end
temp = y(1:L1-1);
y(1:L1-1) = y(1:L1-1) + y(L2+(1:L1-1));
y(L2+(1:L1-1)) = y(L2+(1:L1-1)) + temp;
%compare to direct convolution
y2 = conv(h,x);
plot(y,'o-');hold
plot(y2,'r.--')
legend('circular','direct')passive Sum difference method for finding bearing
clc
clear all
close all
fs=20e7;
T=1/fs;
t=0:T:0.0000002;
f=10e6;
scanagle_deg=60; % Range covered by seeker antenna
step_deg=6;
% for ploting original arrays of sum, diff and ratio
theta=-deg2rad(scanagle_deg):deg2rad(step_deg):deg2rad(scanagle_deg);
for i=1:length(theta);
ant1=1*sin(2*pi*f*t); %Antenna 1
ant2=1*sin(2*pi*f*t+theta(i)); %Antenna 2
[my,mx]=max(ant1);
sum_ant(i)=ant1(mx)+ant2(mx); %Sum of antennas
diff_ant(i)=ant1(mx)-ant2(mx); %diff of antennas
ratio_ant(i)=diff_ant(i)/sum_ant(i); %ratio
end
% subplot(311)
% plot(t,ant1,t,ant2)
% subplot(211)
% plot(rad2deg(theta),sum_ant,rad2deg(theta),diff_ant)
% subplot(212)
% plot(rad2deg(theta),ratio_ant)
%
% figure
% subplot(211)
% plot(rad2deg(theta),sumn_ant,rad2deg(theta),diffn_ant)
% subplot(212)
% plot(rad2deg(theta),ration_ant)
%%%%%%%%%%%%%%%%Recieved signal at antenna for DOA%%%%%%%%%%%%%%%%%%%%%%%%%%
A_ant1=2; % Amplitude of antenna 1
A_ant2=3; % Amplitude of antenna 2
DOA_angle_target=20; % DOA of TARGET
sim_ant11=A_ant1*sin(2*pi*f*t); % Simulated antenna 1
sim_ant22=A_ant2*sin(2*pi*f*t+deg2rad(DOA_angle_target)); % Simulated antenna 2
sim_ant1=sim_ant11/max(sim_ant11); %normalization
sim_ant2=sim_ant22/max(sim_ant22);
[smy,smx]=max(sim_ant1);
%%%%also take for same sample value of data
sim_sum_ant=sim_ant1(smx)+sim_ant2(smx);
sim_diff_ant=sim_ant1(smx)-sim_ant2(smx);
sim_ratio_ant=sim_diff_ant/sim_sum_ant;
%%%%%%%%%%%%% Polynomial fitting of ratio_ant of 6th order%%%%%%%%%%%%
% Error in DOA is obtained because of this
% if polynomial fitting is removed or more accurate results are obtained
% DOA of target can be found with greater accuracy
% Polynomial was obtained though curve fitting toolbox
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p1 = 1.047e-008;
p2 = -6.907e-007;
p3 = 2.383e-005;
p4 = -0.0004914;
p5 = 0.008555;
p6 = -0.09626;
p7 = 0.4215;
x=0;
for ii=1:2200
x=x+0.01;
fitted_model_rat_ant(ii) = p1*x^6 + p2*x^5 + p3*x^4 + p4*x^3 + p5*x^2 + p6*x + p7;
end
theta_new=-scanagle_deg:.0545:scanagle_deg; % Theta for ploting fitted model
figure
plot(theta_new(1:2200),fitted_model_rat_ant)
c1=ceil(sim_ratio_ant*7000); % For comparison of sim_ratio ant and model
c2=ceil(fitted_model_rat_ant*7000); % Different threshold Threshold can be chosen
[r,c,v]= find(c2==c1);
detected_theta=(c.*0.0545)-60 %theta from curve fitting model
if(A_ant1>A_ant2) % condition for checking which angle was correct
correct_theta=detected_theta(1)
elseif(A_ant1<A_ant2)
correct_theta=detected_theta(2)
else
correct_theta=detected_theta
endvectors_intr.asm - External Interruption Configuration
*Vectors_intr.asm Vector file for interrupt INT11
.global _vectors ;global symbols
.global _c_int00
.global _vector1
.global _vector2
.global _vector3
.global _c_int04 ; símbolo para EXT_INT4
.global _vector5
.global _vector6
.global _vector7
.global _vector8
.global _vector9
.global _vector10
.global _c_int11 ;for INT11
.global _vector12
.global _vector13
.global _vector14
.global _vector15
.ref _c_int00 ;entry address
VEC_ENTRY .macro addr ;macro for ISR
STW B0,*--B15
MVKL addr,B0
MVKH addr,B0
B B0
LDW *B15++,B0
NOP 2
NOP
NOP
.endm
_vec_nmi:
B NRP
NOP 5
_vec_dummy:
B IRP
NOP 5
.sect ".vecs" ;aligned IST section
.align 1024
_vectors:
_vector0: VEC_ENTRY _c_int00 ;RESET
_vector1: VEC_ENTRY _vec_nmi ;NMI
_vector2: VEC_ENTRY _vec_dummy ;RSVD
_vector3: VEC_ENTRY _vec_dummy
_vector4: VEC_ENTRY _c_int04 ;INT04 Externa
_vector5: VEC_ENTRY _vec_dummy
_vector6: VEC_ENTRY _vec_dummy
_vector7: VEC_ENTRY _vec_dummy
_vector8: VEC_ENTRY _vec_dummy
_vector9: VEC_ENTRY _vec_dummy
_vector10: VEC_ENTRY _vec_dummy
_vector11: VEC_ENTRY _c_int11 ;ISR address
_vector12: VEC_ENTRY _vec_dummy
_vector13: VEC_ENTRY _vec_dummy
_vector14: VEC_ENTRY _vec_dummy
_vector15: VEC_ENTRY _vec_dummyTI F281x SCI Configuration
/**********************************************************************/
#include "DSP281x_Device.h"
/**********************************************************************
* Function: InitSci()
*
* Description: This function initializes F281x SCI. In this case the desired rate
* is 9600 bps, 8 databits, 1 stop bit and no parity.
**********************************************************************/
void InitSci(void)
{
/*** SCI reset */
ScibRegs.SCICTL1.bit.SWRESET= 0;
/*** Setting the baud rate to 9600 bps assuming LSPCLK=150 Mhz */
ScibRegs.SCIHBAUD = 0x07;
ScibRegs.SCILBAUD = 0xA0;
/*** Configuration for no stop bits, no parity and 8 data bits */
ScibRegs.SCICCR.all = 0x07;
// bit 7 0: STOP BITS
// bit 6 0: PARITY
// bit 5 0: PARITY ENABLE
// bit 4 0: LOOP BACK ENA
// bit 3 0: ADDR/IDLE MODE
// bit 2-0 111: SCI CHAR
/*** Transmitter and receiver enabled */
ScibRegs.SCICTL1.all = 0x03;
// bit 7 0: Reserved
// bit 6 0: RX ERR INT
// bit 5 0: SW RESET
// bit 4 0 Reserved
// bit 3 0: TXWAKE
// bit 2 0: SLEEP
// bit 1 1: TXENA, Habilitación del transmisor
// bit 0 1: RXENA, Habilitación del receptor
/*** Tx and Rx interruption enabled */
ScibRegs.SCICTL2.all = 0x0003;
// bit 7 0: TXRDY
// bit 6 0: TX EMPTY
// bit 5-2 0's: Reserved
// bit 1 1: RX/BK INT
// bit 0 1: TX INT
/***Emulation mode configuration */
ScibRegs.SCIPRI.all = 0x0010;
// bit 7-5 0's: Reserved
// bit 4-3 10: Emulation Mode
// bit 2-0 0's: Reserved
/*** SCI Reset */
ScibRegs.SCICTL1.bit.SWRESET= 1;
/*** SCI interruptions enabled in PIE */
PieCtrlRegs.PIEIER9.bit.INTx4= 1;
PieCtrlRegs.PIEIER9.bit.INTx3= 1;
IER |= 0x0100;
}
/**End of file*/CRC calculation
unsigned char crc(unsigned char data[],int L)
{
unsigned char P=0x83;
unsigned char *data_ptr;
unsigned char crc_result;
int bit;
crc_result=0;
for(data_ptr=data;data_ptr<data+L;data_ptr++)
{
for(bit=7;bit>=0;bit--)
{
crc_result=(crc_result<<1)+((*data_ptr&0x1<<bit)>>bit);
if (crc_result>=128)
{
crc_result=crc_result^P;
}
}
}
//multiply by 2^n-k;
for(bit=6;bit>=0;bit--)
{
crc_result=(crc_result<<1);
if (crc_result>=128)
{
crc_result=crc_result^P;
}
}
return crc_result;
}Image denoise: denoising using soft threshold
function y = soft(x,T)
%function used to denoise a noisy image with given soft threshold
%x = noisy image
%T = threshold
y = max(abs(x) - T, 0);
y = y./(y+T) .* x;Interleave Zeros
function out = interleave_zeros(a)
% Interleave vector with zeros.
% Only 1D matrices are supported.
%
% Usage: b = INTERLEAVE_ZEROS(a);
%
% A is the vector that will have zeros inserted
% B will contain the vector set A and the interleaved
% zeros (twice the length)
%
% Author: sparafucile17 2/13/04
my_zeros = zeros(1, length(a));
out = matintrlv([a my_zeros],2,length(a));Zero Crossing Counter
function count = zero_crossings(x)
% Count the number of zero-crossings from the supplied time-domain
% input vector. A simple method is applied here that can be easily
% ported to a real-time system that would minimize the number of
% if-else conditionals.
%
% Usage: COUNT = zero_crossings(X);
%
% X is the input time-domain signal (one dimensional)
% COUNT is the amount of zero-crossings in the input signal
%
% Author: sparafucile17 06/27/04
%
% initial value
count = 0;
% error checks
if(length(x) == 1)
error('ERROR: input signal must have more than one element');
end
if((size(x, 2) ~= 1) && (size(x, 1) ~= 1))
error('ERROR: Input must be one-dimensional');
end
% force signal to be a vector oriented in the same direction
x = x(:);
num_samples = length(x);
for i=2:num_samples
% Any time you multiply to adjacent values that have a sign difference
% the result will always be negative. When the signs are identical,
% the product will always be positive.
if((x(i) * x(i-1)) < 0)
count = count + 1;
end
endML Estimation of Several Frequencies Using the Nonuniform FFT
function st = Demo
% Author: J. Selva.
% Date: July 2011.
% Version: 0.
st = [];
n = input('Options: 1 perform single estimation, 2 estimate RMS error: ');
switch n
case 1
st = TestSingleTrial;
case 2
TestRMSError;
end
function st = TestSingleTrial
%Scenario with 5 frequencies.
M = 2048;
SNRInSpecdB = 27; %Signal-to-noise ratio after correlation.
Nf = 5;
df = 1/M; %Minimum separation among frequencies.
MaxAmp = 1; %Maximum amplitude.
PFA = 10.^-5;
nVar = M*MaxAmp^2/10^(SNRInSpecdB/10);
nThres = NoiseThres(M,nVar,PFA);
MinAmp = 2*nThres/M;
[vf,va] = RandomFreqAndAmpVectors(Nf,df,MinAmp,MaxAmp);
y = DataRealization(vf,va,M,nVar);
p = Preprocessing(y,nVar,PFA);
disp('True frequencies: ')
disp(vf)
disp('True amplitudes: ')
disp(va)
disp(['Post-correlation Signal-to-noise ratio: ',...
num2str(SNRInSpecdB),' dB'])
s = input('To compute ML estimate, use (S)erial or (M)anual method? ','s');
if strcmp(lower(s),'s')
st = MethodSerial1(p);
else
st = MethodManual1(p);
end
disp('ML frequency estimates: ')
disp(st.Freqs)
disp('ML amplitude estimates: ')
disp(st.Amps)
function TestRMSError
%Estimates the RMS error and computes the Cramer-Rao bounds.
format compact
format short g
M = 2048;
SNRInSpecdB = 27;
Nf = 5;
df = 1/M;
MaxAmp = 1;
PFA = 10.^-5;
nVar = M*MaxAmp^2/10^(SNRInSpecdB/10);
nThres = NoiseThres(M,nVar,PFA);
MinAmp = 2*nThres/M;
[vf,va] = RandomFreqAndAmpVectors(Nf,df,MinAmp,MaxAmp);
crb = CRBBound(vf,va,M,nVar);
Ntr = 0;
S2 = zeros(Nf,1);
for k = 1:1e5
y = DataRealization(vf,va,M,nVar);
p = Preprocessing(y,nVar,PFA);
st = MethodSerial1(p);
if length(st.Freqs) == 5
Ntr = Ntr + 1;
S2 = S2 + (vf-st.Freqs).^2;
end
if mod(k,100) == 0
disp(['Number or trials: ', num2str(k)])
disp([' Frequencies: ' num2str(vf.')])
disp(['Deviation CRB bounds (dB): ' num2str(20*log10(crb.'))])
disp([' RMS errors (dB): ' num2str(10*log10(S2(:).'/Ntr))])
disp(' ');
end
end
function [vf,va] = RandomFreqAndAmpVectors(Nf,MinFreqSep,MinAmp,MaxAmp)
%Constructs random vectors of frequencies and amplitudes. The module of the largest
%amplitude is MaxAmp. There are no amplitudes with module below MinAmp, and the minimum
%separation among frequencies is MinFreqSep.
vf = sort(rand(Nf,1));
vf = vf + MinFreqSep*(cumsum(ones(Nf,1))-1);
vf = vf/(1+(Nf-1)*MinFreqSep) - 0.5;
pr = rand(Nf,1);
ind = floor(rand(1)/Nf)+1;
pr(ind) = 1;
va = (MinAmp+(MaxAmp-MinAmp)*pr) .* exp(1i*2*pi*rand(Nf,1));
function y = DataRealization(vf,a,M,nVar)
m1 = FirstIndex(M);
y = exp(1i*2*pi*(m1:m1+M-1).'*vf(:).')*a(:)+randn(M,2)*[1 ; 1i]*sqrt(nVar/2);
function p = Preprocessing(y,nVar,PFA)
%Performs the preprocessing. Fundamentally this consists in computing the FFT inside
%"Correlation.m".
%The input parameters are
%
% y: Data vector.
% nVar: Noise variance estimate.
% PFA: Probability of seeing a frequency when there isn't any. Usually a very low value
% like 1e-5.
p = [];
p.M = length(y);
p.nVar = nVar;
p.PFA = PFA;
p.nThres = NoiseThres(p.M,p.nVar,p.PFA); %Noise threshold.
icf = Correlation(y,p.nThres); %Struct with data for interpolating the correlation.
p.cfd = Cost(real(y(:)'*y(:)),icf); %Struct with data for computing the ML cost function ...
%and its derivatives.
p.freqThres = 1/(10*p.M); %If two frequencies differ in less than ...
%freqThres in the iterative search, then one of them is
%discarded. This prevents ill-conditioning.
p.freqStopThres = 1/(p.M*1e4); %Threshold for stopping the ML cost function minimization.
function st = MethodSerial1(p)
%Adds frequencies and computes ML estimates sequentially.
vfEst = [];
ModifiedEstimation = logical(1);
while ModifiedEstimation
[vf1,va1] = Correlation(p.cfd.icf,'gridMainLocalMaxima',vfEst);
ModifiedEstimation = ~isempty(vf1);
if ModifiedEstimation
vfEst(end+1) = vf1(1);
[listv,listLv,listCL,ErrRatio] = LookForMinimum(p.cfd,vfEst,...
10.^-8,p.freqThres,p.freqStopThres);
vfEst = sort(listv(:,end));
end
end
st = Report(p,vfEst);
function st = MethodManual1(p,varargin)
%Allows the user to set the initial frequency estimates on a plot.
if nargin == 1
vfEst = [];
else
vfEst = varargin{1};
end
vfEst = vfEst(:);
figure
subplot(2,1,1)
Correlation(p.cfd.icf,'PlotPeriodogramFit',vfEst);
yl = get(gca,'YLim');
subplot(2,1,2)
Correlation(p.cfd.icf,'PlotPeriodogramFitError',vfEst);
while 1
pr1 = input('A/a -> Add freq., D/d -> Delete freq, X/x -> exit: ','s');
switch lower(pr1)
case 'a'
[f1,pr] = ginput(1);
vfEst = sort([vfEst; f1]);
[listv,listLv,listCL,ErrRatio] = LookForMinimum(p.cfd,vfEst,10.^-8,...
p.freqThres,p.freqStopThres);
vfEst = sort(listv(:,end));
subplot(2,1,1)
Correlation(p.cfd.icf,'PlotPeriodogramFit',vfEst);
subplot(2,1,2)
Correlation(p.cfd.icf,'PlotPeriodogramFitError',vfEst);
case 'd'
[f1,pr] = ginput(1);
[pr,ind] = min(abs(f1-vfEst));
vfEst(ind) = [];
subplot(2,1,1)
Correlation(p.cfd.icf,'PlotPeriodogramFit',vfEst);
subplot(2,1,2)
Correlation(p.cfd.icf,'PlotPeriodogramFitError',vfEst);
case 'x'
st = Report(p,vfEst);
return
end
end
function thres = NoiseThres(M,nVar,PFA)
%Computes the noise threshold. The method will discard any correlation local maximum below
%this threshold.
thres = sqrt((M*nVar/2)*icdf('chi2',(1-PFA)^(1/M),2));
function m1 = FirstIndex(M)
m1 = -ceil(M/2);
function varargout = Correlation(varargin)
%Computes the correlation value and its derivatives by calling "Barycentric". It does
%other things like plotting.
if ~isstruct(varargin{1})
st = [];
temp = {[],[],10.^-15,1.5};
[temp{1:nargin}] = deal(varargin{:});
[a,st.nThres,epsilon,OvFMin] = deal(temp{:});
a = a(:).';
st.M = length(a);
st.K = 2^nextpow2(st.M*OvFMin);
st.bary = Barycentric(0,SampleFourierCorrelation1DVer1(a,st.K),...
-2*FirstIndex(st.M),1/st.K,whatP(1/st.K,-2*FirstIndex(st.M),epsilon));
st.IndActive = find(abs(st.bary.a) > st.nThres);
st.IndActive = unique([st.IndActive-1,st.IndActive,st.IndActive+1]);
if st.IndActive(1) == 0
st.IndActive(1) = [];
end
if st.IndActive(end) == st.K+1
st.IndActive(end) = [];
end
st.aTActive = st.bary.a(st.IndActive);
st.vfActive = centreFreq((st.IndActive-1)/st.K);
varargout = {st};
return
end
[st,msg] = deal(varargin{1:2});
switch msg
case 'value'
varargout = {Barycentric(st.bary,varargin{3:end})};
case 'gridValueInRange'
[ind1,ind2] = deal(varargin{3:4});
varargout = {st.bary.a(mod(ind1:ind2,st.K)+1)};
case 'ampEstimates'
vf = varargin{3};
if length(vf) == 0
varargout = {[]};
end
varargout = {RegGeoSum(st.M,vf,vf,0,0) \ Barycentric(st.bary,vf)};
case 'gridValue'
lInd = varargin{3};
if nargin == 4
vf = varargin{4};
else
vf = [];
end
if isempty(vf)
varargout = {st.bary.a(mod(lInd,st.K)+1)};
else
vf1 = centreFreq(lInd/st.K);
VRaa = RegGeoSum(st.M,vf,vf,0,0);
VRa1a = RegGeoSum(st.M,vf1,vf,0,0);
pr = VRa1a*(VRaa\Barycentric(st.bary,vf));
pr = reshape(pr,size(lInd));
varargout = {st.bary.a(mod(lInd,st.K)+1)-pr};
end
case 'gridFit'
lInd = varargin{3};
if nargin == 4
vf = varargin{4};
else
vf = [];
end
if isempty(vf)
varargout = {zeros(size(lInd))};
else
vf1 = centreFreq(lInd/st.K);
VRaa = RegGeoSum(st.M,vf,vf,0,0);
VRa1a = RegGeoSum(st.M,vf1,vf,0,0);
pr = VRa1a*(VRaa\Barycentric(st.bary,vf));
pr = reshape(pr,size(lInd));
varargout = {pr};
end
case 'plotErrorIndB'
if nargin == 2
vf = [];
else
vf = varargin{3};
end
s = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1,vf)));
gf = centreFreq((0:st.K-1)/st.K);
[gf,ord] = sort(gf);
s = s(ord);
plot(gf,s,[-0.5,0.5],20*log10(abs(st.nThres)*[1,1]))
Ms = max(s);
set(gca,'YLim',[Ms-50,Ms+10])
xlabel('Frequency (f)')
ylabel('Error spectrum (dB)')
grid on
varargout = {};
case 'PlotPeriodogramFit'
vfEst = varargin{3};
if isempty(vfEst)
s0 = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1)));
gf = centreFreq((0:st.K-1)/st.K);
[gf,ord] = sort(gf);
s0 = s0(ord);
plot(gf,s0)
hold on
nt = 20*log10(abs(st.nThres));
plot([-0.5,0.5],nt([1,1]).',...
'LineWidth',2,'LineStyle','--','Color',[1 0 0])
hold off
text(-0.49,nt+2,'Noise threshold')
xlabel('Frequency (f)')
ylabel('Amplitude (dB)')
title('Periodogram data, fit, and est. frequencies (stems)')
return
end
vaEst = Correlation(st,'ampEstimates',vfEst);
stem(vfEst,20*log10(st.M*abs(vaEst)),'g')
hold on
s0 = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1)));
gf = centreFreq((0:st.K-1)/st.K);
[gf,ord] = sort(gf);
s0 = s0(ord);
plot(gf,s0,'b')
hold on
nt = 20*log10(abs(st.nThres));
plot([-0.5,0.5],nt([1,1]).',...
'LineWidth',2,'LineStyle','--','Color',[1 0 0])
s1 = 20*log10(abs(Correlation(st,'gridFit',0:st.K-1,vfEst)));
s1 = s1(ord);
plot(gf,s1,'r')
Ms = max(s0);
set(gca,'YLim',[Ms-50,Ms+10])
xlabel('Frequency (f)')
ylabel('Amplitude (dB)')
title('Periodogram fit')
hold off
varargout = {};
case 'PlotPeriodogramFitError'
if nargin == 2
vfEst = [];
else
vfEst = varargin{3};
end
vaEst = Correlation(st,'ampEstimates',vfEst);
s0 = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1,vfEst)));
gf = centreFreq((0:st.K-1)/st.K);
[gf,ord] = sort(gf);
s0 = s0(ord);
nt = 20*log10(abs(st.nThres));
plot(gf,s0)
hold on
plot([-0.5,0.5],nt([1,1]).',...
'LineWidth',2,'LineStyle','--','Color',[1 0 0])
hold off
Ms = max(s0);
set(gca,'YLim',[Ms-50,Ms+10])
xlabel('Frequency (f)')
ylabel('Amplitude (dB)')
title('Periodogram error')
varargout = {};
case 'gridValueGivenAmps'
[lInd,vf,aEst] = deal(varargin{3:5});
vf1 = centreFreq(lInd/st.K);
VRa1a = RegGeoSum(st.M,vf1,vf,0,0);
pr1 = st.bary.a(mod(lInd,st.K)+1);
pr1 = reshape(pr1,size(lInd));
pr2 = VRa1a*aEst;
pr2 = reshape(pr2,size(lInd));
varargout = {pr1-pr2};
case 'gridLocalMaxima'
if nargin == 2
pr = [-2,-1,0].';
pr = abs(st.bary.a(mod(st.IndActive([1;1;1],:) + ...
pr(:,ones(1,length(st.IndActive))),st.K)+1));
pr = pr(1,:)<=pr(2,:) & pr(3,:)<=pr(2,:);
vf = st.vfActive(pr);
va = st.aTActive(pr);
varargout = {vf,va};
return
end
vf = varargin{3};
if length(vf) == 0
[varargout{1:2}] = Correlation(st,'gridLocalMaxima');
return
end
amp = Correlation(st,'ampEstimates',vf);
tmp = abs(Correlation(st,...
'gridValueGivenAmps',st.IndActive-1,vf,amp)) >= st.nThres;
ind = st.IndActive(tmp);
pr = [-2,-1,0].';
ind1 = ind([1,1,1].',:) + pr(:,ones(1,length(ind)));
pr = abs(Correlation(st,'gridValueGivenAmps',ind1,vf,amp));
pr = pr(1,:)<=pr(2,:) & pr(3,:)<=pr(2,:);
ind = ind(pr)-1;
if isempty(ind)
varargout = {[],[]};
else
varargout = {centreFreq(ind/st.K),Correlation(st,...
'gridValueGivenAmps',ind,vf,amp)};
end
case 'gridMainLocalMaxima'
[vf,va] = Correlation(st,'gridLocalMaxima',varargin{3:end});
if length(vf) == 0
varargout = {[],[]};
return
end
[varargout{1:2}] = selectLocalMaximaVer3(vf,va,st.nThres,st.M);
case 'gridMaximum'
pr = Correlation(st,'gridValue',0:st.K-1,varargin{3:end});
[ma,ind] = max(abs(pr));
varargout = {(ind-1)/st.K,ma};
end
function tab = SampleFourierCorrelation1DVer1(a,K)
M = length(a);
m1 = FirstIndex(M);
a = a(:).';
tab = fft([a(-m1+1:end),zeros(1,K-M),a(1:-m1)]);
function out = centreFreq(f)
out = f + ceil(-0.5-f);
function P = whatP(T,B,epsilon)
P = ceil(acsch(epsilon)/(pi*(1-B*T)));
function d = dist(x,y)
if x<y
d = min([y-x,x+1-y]);
else
d = min([x-y,y+1-x]);
end
function v = sincMask(x,y,M)
v = zeros(size(x));
for k = 1:length(x)
x1 = max([dist(x(k),y)-1/(2*M),0]);
v(k) = 1/max([M*sin(pi*x1),1]);
end
function v = slBoundVer2(x,vam,vfe,nThres,M1)
vam = vam(:);
vfe = vfe(:);
sll = sincMask(vfe,x,M1).'*vam;
v = max([sll+nThres,1.1*sll]);
function varargout = selectLocalMaximaVer3(vf,va,nThres,M)
va1 = abs(va(:));
[pr,ord] = sort(va1);
ord = flipud(ord);
vf1 = vf(ord);
va1 = va1(ord);
n1 = length(va1);
if va1(1) < nThres
varargout = {[],[]};
return
end
va2 = va1(1);
va2Abs = va1(1);
vf2 = vf1(1);
for n=2:n1
if va1(n)<nThres
continue
end
if va1(n) < slBoundVer2(vf1(n),va2Abs,vf2,nThres,M)
continue
end
vf2(end+1) = vf1(n);
va2(end+1) = va1(n);
va2Abs(end+1) = abs(va(ord(n)));
end
varargout = {vf2,va2};
function out = Barycentric(varargin)
%Performs barycentric interpolation using the method in [Bary].
if ~isstruct(varargin{1})
st = [];
[st.n0,st.a,st.B,st.T,st.P] = deal(varargin{:});
st.N = length(st.a);
st.alpha = baryWeights(st.T,st.B,st.P);
st.alpha = st.alpha([1:st.P,st.P+2:end,st.P+1]);
st.vt = (-st.P:st.P)*st.T;
st.vt = st.vt([1:st.P,st.P+2:end,st.P+1]);
out = st;
else
[st,v] = deal(varargin{1:2});
if nargin == 3
nd = varargin{3};
else
nd = 0;
end
v = v(:);
Nv = numel(v);
n = floor(v/st.T+0.5);
u = v - n*st.T;
n1 = n - st.n0;
pr = [-st.P:-1,1:st.P,0];
da = st.a(mod(n1(:,ones(1,2*st.P+1)) + pr(ones(Nv,1),:),st.N)+1);
out = DerBarycentricInterp3Vec(st.alpha,da,st.vt,u,nd);
end
function w = baryWeights(T,B,P)
vt = (-P:P)*T;
g = APPulse(vt,1/T-B,T*(P+1));
gam = gamma(vt/T+P+1) .* gamma(-vt/T+P+1) .* g;
N = length(vt);
LD = ones(1,N);
for k = 1:N
LD(k) = prod(vt(k)-vt(1:k-1))* prod(vt(k)-vt(k+1:N));
end
w = gam ./ LD;
w = w / max(abs(w));
function v = APPulse(t,B,TSL)
v = real(sinc(B*sqrt(t.^2-TSL^2)))/real(sinc(1i*pi*B*TSL));
function out = DerBarycentricInterp3Vec(alpha,zL,t,tauL,nd)
NtauL = size(tauL,1);
LBtau = 200; %Block size for vectorized processing. Change this to adjust the memory
%usage.
NBlock = 1+floor(NtauL/LBtau);
t = t(:).';
Nt = length(t);
out = zeros(NtauL,nd+1);
for r = 0:NBlock-1
ind1 = 1+r*LBtau;
ind2 = min([(r+1)*LBtau,NtauL]);
Ntau = ind2-ind1+1;
z = zL(ind1:ind2,:);
tau = tauL(ind1:ind2);
vD = zeros(Ntau,1);
LF = [1,zeros(1,Nt-1)];
LF = LF(ones(Ntau,1),:);
for k = 1:Nt-1
vD = vD .* (tau-t(k))+alpha(k)*LF(:,k);
LF(:,k+1) = LF(:,k) .* (tau-t(k));
end
vD = vD .* (tau-t(Nt)) + alpha(Nt) * LF(:,Nt);
for kd = 0:nd
pr = z(:,end); z1 = z-pr(:,ones(1,Nt)); cN = zeros(Ntau,1);
for k = 1:Nt-1
cN = cN .* (tau-t(k)) + z1(:,k) .* alpha(k) .* LF(:,k);
end
cN = cN ./ vD;
ztau = z(:,end)+(tau-t(end)) .* cN;
out(ind1:ind2,kd+1) = ztau;
if kd < nd
pr1 = ztau;
pr1 = z(:,1:end-1) - pr1(:,ones(1,Nt-1));
pr2 = t(1:end-1);
pr2 = pr2(ones(Ntau,1),:)-tau(:,ones(1,Nt-1));
z = [pr1 ./ pr2, cN] * (kd+1);
end
end
end
function varargout = Cost(varargin)
%Computes the ML cost function differentials by calling "Correlation" and "RegGeoSum".
if ~isstruct(varargin{1})
st = [];
[st.EnergyRyy,st.icf] = deal(varargin{:});
varargout = {st};
return
end
st = varargin{1};
if nargin > 2
vf = varargin{3};
else
vf = [];
end
switch varargin{2}
case 'value'
vf = varargin{3};
Rss = RegGeoSum(st.icf.M,vf,vf,0,0);
Rsd = RegGeoSum(st.icf.M,vf,vf,0,1);
Rdd = RegGeoSum(st.icf.M,vf,vf,1,1);
pr = Correlation(st.icf,'value',vf,1);
Rsy = pr(:,1);
Rdy = pr(:,2);
pr = Rss\[Rsy,Rsd];
MRsy = pr(:,1);
MRsd = pr(:,2:end);
varargout = {st.EnergyRyy-RealTrace(Rsy',MRsy),-2*RealDiag(MRsy,Rdy'-MRsy'*Rsd),...
2*RealHadamard((Rdd'-Rsd'*MRsd).',MRsy*MRsy')};
case 'value0'
if length(vf) == 0
varargout = {st.EnergyRyy};
return
end
vf = varargin{3};
Rss = RegGeoSum(st.icf.M,vf,vf,0,0);
Rsy = Correlation(st.icf,'value',vf);
[varargout{1}] = st.EnergyRyy-RealTrace(Rsy',Rss\Rsy);
end
function v = RealHadamard(A,B)
v = real(A).*real(B)-imag(A).*imag(B);
function C = RealDiag(A,B)
% Same as real(diag(A*B)) but using less flops.
C = RealHadamard(A.',B);
if size(C,1) == 1
C = C.';
else
C = sum(C).';
end
function C = RealTrace(A,B)
% Same as real(trace(A*B)) but using less flops.
C=RealHadamard(A.',B);
C=sum(C(:));
function [listv,listLv,listCL,ErrRatio] = ...
LookForMinimum(L,v0,nThres,freqThres,freqStopThres,varargin)
%Executes a descent algorithm using a method similar to that in [ML].
if ~isempty(varargin)
nItMax = varargin{1};
else
nItMax = 20;
end
v0 = v0(:);
Nf = numel(v0);
listv = [v0,zeros(Nf,nItMax)];
listLv = [Cost(L,'value0',v0),zeros(1,nItMax)];
listCL = [NaN,zeros(1,nItMax)];
succ = 1;
nIt = 1;
while succ && nIt <= nItMax
[succ,v,Lv,CL] = NewtonIt(L,listv(:,nIt),freqThres);
if succ
nIt = nIt + 1;
listv(:,nIt) = v;
listLv(nIt) = Lv;
listCL(nIt) = CL;
else
break
end
if max(abs(listv(:,nIt)-listv(:,nIt-1))) < freqStopThres
break
end
end
ErrRatio = (Cost(L,'value0')-listLv(end))/listLv(end);
listv(:,nIt+1:end) = [];
listLv(nIt+1:end) = [];
listCL(nIt+1:end) = [];
function [v,vf] = validVector(vf,thres)
vf = sort(vf(:));
if length(vf)==1
v = 1;
return
end
m = vf(2:end)-vf(1:end-1) < thres;
if any(m)
v = 0;
vf([logical(0); m]) = [];
else
v = 1;
end
function [succ,v1,Lv1,nCutsLeft] = NewtonIt(L,v0,thres)
Lv1 = [];
nCutsLeft = [];
[vv,v1] = validVector(v0,thres);
if ~vv
succ = 0;
return
end
[Lv0,g0,h0] = Cost(L,'value',v0);
d = -h0\g0;
nCutsLeft = 20;
while nCutsLeft > 0
[vv,v1] = validVector(v0+d,thres);
if ~vv
succ = 0;
v1 = v0;
end
Lv1 = Cost(L,'value0',v0+d);
if Lv1 < Lv0
succ = 1;
v1 = v0 + d;
break;
end
nCutsLeft = nCutsLeft - 1;
d = d/2;
end
if nCutsLeft == 0
succ = 0;
end
function out = RegGeoSum(varargin)
%Computes the value of a geometric sum or its differentials using the exact formula. It
%takes care of the preventable singularity at zero.
switch nargin
case 2
[M,vf] = deal(varargin{:});
nd = 0;
msg = 'value';
case 3
[M,vf,nd] = deal(varargin{:});
msg = 'value';
case 5
[M,vf1,vf2,nd1,nd2] = deal(varargin{:});
msg = 'corr';
end
switch msg
case 'value'
out = RegGeoSum1(M,vf,nd);
case 'corr'
vf1 = vf1(:);
Nf1 = length(vf1);
vf2 = vf2(:);
Nf2 = length(vf2);
out = RegGeoSum1(M,-vf1(:,ones(1,Nf2))+vf2(:,ones(1,Nf1)).',nd1+nd2)*(-1)^nd1;
end
function v = RegGeoSum1(M,vf,nd)
m1 = -ceil(M/2);
ind = abs(vf)>eps^(1/3);
vf1 = vf(ind);
v = zeros(size(vf));
switch nd
case 0
v(ind) = -(exp(1i*2*pi*m1*vf1)-exp(1i*2*pi*(m1+M)*vf1))./(-1+exp(1i*2*pi*vf1));
v(~ind) = M;
case 1
v(ind) = (-2*1i*pi*( exp(1i*2*pi*vf1*(1 + m1))*(-1 + m1) - exp(1i*2*pi*vf1*m1)*m1 ...
-exp(1i*pi*2*vf1*(1 + M + m1))*(-1 + M + m1) + ...
exp(1i*2*pi*vf1*(M + m1))*(M + m1)))./(-1 + exp(1i*2*pi*vf1)).^2;
v(~ind) = pi*1i*M*(-1+M+2*m1);
case 2
v(ind) = (4*(exp(1i*2*pi*vf1*(2 + m1))*(-1 + m1)^2 + exp(1i*2*pi*vf1*m1)*m1^2 - ...
exp(1i*2*pi*vf1*(2 + M + m1))*(-1 + M + m1)^2 - ...
exp(1i*2*pi*vf1*(M + m1))*(M + m1)^2 + ...
exp(1i*2*pi*vf1*(1 + m1))*(1 - 2*(-1 + m1)*m1) + ...
exp(1i*2*pi*vf1*(1 + M + m1))* ...
(-1 + 2*M^2 + 2*(-1 + m1)*m1 + M*(-2 + 4*m1)))*pi^2)./ ...
(-1 + exp(1i*2*pi*vf1)).^3;
v(~ind) = -(2*pi^2/3)*M*(1+6*m1^2+6*m1*(-1+M)-3*M+2*M^2);
end
function st = Report(p,vfEst)
st = [];
st.Freqs = vfEst;
st.Amps = Correlation(p.cfd.icf,'ampEstimates',vfEst);
function FreqsDevBounds = CRBBound(vf,va,M,nVar)
%Computes the Cramer-Rao bound of the frequencies.
va = va(:);
vf = vf(:);
Raa = RegGeoSum(M,vf,vf,0,0);
Rad = RegGeoSum(M,vf,vf,0,1);
Rdd = RegGeoSum(M,vf,vf,1,1);
FreqsDevBounds = sqrt((nVar/2)*diag(inv(real((conj(va)*va.') .* (Rdd-Rad'*(Raa\Rad))))));
function st = MLEstimate(y,nVar,PFA)
p = Preprocessing(y,nVar,PFA);
st = MethodSerial1(p);
st.CRBBoundsOnFreqDevs = CRBBound(st.Freqs,st.Amps,length(y),nVar);Phaser audio effect
/*
Phaser audio effect:
Xin Yout
-------------------------------[dir_mix]--------------------->[+]--------->
| ^
| |
|-->[VNS1]-->[VNS2]-->[VNS3]...-->[VNSN]-->[pha_mix]------
^ ^ ^ ^
| | | |
|--------|--------|---...------
|
[LFO]
VNS = Variable notch stage
*/
#include "br_iir.h"
#include "Phaser.h"
/*This defines the phaser stages
that is the number of variable notch blocks
*/
#define PH_STAGES 20
static short center_freq; /*Center frequency counter*/
static short samp_freq; /*Sampling frequency*/
static short counter; /*Smaple counter*/
static short counter_limit; /*Smaple counter limit*/
static short control; /*LFO Control*/
static short max_freq; /*Maximum notch center frequency*/
static short min_freq; /*Minimum notch center frequency*/
static double pha_mix; /*Filtered signal mix*/
static short f_step; /*Sweep frequency step*/
static double dir_mix; /*Direct signal mix*/
static struct br_filter H[PH_STAGES]; /*Array of notch filters stages*/
/*
This funtion initializes the phaser control variables
and the variable notch filter coefficients array
*/
void Phaser_init(short effect_rate,short sampling,short maxf,short minf,short Q,double gainfactor,double pha_mixume,short freq_step, double dmix) {
/*Initialize notch filter coefficients set array*/
br_iir_init(sampling,gainfactor,Q,freq_step, minf);
/*Initializes the phaser control variables*/
center_freq = 0;
samp_freq = sampling;
counter = effect_rate;
control = 0;
counter_limit = effect_rate;
/*Convert frequencies to integer indexes*/
min_freq = 0;
max_freq = (maxf - minf)/freq_step;
pha_mix = pha_mixume;
f_step = freq_step;
dir_mix = dmix;
}
/*
This function does the actual phasing processing
1. It takes the input sample and pass it trough the
cascaded notch filter stages
2. It takes tha output of the cascaded notch filters
and scales it, scales the input sample and generate
the output effect sample.
*/
double Phaser_process(double xin) {
double yout;
int i;
yout = br_iir_filter(xin,&H[0]);
for(i = 1; i < PH_STAGES; i++) {
yout = br_iir_filter(yout,&H[i]);
}
yout = dir_mix*xin + pha_mix*yout;
return yout;
}
/*
This function makes vary the center notch frequency
in all the cascaded notch filter stages by a simulated
triangle wave LFO that goes up and down
*/
void Phaser_sweep(void) {
int i;
if (!--counter) {
if (!control) {
center_freq+=f_step;
if (center_freq > max_freq) {
control = 1;
}
}
else if (control) {
center_freq-=f_step;
if (center_freq == min_freq) {
control = 0;
}
}
for(i = 0; i < PH_STAGES; i++) {
br_iir_setup(&H[i],center_freq);
}
counter = counter_limit;
}
}
/************
Phaser.h
***********/
#ifndef __PHASER_H__
#define __PHASER_H__
extern void Phaser_init(short effect_rate,short sampling,short maxf,short minf,short Q,double gainfactor,double mix_volume,short freq_step, double dmix);
extern double Phaser_process(double xin);
extern void Phaser_sweep(void);
#endifVariable 2nd order band pass IIR filter
#include "bp_iir.h"
#include <math.h>
#define BP_MAX_COEFS 120
#define PI 3.1415926
/*This is an array of the filter parameters object
defined in the br_iir.h file*/
static struct bp_coeffs bp_coeff_arr[BP_MAX_COEFS];
/*This initialization function will create the
band pass filter coefficients array, you have to specify:
fsfilt = Sampling Frequency
gb = Gain at cut frequencies
Q = Q factor, Higher Q gives narrower band
fstep = Frequency step to increase center frequencies
in the array
fmin = Minimum frequency for the range of center frequencies
*/
void bp_iir_init(double fsfilt,double gb,double Q,short fstep, short fmin) {
int i;
double damp;
double wo;
damp = gb/sqrt(1 - pow(gb,2));
for (i=0;i<BP_MAX_COEFS;i++) {
wo = 2*PI*(fstep*i + fmin)/fsfilt;
bp_coeff_arr[i].e = 1/(1 + damp*tan(wo/(Q*2)));
bp_coeff_arr[i].p = cos(wo);
bp_coeff_arr[i].d[0] = (1-bp_coeff_arr[i].e);
bp_coeff_arr[i].d[1] = 2*bp_coeff_arr[i].e*bp_coeff_arr[i].p;
bp_coeff_arr[i].d[2] = (2*bp_coeff_arr[i].e-1);
}
}
/*This function loads a given set of band pass filter coefficients acording to a center frequency index
into a band pass filter object
H = filter object
ind = index of the array mapped to a center frequency
*/
void bp_iir_setup(struct bp_filter * H,int ind) {
H->e = bp_coeff_arr[ind].e;
H->p = bp_coeff_arr[ind].p;
H->d[0] = bp_coeff_arr[ind].d[0];
H->d[1] = bp_coeff_arr[ind].d[1];
H->d[2] = bp_coeff_arr[ind].d[2];
}
/*This function loads a given set of band pass filter coefficients acording to a center frequency index
into a band pass filter object
H = filter object
ind = index of the array mapped to a center frequency
*/
double bp_iir_filter(double yin,struct bp_filter * H) {
double yout;
H->x[0] = H->x[1];
H->x[1] = H->x[2];
H->x[2] = yin;
H->y[0] = H->y[1];
H->y[1] = H->y[2];
H->y[2] = H->d[0]* H->x[2] - H->d[0]* H->x[0] + (H->d[1]* H->y[1]) - H->d[2]* H->y[0];
yout = H->y[2];
return yout;
}
/******************
bp_iir.h
**********************/
#ifndef __BP_IIR_H__
#define __BP_IIR_H__
struct bp_coeffs{
double e;
double p;
double d[3];
};
struct bp_filter{
double e;
double p;
double d[3];
double x[3];
double y[3];
};
extern void bp_iir_init(double fsfilt,double gb,double Q,short fstep, short fmin);
extern void bp_iir_setup(struct bp_filter * H,int index);
extern double bp_iir_filter(double yin,struct bp_filter * H);
#endifDigital Comb Filter
% Usage: [B,A] = COMB(order, scalar);
%
% ORDER is the number of samples delayed prior to add
% SCALAR is the coefficient that will be applied to
% the delayed signal path at the final summation block.
%
% Note, there are two types of comb filters. A DC-blocker and a DC-passer.
% To get a DC-Blocker (tooth at DC), pass in a -1 for the scalar.
% To get a DC-Passer (+6dB at DC), pass in a +1 for the scalar.
%
% By default, if the scalar is not passed, a DC-Passer is assumed.
%
% Author: sparafucile17 03/16/04
% Validate that the proper argument count was supplied
error(nargchk(1, 2, nargin));
% Use scalar is passed as an argument, otherwise assume scalar=1;
if (nargin == 1)
scalar = 1;
else
scalar = varargin{1};
end
% Input has zeros to simulate a single delay of N samples
a = [ 1 ];
b = [ 1 zeros(1, N-1) scalar*1];Equal Loudness Curves (ISO226)
function [spl, freq] = iso226(phon);
%
% Generates an Equal Loudness Contour as described in ISO 226
%
% Usage: [SPL FREQ] = ISO226(PHON);
%
% PHON is the phon value in dB SPL that you want the equal
% loudness curve to represent. (1phon = 1dB @ 1kHz)
% SPL is the Sound Pressure Level amplitude returned for
% each of the 29 frequencies evaluated by ISO226.
% FREQ is the returned vector of frequencies that ISO226
% evaluates to generate the contour.
%
% Desc: This function will return the equal loudness contour for
% your desired phon level. The frequencies evaulated in this
% function only span from 20Hz - 12.5kHz, and only 29 selective
% frequencies are covered. This is the limitation of the ISO
% standard.
%
% In addition the valid phon range should be 0 - 90 dB SPL.
% Values outside this range do not have experimental values
% and their contours should be treated as inaccurate.
%
% If more samples are required you should be able to easily
% interpolate these values using spline().
%
% Author: sparafucile17 03/01/05
% /---------------------------------------\
%%%%%%%%%%%%%%%%% TABLES FROM ISO226 %%%%%%%%%%%%%%%%%
% \---------------------------------------/
f = [20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 800 ...
1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000 12500];
af = [0.532 0.506 0.480 0.455 0.432 0.409 0.387 0.367 0.349 0.330 0.315 ...
0.301 0.288 0.276 0.267 0.259 0.253 0.250 0.246 0.244 0.243 0.243 ...
0.243 0.242 0.242 0.245 0.254 0.271 0.301];
Lu = [-31.6 -27.2 -23.0 -19.1 -15.9 -13.0 -10.3 -8.1 -6.2 -4.5 -3.1 ...
-2.0 -1.1 -0.4 0.0 0.3 0.5 0.0 -2.7 -4.1 -1.0 1.7 ...
2.5 1.2 -2.1 -7.1 -11.2 -10.7 -3.1];
Tf = [ 78.5 68.7 59.5 51.1 44.0 37.5 31.5 26.5 22.1 17.9 14.4 ...
11.4 8.6 6.2 4.4 3.0 2.2 2.4 3.5 1.7 -1.3 -4.2 ...
-6.0 -5.4 -1.5 6.0 12.6 13.9 12.3];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Error Trapping
if((phon < 0) | (phon > 90))
disp('Phon value out of bounds!')
spl = 0;
freq = 0;
else
%Setup user-defined values for equation
Ln = phon;
%Deriving sound pressure level from loudness level (iso226 sect 4.1)
Af=4.47E-3 * (10.^(0.025*Ln) - 1.15) + (0.4*10.^(((Tf+Lu)/10)-9 )).^af;
Lp=((10./af).*log10(Af)) - Lu + 94;
%Return user data
spl = Lp;
freq = f;
endVariable 2nd order Notch IIR filter
#include "br_iir.h"
#include <math.h>
#define BR_MAX_COEFS 120
#define PI 3.1415926
/*This is an array of the filter parameters object
defined in the br_iir.h file*/
static struct br_coeffs br_coeff_arr[BR_MAX_COEFS];
/*This initialization function will create the
notch filter coefficients array, you have to specify:
fsfilt = Sampling Frequency
gb = Gain at cut frequencies
Q = Q factor, Higher Q gives narrower band
fstep = Frequency step to increase center frequencies
in the array
fmin = Minimum frequency for the range of center frequencies
*/
void br_iir_init(double fsfilt,double gb,double Q,double fstep, short fmin) {
int i;
double damp;
double wo;
damp = sqrt(1 - pow(gb,2))/gb;
for (i=0;i<BR_MAX_COEFS;i++) {
wo = 2*PI*(fstep*i + fmin)/fsfilt;
br_coeff_arr[i].e = 1/(1 + damp*tan(wo/(Q*2)));
br_coeff_arr[i].p = cos(wo);
br_coeff_arr[i].d[0] = br_coeff_arr[i].e;
br_coeff_arr[i].d[1] = 2*br_coeff_arr[i].e*br_coeff_arr[i].p;
br_coeff_arr[i].d[2] = (2*br_coeff_arr[i].e-1);
}
}
/*This function loads a given set of notch filter coefficients acording to a center frequency index
into a notch filter object
H = filter object
ind = index of the array mapped to a center frequency
*/
void br_iir_setup(struct br_filter * H,int ind) {
H->e = br_coeff_arr[ind].e;
H->p = br_coeff_arr[ind].p;
H->d[0] = br_coeff_arr[ind].d[0];
H->d[1] = br_coeff_arr[ind].d[1];
H->d[2] = br_coeff_arr[ind].d[2];
}
/*
This function does the actual notch filter processing
yin = input sample
H = filter object
*/
double br_iir_filter(double yin,struct br_filter * H) {
double yout;
H->x[0] = H->x[1];
H->x[1] = H->x[2];
H->x[2] = yin;
H->y[0] = H->y[1];
H->y[1] = H->y[2];
H->y[2] = H->d[0]*H->x[2] - H->d[1]*H->x[1] + H->d[0]*H->x[0] + H->d[1]*H->y[1] - H->d[2]*H->y[0];
yout = H->y[2];
return yout;
}
/******************/
#ifndef __BR_IIR_H__
#define __BR_IIR_H__
/*
Notch filter coefficients object
*/
struct br_coeffs {
double e;
double p;
double d[3];
};
/*
Notch filter object
*/
struct br_filter {
double e;
double p;
double d[3];
double x[3];
double y[3];
};
extern void br_iir_init(double fsfilt,double gb,double Q,double fstep, short fmin);
extern void br_iir_setup(struct br_filter * H,int index);
extern double br_iir_filter(double yin,struct br_filter * H);
#endifFarrow 2-D Image Interpolation
function Demo
%Author: J. Selva.
%Date: June 2011.
close all
format short g
format compact
T = 1; %Sampling period in both dimensions.
B = 1/(1.223*T); %This is the two-sided bandwidth in both dimensions. It is a typical ...
%value in SAR image co-registration.
P = 18; %Kernel semi-length.
Q = 10; %Polynomial order in Farrow structure.
%The next call to "ChebApproxPolys" computes all Farrow coefficients.
t0 = -P*T-T/2; %The parameters t0, Nt, 2*P+1 specify the 2*P+1 intervals in which Knab's ...
%pulse is to be approximated by polynomials.
Nt = 2*P+1;
Dt = T;
IOut = [-T/2,T/2];
%Perform Chebyshev interpolation. The parameters T, B and P are passed to KnabAPPulse.
K = ChebApproxPolys('KnabAPPulse',t0,Nt,Dt,IOut,40,Q,T,B,P).';
%Evaluate the test image in a regular grid with spacing 1.
Lx = 100;
Ly = 100;
x0 = (0:Lx-1)';
x0 = x0(:,ones(1,Lx));
y0 = x0.';
x0 = x0(:);
y0 = y0(:);
disp('Computing image in a regular grid...')
tic;
Image = reshape(TestImage5(x0,y0,B),[Lx,Ly]);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
clear x0 y0
%Specify NP points on the image at random.
NP = 10000;
x = P*T+rand(1,NP)*(Lx-2*P)*T;
y = P*T+rand(1,NP)*(Ly-2*P)*T;
%Evaluate image directly.
disp('Computing image at selected positions exactly...')
tic;
vRef = TestImage5(x,y,B);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
%Interpolate the image using the 2D-Farrow structure.
disp('Computing image at selected positions using Farrow 2-D...')
tic;
v = Farrow2DInterp(Image,K,x,y);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
%Show the error.
disp(['Maximum interpolation error: ', ...
num2str(20*log10(max(abs(vRef-v))/max(abs(vRef)))),' dB'])
disp(' ')
disp('Display image and interpolation error:')
disp(' ')
%Evaluate the test image in a regular grid with spacing 1.
Lx = 100;
Ly = Lx;
x1 = (P*T:0.5*T:(Lx-P)*T)';
x1 = x1(:,ones(1,length(x1)));
y1 = x1.';
x1 = x1(:);
y1 = y1(:);
tic;
disp('Computing image in over-sampled grid ...')
vRef = TestImage5(x1,y1,B);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
tic;
disp('Computing signal in over-sampled grid using 2-D Farrow ...')
v = Farrow2DInterp(Image,K,x1,y1);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
x1 = (P*T:0.5*T:(Lx-P)*T)';
y1 = x1.';
disp(['Maximum interpolation error: ', ...
num2str(20*log10(max(abs(vRef-v))/max(abs(vRef)))),' dB'])
vRef = reshape(vRef,[length(x1),length(y1)]);
v = reshape(v,[length(x1),length(y1)]);
mesh(x1,y1,abs(vRef))
xlabel('x coordinate')
ylabel('y coordinate')
zlabel('Image''s abs. value')
title('Test Image (abs. value)')
figure
mesh(x1,y1,20*log10(abs(v-vRef)/max(abs(vRef(:)))))
xlabel('x coordinate')
ylabel('y coordinate')
zlabel('Interp. error (dB)')
title('Interpolation error (dB)')
function P = ChebApproxPolys(Func,t0,Nt,Dt,IOut,NPoints,Q,varargin)
%This function constructs Chebyshev interpolation polynomials in a set of consecutive ...
%intervals.
%
% Func is the function to be interpolated.
% t0 is the first instant of the first interval.
% Nt is the number of intervals.
% IOut is the interval in which the output polynomials will be evaluated.
% Dt is the interval width. The intervals are [t0,t0+Dt], [t0+Dt,t0+2*Dt],....
% [t0+(Nt-1)*Dt,t0+Nt*Dt].
% NPoints is the number of points in which Func is evaluated, for each interval.
% Q is the output polynomial order.
% varargin contains extra parameters which are passed to varargin.
%Construct matrix that maps function values to polynomial coefficients.
%The algorithm is that in section 9.8 of
%
% W. H. Press et al., Numerical Recipes in C. Cambridge University
% Press, 1997.
M = MapPolyMatrix([-1,1],IOut,Q) * ChebCoefMatrix(Q) * ...
BlockChebInvMatrix(NPoints,Q);
% Get Chebyshev roots.
r = ChebRoots(NPoints);
r = r(:);
P = zeros(Q,Nt);
%Compute polynomial for each interval.
for k = 0:Nt-1,
v1 = feval(Func,MapIaToIb(r,[-1,1],t0+[k,k+1]*Dt),varargin{:});
P(:,k+1) = double(M * v1(:));
end
function M = MapPolyMatrix(IIn,IOut,Q)
%For a polynomial with Q coefficients which is evaluated in interval IIn, this function ...
%gives the matrix that transforms its coefficients, so that it can be evaluated in the ...
%associated points in interval IOut.
ayx = (IIn(1)-IIn(2))/(IOut(1)-IOut(2));
byx = IIn(1)-ayx*IOut(1);
if byx == 0
M = diag(ayx.^(0:Q-1));
else
M = zeros(Q);
for q = 0:Q-1
for r = 0:q
M(r+1,q+1) = nchoosek(q,r) * ayx^r * byx^(q-r);
end
end
end
function M = ChebCoefMatrix(N)
M = zeros(N);
M(1) = 1;
if N > 1
M(2,2) = 1;
end
for n = 3:N,
M(2:n,n) = 2*M(1:n-1,n-1);
M(1:n-2,n) = M(1:n-2,n) - M(1:n-2,n-2);
end
function M = BlockChebInvMatrix(N,Q)
r = ChebRoots(N);
r= r(:);
M = zeros(N,Q);
if Q >=1
M(:,1) = 1;
end
if Q >= 2
M(:,2) = r;
end
for q = 3:Q,
M(:,q) = 2 * r .* M(:,q-1) - M(:,q-2);
end
M = M.'/N;
M(2:end,:) = M(2:end,:) * 2;
function R = ChebRoots(N)
R = cos(pi*((1:N)-1/2)/N);
function y = MapIaToIb(x,Ia,Ib)
%Linear transformation that maps interval Ia onto interval Ib.
y = (Ib(1)+Ib(2))/2 + (x-(Ia(1)+Ia(2))/2) * ((Ib(2)-Ib(1))/(Ia(2)-Ia(1)));
function [g,BoundAP] = KnabAPPulse(t,T,B,P)
%Knab interpolation pulse in
%
% J. J. Knab, 'Interpolation of band-limited functions using the Approximate
% Prolate series', IEEE Transactions on Information Theory, vol. IT-25,
% no. 6, pp. 717-720, Nov 1979.
%
%Grid points -P*T:T:P*T.
%Interpolation interval: [-T/2,T/2].
g = sinc(t/T) .* sinc((1/T-B)*sqrt(t.^2-(P*T)^2)) / sinc(j*(1-B*T)*P);
if nargout == 2
BoundAP = 1/sinh(P*pi*(1-B*T));
end
function v = TestImage5(x,y,B)
%Test image formed by a sum of random sincs.
Seed = rand('state');
rand('state',362436069);
L = 100; %Image size in both dimensions.
Np = 5000;
p = rand(Np,2)*L;
A = exp(j*2*pi*rand(Np,1));
rand('state',Seed);
v = zeros(size(x));
for k = 1:length(x(:))
v(k) = A(:).' * (sinc(B*(x(k)-p(:,1))) .* sinc(B*(y(k)-p(:,2))));
end
function vx = Farrow2DInterp(Image,Kernel,x,y)
%Farrow interpolation in two dimensions using interlaced FFTs.
vxSize = size(x);
[LK,NK] = size(Kernel);
[Lx,Ly] = size(Image);
nx = round(x);
ny = round(y);
ux = x-nx;
uy = y-ny;
np = 1+nx+ny*Lx;
np = np(:);
clear x y nx ny
P = (LK-1)/2;
NFFTx = 2^ceil(log2(Lx+LK-1));
NFFTy = 2^ceil(log2(Ly+LK-1));
KernelFx = fft(Kernel,NFFTx);
KernelFy = fft(Kernel,NFFTy);
ImageFx = fft(Image,NFFTx);
vx = zeros(size(np))+j*zeros(size(np));
for kx = NK:-1:1,
ImageCx = ifft(ImageFx .* KernelFx(:,kx(ones(1,Ly))),NFFTx);
ImageCx = ImageCx(1+P:Lx+LK-1-P,:);
ImageCxFy = fft(ImageCx,NFFTy,2);
vy = zeros(size(np))+j*zeros(size(np));
for ky = NK:-1:1,
ImageCxCy = ifft(ImageCxFy .* KernelFy(:,ky(ones(1,Lx))).',NFFTy,2);
ImageCxCy = ImageCxCy(:,1+P:Ly+LK-1-P);
vy(:) = vy(:) .* uy(:) + ImageCxCy(np);
end
vx(:) = vx(:) .* ux(:) + vy(:);
end
vx = reshape(vx,vxSize);Shelving Filter Design
function [b, a] = shelving(G, fc, fs, Q, type)
%
% Derive coefficients for a shelving filter with a given amplitude and
% cutoff frequency. All coefficients are calculated as described in
% Zolzer's DAFX book (p. 50 -55).
%
% Usage: [B,A] = shelving(G, Fc, Fs, Q, type);
%
% G is the logrithmic gain (in dB)
% FC is the center frequency
% Fs is the sampling rate
% Q adjusts the slope be replacing the sqrt(2) term
% type is a character string defining filter type
% Choices are: 'Base_Shelf' or 'Treble_Shelf'
%
% Author: sparafucile17 08/22/05
%
%Error Check
if((strcmp(type,'Base_Shelf') ~= 1) && (strcmp(type,'Treble_Shelf') ~= 1))
error(['Unsupported Filter Type: ' type]);
end
K = tan((pi * fc)/fs);
V0 = 10^(G/20);
root2 = 1/Q; %sqrt(2)
%Invert gain if a cut
if(V0 < 1)
V0 = 1/V0;
end
%%%%%%%%%%%%%%%%%%%%
% BASE BOOST
%%%%%%%%%%%%%%%%%%%%
if(( G > 0 ) & (strcmp(type,'Base_Shelf')))
b0 = (1 + sqrt(V0)*root2*K + V0*K^2) / (1 + root2*K + K^2);
b1 = (2 * (V0*K^2 - 1) ) / (1 + root2*K + K^2);
b2 = (1 - sqrt(V0)*root2*K + V0*K^2) / (1 + root2*K + K^2);
a1 = (2 * (K^2 - 1) ) / (1 + root2*K + K^2);
a2 = (1 - root2*K + K^2) / (1 + root2*K + K^2);
%%%%%%%%%%%%%%%%%%%%
% BASE CUT
%%%%%%%%%%%%%%%%%%%%
elseif (( G < 0 ) & (strcmp(type,'Base_Shelf')))
b0 = (1 + root2*K + K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);
b1 = (2 * (K^2 - 1) ) / (1 + root2*sqrt(V0)*K + V0*K^2);
b2 = (1 - root2*K + K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);
a1 = (2 * (V0*K^2 - 1) ) / (1 + root2*sqrt(V0)*K + V0*K^2);
a2 = (1 - root2*sqrt(V0)*K + V0*K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);
%%%%%%%%%%%%%%%%%%%%
% TREBLE BOOST
%%%%%%%%%%%%%%%%%%%%
elseif (( G > 0 ) & (strcmp(type,'Treble_Shelf')))
b0 = (V0 + root2*sqrt(V0)*K + K^2) / (1 + root2*K + K^2);
b1 = (2 * (K^2 - V0) ) / (1 + root2*K + K^2);
b2 = (V0 - root2*sqrt(V0)*K + K^2) / (1 + root2*K + K^2);
a1 = (2 * (K^2 - 1) ) / (1 + root2*K + K^2);
a2 = (1 - root2*K + K^2) / (1 + root2*K + K^2);
%%%%%%%%%%%%%%%%%%%%
% TREBLE CUT
%%%%%%%%%%%%%%%%%%%%
elseif (( G < 0 ) & (strcmp(type,'Treble_Shelf')))
b0 = (1 + root2*K + K^2) / (V0 + root2*sqrt(V0)*K + K^2);
b1 = (2 * (K^2 - 1) ) / (V0 + root2*sqrt(V0)*K + K^2);
b2 = (1 - root2*K + K^2) / (V0 + root2*sqrt(V0)*K + K^2);
a1 = (2 * ((K^2)/V0 - 1) ) / (1 + root2/sqrt(V0)*K + (K^2)/V0);
a2 = (1 - root2/sqrt(V0)*K + (K^2)/V0) / (1 + root2/sqrt(V0)*K + (K^2)/V0);
%%%%%%%%%%%%%%%%%%%%
% All-Pass
%%%%%%%%%%%%%%%%%%%%
else
b0 = V0;
b1 = 0;
b2 = 0;
a1 = 0;
a2 = 0;
end
%return values
a = [ 1, a1, a2];
b = [ b0, b1, b2];
