% ----------------------------------------------------
% Title: Discrete Wavelet Transform
% by means of a linear convolution matrix
%
% Author: David Valencia
% UPIITA IPN
%
% Posted at: http://www.dsprelated.com/showcode/11.php
%
% Description:
% Computes the Discrete Wavelet Transform
% of n levels (and its branches) from a n sample input.
% Base filters can be of any lenght.
%
% Generalized DWT filter bank computed via the
% convolution matrix method (linear convolution)
%
% Dependencies:
%
% formfilter.m
% http://www.dsprelated.com/showcode/12.php
% upsample2.m
% http://www.dsprelated.com/showcode/10.php
%
% Revisions:
% v1.0a Commented and translated in English
%
% ----------------------------------------------------
close all; clear all; clc;
disp('== GENERALIZED MATRICIAL DWT ==')
%% STEP 1 - Base Filters
% Select type of base filters
typeofbasis = 'o';
typbior = 'bior2.2';
typor = 'db3';
% Obtain base filters
if(typeofbasis == 'b')
[Rf,Df] = biorwavf(typbior);
[h0,h1,g0,g1] = biorfilt(Df,Rf);
elseif (typeofbasis == 'o')
[h0,h1,g0,g1] = wfilters(typor);
end;
%% STEP 2 - Parameter configuration
% One can declare or recover an input vector from another
% program, but here an example vector is built for testing.
N = 16; %Number of input samples
x = (1:N)';
L = length(h0); %Base filter lenght
n_stages = 2; %Of the low-pass stage
n_branches = 2^n_stages; %Number of branches
dec_fact = n_branches; %Decimating factor
% L = Basic filters lenght (h0 ó h1)
% N = Input vector length
% n_stages = number of stages, it generates (2^n_stages) branches
hx = formfilter(n_stages, 1, h0, h1);
Lhx = length(hx);
fprintf('\nLinear convolution matrix processing\n');
%% STEP 3: Build analysis matrices
% -------------------------------------------------------------
% High-pass analysis matrix (single scale) [linear convolution]
% As the DWT has a simple high-pass branch, the matrix is easy
% to build with the base filter.
% Append zeros to each side of the vector to use it as a sliding
% part to build the matrix
movil_1 = [zeros(1,N-1),fliplr(h1),zeros(1,N-1)];
lm = length(movil_1);
n_rows = ceil((L+N-1)/2); %Compute the number of rows
W1 = zeros(n_rows,N); % Allocate memory for the high pass matrix
dec_factor = 2; %Decimate factor for the high-pass matrix
% Build matrix (with sequential shifting)
for i = 0:(n_rows-1)
W1(i+1,:) = movil_1(lm-(i*dec_factor)-N+1 : lm-(i*dec_factor)); %Rama pasaaltas
end
disp('High-pass analysis matrix: ');
W1
disp('High-pass branch results');
y1 = W1 * x
% -------------------------------------------------------------
% Low-pass several stage analysis matrix [linear convolution]
% Get an equivalent filter to get its length
hx = formfilter(n_stages, 1, h0, h1);
Lhx = length(hx);
% Compute the number of rows per branch
rowsperfilter = ceil((N+Lhx)/dec_fact);
% Allocate memory for the low-pass matrix
W0 = zeros((n_branches/2)*rowsperfilter,N);
tic; %Start cronometer
% Build low pass filter
for i = 0:(n_branches/2)-1
% Get the equivalent filter depending on the number of stages
if n_stages < 3
hx = formfilter(n_stages, (n_branches/2) - i, h0, h1);
else
hx = formfilter(n_stages, (n_branches/2) - i - 1, h0, h1);
end
% Append zeros to the vector
movil = [zeros(1,N-1),fliplr(hx),zeros(1,N-1)];
lm = length(movil);
% Shift and add the vector to the matrix
for j = 0:rowsperfilter-1
W0(i*rowsperfilter + j + 1,:) = movil(lm-(dec_fact*j)-N + 1 : lm-(dec_fact*j));
end
end
disp('Low-pass many stage analysis matrix: ');
W0
disp('Low-pass filter bank results: ');
y0 = W0*x
%% STEP 4: Recover signal
disp('Recovered signal: ');
xn1 = (W1')*y1;
xn0 = (W0')*y0;
xn = xn1+xn0
%% STEP 5: Compute error
err = (x - xn).^2;
disp('Error :');
err = (sum(err)/N);
disp(err);
fprintf('\nProcessing time = %d seg \n',toc);
fprintf('\nProgram end >>> \n');
//Design of FIR Filter using Frquency Sampling Technique
//Low Pass Filter Design
//Cutoff Frequency Wc = pi/2
//M = Filter Lenth = 7
clear;
clc;
M = 7;
N = ((M-1)/2)+1
wc = %pi/2;
for k =1:M
w(k) = ((2*%pi)/M)*(k-1);
if (w(k)>=wc)
k-1
break
end
end
for i = 1:k-1
Hr(i) = 1;
G(i) = ((-1)^(i-1))*Hr(i);
end
for i = k:N
Hr(i) = 0;
G(i) = ((-1)^(i-1))*Hr(i);
end
h = zeros(1,M);
for n = 1:M
for k = 2:N
h(n) = G(k)*cos((2*%pi/M)*(k-1)*((n-1)+(1/2)))+h(n);
end
h(n) = (1/M)* (G(1)+2*h(n));
end
[hzm,fr]=frmag(h,256);
hzm_dB = 20*log10(hzm)./max(hzm);
plot(2*fr,hzm_dB)
xtitle('Frequency Response of LPF with Normalized cutoff =0.5','Normalized Frequency W ------>','Magnitude Response H(w)---->');
xgrid(1);
// Hamming Weight and Hamming Distance
//H(7,4)
//Code Word Length = 7, Message Word length = 4, Parity bits =3
close;
clc;
//Getting Code Words
code1 = input('Enter the first code word');
code2 = input('Enter the second code word');
Hamming_Distance = 0;
for i = 1:length(code1)
Hamming_Distance =Hamming_Distance+xor(code1(i),code2(i));
end
disp(Hamming_Distance,'Hamming Distance')
//Result
//Enter the first code word [0,1,1,1,0,0,1]
//Enter the second code word[1,1,0,0,1,0,1]
//Hamming Distance 4.
%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
%Sampling Rate
Fs = 48000;
%Analog C-weighting filter according to IEC/CD 1672.
1672.
f1 = 20.598997;
f4 = 12194.217;
C1000 = 0.0619;
pi = 3.14159265358979;
NUM = [ (2*pi*f4)^2*(10^(C1000/20)) 0 0 ];
DEN = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]);
%Bilinear transformation of analog design to get the digital [b,a] = bilinear(NUM,DEN,Fs);
function y = peak2peak(signal)
% Return the peak-to-peak amplitude of the supplied signal. This is the
% same as max(signal) minus min(signal).
%
% Usage: y = PEAK2PEAK(signal);
%
% SIGNAL is your one-dimensional input array
%
% Author: sparafucile17
% Input must have some length
if(length(signal) == 1)
error('ERROR: input signal must have more than one element');
end
% This function only supports one-dimensional arrays
if((size(signal, 2) ~= 1) && (size(signal, 1) ~= 1))
error('ERROR: Input must be one-dimensional');
end
% Find the peak and return it
min_sig = min(signal);
max_sig = max(signal);
% Answer
y = max_sig - min_sig;
function [per_error] = percent_error(measured, actual)
% Creates the percent error between measured value and actual value
%
% Usage: percent_error(MEASURED, ACTUAL);
%
% Measured is the your result
% Actual is the value that your result should be
%
% Author: sparafucile17
per_error = abs(( (measured - actual) ./ actual ) * 100);
function [time,amp] = minima(signal, len)
%
% Finds the local minimum values of the given signal.
%
% Usage: [IND,AMP] = MINIMA(X, N);
%
% N is the number of points that need to be evaluated
% (Normally equal to LENGTH(X))
% X is the data
% IND contains the indices of X that contain the minima data
% AMP contains the minima amplitudes
%
% Author: sparafucile17 10/02/2003
%Initialize data
index = 1;
prev = signal(1);
local_time = 0;
local_amp = 0;
prev_slope = 1; %allow a maxima point at the second sample
%prevent length error
if(len > length(signal))
len = length(signal)
end
%Loop through data and find local minimums
for loop_i=2:len
cur = signal(loop_i);
slope = cur - prev;
if((cur < prev) & (prev_slope > 0)) %Positive slope and amplitude dips
local_time(index) = loop_i-1;
local_amp(index) = prev;
index = index+1;
end
prev = cur;
prev_slope = slope;
end
%After loop assign data to output variables
time = local_time;
amp = local_amp;
function [time,amp] = maxima(signal, thresh, len)
% Finds the local maximum values of the given signal.
%
% Usage: [IND,AMP] = MAXIMA(X, THRESH, N);
%
% X is the data
% THRESH is the threshold that signal must be greater
% than in order to be counted. (Default = 0.0)
% N is the number of points that need to be evaluated
% (Defaults = LENGTH(X))
% IND contains the indices of X that contain the maxima data
% AMP contains the maxima amplitudes
%
% Author: sparafucile17 10/02/2003
if(exist('len') == 0)
len = length(signal);
end
if(exist('thresh') == 0)
thresh = 0.0;
end
%Initialize data
index = 1;
prev = signal(1);
local_time = 0;
local_amp = 0;
prev_slope = 1; %allow a maxima point at the second sample
%prevent length error
if(len > length(signal))
len = length(signal)
end
%Loop through data and find local maximums
for loop_i=2:len
cur = signal(loop_i);
slope = cur - prev;
if((cur < prev) & (prev_slope > 0) & (cur > thresh)) %Positive slope and amplitude dips
local_time(index) = loop_i-1;
local_amp(index) = prev;
index = index+1;
end
prev = cur;
prev_slope = slope;
end
%After loop assign data to output variables
time = local_time;
amp = local_amp;
function out = create_pink_noise(Fs, Sec, Amp)
% Creates a pink noise signal and saves it as a wav file
%
% Usage: create_noise(Fs, Sec, Amp);
%
% Fs is the desired sampling rate
% Sec is the duration of the signal in seconds
% Amp is the amplitude in dB of the signal (0dB to -144dB)
%
% Author: sparafucile17 06/14/02
%error trapping
if((Amp > 0) || (Amp < -144))
error('Amplitude is not within the range of 0dB to -144dB');
end
%Create Whitenoise
white_noise = randn((Fs*Sec)+1,1);
%Apply weighted sum of first order filters to approximate a -10dB/decade
%filter. This is Paul Kellet's "refined" method (a.k.a instrumentation
%grade) It is accurate to within +/-0.05dB above 9.2Hz
b=zeros(7,1);
for i=1:((Fs*Sec)+1)
b(1) = 0.99886 * b(1) + white_noise(i) * 0.0555179;
b(2) = 0.99332 * b(2) + white_noise(i) * 0.0750759;
b(3) = 0.96900 * b(3) + white_noise(i) * 0.1538520;
b(4) = 0.86650 * b(4) + white_noise(i) * 0.3104856;
b(5) = 0.55000 * b(5) + white_noise(i) * 0.5329522;
b(6) = -0.7616 * b(6) - white_noise(i) * 0.0168980;
pink_noise(i) = b(1) + b(2) + b(3) + b(4) + b(5) + b(6) + b(7) + white_noise(i) * 0.5362;
b(7) = white_noise(i) * 0.115926;
end
%Normalize to +/- 1
if(abs(min(pink_noise)) > max(pink_noise))
pink_noise = pink_noise / abs(min(pink_noise));
else
pink_noise = pink_noise / max(pink_noise);
end
%Normalize to prevent positive saturation (We can't represent +1.0)
pink_noise = pink_noise /abs(((2^31)-1)/(2^31));
%Scale signal to match desired level
pink_noise = pink_noise * 10^(Amp/20);
%Output noise signal
out = pink_noise(1:end-1);
%Sampling Rate
Fs = 48000;
%Analog C-weighting filter according to IEC/CD 1672.
1672.
f1 = 20.598997;
f4 = 12194.217;
C1000 = 0.0619;
pi = 3.14159265358979;
NUM = [ (2*pi*f4)^2*(10^(C1000/20)) 0 0 ];
DEN = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]);
%Bilinear transformation of analog design to get the digital [b,a] = bilinear(NUM,DEN,Fs);
%Sampling Rate
Fs = 48000;
%Analog A-weighting filter according to IEC/CD 1672.
f1 = 20.598997;
f2 = 107.65265;
f3 = 737.86223;
f4 = 12194.217;
A1000 = 1.9997;
pi = 3.14159265358979;
NUM = [ (2*pi*f4)^2*(10^(A1000/20)) 0 0 0 0 ];
DEN = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]);
DEN = conv(conv(DEN,[1 2*pi*f3]),[1 2*pi*f2]);
%Bilinear transformation of analog design to get the digital filter.
[b,a] = bilinear(NUM,DEN,Fs);
/*************** AutoWah.c ********************************/
#include "bp_iir.h"
#include "AutoWah.h"
static short center_freq;
static short samp_freq;
static short counter;
static short counter_limit;
static short control;
static short max_freq;
static short min_freq;
static short f_step;
static struct bp_filter H;
/*
This is the auto wah effect initialization function.
This initializes the band pass filter and the effect control variables
*/
void AutoWah_init(short effect_rate,short sampling,short maxf,short minf,short Q,double gainfactor,short freq_step) {
double C;
/*Process variables*/
center_freq = 0;
samp_freq = sampling;
counter = effect_rate;
control = 0;
/*User Parametters*/
counter_limit = effect_rate;
/*Convert frequencies to index ranges*/
min_freq = 0;
max_freq = (maxf - minf)/freq_step;
bp_iir_init(sampling,gainfactor,Q,freq_step,minf);
f_step = freq_step;
}
/*
This function generates the current output value
Note that if the input and output signal are integer
unsigned types, we need to add a half scale offset
*/
double AutoWah_process(int xin) {
double yout;
yout = bp_iir_filter(xin,&H);
#ifdef INPUT_UNSIGNED
yout += 32767;
#endif
return yout;
}
/*
This function will emulate a LFO that will vary according
to the effect_rate parameter set in the AutoWah_init function.
*/
void AutoWah_sweep(void) {
double yout;
if (!--counter) {
if (!control) {
bp_iir_setup(&H,(center_freq+=f_step));
if (center_freq > max_freq) {
control = 1;
}
}
else if (control) {
bp_iir_setup(&H,(center_freq-=f_step));
if (center_freq == min_freq) {
control = 0;
}
}
counter = counter_limit;
}
}
/*************** AutoWah.h ****************************/
#ifndef __AUTOWAH_H__
#define __AUTOWAH_H__
extern void AutoWah_init(short effect_rate,short sampling,short maxf,short minf,short Q,double gainfactor,short freq_step);
extern double AutoWah_process(int xin);
extern void AutoWah_sweep(void);
#endif
/************** Usage Example **************************/
#include "AutoWah.h"
void main(void) {
short xin;
short yout;
AutoWah_init(2000, /*Effect rate 2000*/
16000, /*Sampling Frequency*/
1000, /*Maximum frequency*/
500, /*Minimum frequency*/
4, /*Q*/
0.707, /*Gain factor*/
10 /*Frequency increment*/
);
while(1) {
if (new_sample_flag()) {
/*When there's new sample at your ADC or CODEC input*/
/*Read the sample*/
xin = read_sample();
/*Apply AutoWah_process function to the sample*/
yout =AutoWah_process(xin);
/*Send the output value to your DAC or codec output*/
write_output(yout);
/*Makes the LFO vary*/
AutoWah_sweep();
}
}
}
% ************************************************************
% Spectrum model for GSM signal
% Markus Nentwig, 2011
% based on 3GPP TS 45.004 section 2 "Modulation format for GMSK",
% assuming no additional filtering and an infinite-length
% symbol stream (no slot structure)
% ************************************************************
% ************************************************************
% Parameters
% The "baseline" serves as a gentle reminder that the model
% is only valid near the center frequency.
% For large frequency offsets, one would need more information on
% the particular transmitter used.
% If not available, assume that spectral emission mask requirements
% are met.
% ************************************************************
BW = 285e3;
BW2 = 186e3;
baseline_dB = -76;
% baseline_dB = -999 % disable the constant term
% ************************************************************
% empirical GSM (GMSK narrow-bandwidth pulse) model equation
% ************************************************************
f = (-500e3:1e3:500e3)+1e-3;
gaussPart = exp(-(2*f/BW) .^2);
sincPart = sin(pi*f/BW2) ./ (pi*f/BW2);
flatPart = 10^(baseline_dB/20);
H_dB = 20*log10(abs(gaussPart .* sincPart) + flatPart);
% ************************************************************
% plot the spectrum
% ************************************************************
figure(); grid on; hold on;
h = plot(f/1e6, H_dB, 'b'); set(h, 'linewidth', 2);
% ************************************************************
% plot 'a' GSM spectral emission mask
% note, this is only "an" example
% See 3GPP TS 45.005 section 4.2.1
% "Spectrum due to the modulation and wide band noise"
% section 4.2.1.1
% "General requirements for all types of Base stations and MS"
% note the warning regarding measuring the nominal signal level!
% ************************************************************
maskX_MHz = [0, 100e3, 200e3, 250e3, 400e3, 600e3]/1e6;
maskY_dB = [0.5, 0.5, -30, -33, -60, -60];
h = plot(maskX_MHz, maskY_dB, 'k'); set(h, 'linewidth', 3);
legend('|H(f)|', 'GSM mask example');
xlabel('f/MHz');
ylabel('PSD/dB');
title('GSM spectrum');
# Instantiate the SIIR object. Pass the cutoff frequency
# Fc and the sample rate Fs in Hz. Also define the input
# and output fixed-point type. W=(wl, iwl) where
# wl = word-length and iwl = integer word-length. This
# example uses 23 fraction bits and 1 sign bit.
>>> from siir import SIIR
>>> flt = SIIR(Fstop=1333, Fs=48000, W=(24,0))
# Plot the response of the fixed-point coefficients
>>> plot(flt.hz, 20*log10(flt.h)
# Create a testbench and run a simulation
# (get the simulated response)
>>> from myhdl import Simulation
>>> tb = flt.TestFreqResponse(Nloops=128, Nfft=1024)
>>> Simulation(tb).run()
>>> flt.PlotResponse()
# Happy with results generate the Verilog and VHDL
>>> flt.Convert()
