function [y] = signal_mute(x, index, duration, Fs)
% Return the input vector with a mute that occurs at a specific
% index and with an exponential ramp-down to reduce "pop" sounds.
% The output will start off at 0dB gain and end at -100dB gain.
%
% Usage: y = SIGNAL_MUTE(x, index, duration, Fs);
%
% X is your one-dimensional input array
% INDEX is where in the input signal you want the mute to begin
% DURATION is how long (in seconds) to exponentially ramp down the
% input signal. 100ms is recommended.
% FS is the sample rate
%
% Example:
% You want to mute your signal at 0.5sec with a duration of 100ms
% and with a 48k Fs. (mute complete at 0.6sec)
% y = signal_mute(x, 24000, 0.1, 48000);
%
% Author: sparafucile17 7/29/2003
% Input must have some length
if(length(x) == 1)
error('ERROR: input signal must have more than one element');
end
% This function only supports one-dimensional arrays
if((size(x, 2) ~= 1) && (size(x, 1) ~= 1))
error('ERROR: Input must be one-dimensional');
end
% Make sure there are enough samples to complete the mute
if(length(x) < (index + duration*Fs))
error(['There are not enough samples in X to complete the mute. ' ...
'Either change the mute duration or move the index back.' ]);
end
% Flip vector (temporarily)
if((size(x, 2) ~= 1))
x = x';
flip_back = true;
else
flip_back = false;
end
% Calculate exponential coefficient
dB_atten = -100; %How much attenuation will be at time: index + duration
decayrate = -dB_atten / duration;
coeff = 1.0 - 10^(-decayrate/(20.0*Fs));
% Build the Gain array
gain = [ones(index, 1); zeros((length(x)-index), 1);];
b = [ coeff 0];
a = [ 1 (-1*(1-coeff))];
gain = filter(b,a,gain, [1]);
% Remove minor overshoot at the beginning of gain vector
gain(find(gain > 1)) = 1;
% Apply Mute (gain) to the input signal
y = gain .* x;
% Flip the vector (if required)
if(flip_back == true);
y = y';
end
clc
clear all
close all
fs=20e7;
T=1/fs;
t=0:T:0.0000002;
f=10e6;
scanagle_deg=60;
step_deg=6;
theta=-deg2rad(scanagle_deg):deg2rad(step_deg):deg2rad(scanagle_deg);
for i=1:length(theta);
ant1=2*sin(2*pi*f*t);
ant2=2*sin(2*pi*f*t+theta(i));
[my,mx]=max(ant1);
sum_ant(i)=ant1(mx)+ant2(mx);
diff_ant(i)=ant1(mx)-ant2(mx);
ratio_ant(i)=diff_ant(i)/sum_ant(i);
% if diff_ant(i)==0
% diff_ant(i)=0.1;
% end
% ratio1_ant(i)=sum_ant(i)/diff_ant(i);
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)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sim_ant1=2*sin(2*pi*f*t);
sim_ant2=2*sin(2*pi*f*t+deg2rad(30));
[smy,smx]=max(sim_ant1);
%%%%also take for same sample value of data
sim_sum_ant=sim_ant1(mx)+sim_ant2(mx);
sim_diff_ant=sim_ant1(mx)-sim_ant2(mx);
sim_ratio_ant=sim_diff_ant/sim_sum_ant
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
% ------------------------------------------------------
% Title: formfilter.m
%
% Author: David Valencia
%
% Institution: UPIITA-IPN 2010
%
% for DSPrelated.com
%
% Published in: http://www.dsprelated.com/showcode/12.php
%
% Description: This program receives 2 basic ortogonal
% filters (one is high-pass and one is low-pass)
% the filtering level and the branch number.
% As a result it returns the equivalent filter
% of the specified branch.
%
% Dependencies: upsample2.m
% http://www.dsprelated.com/showcode/10.php
%
% Revision: v1.0a
% - Commented and translated in English
%
% For more details on this code and its implementation
% see the following blog posts:
%
% http://www.dsprelated.com/showarticle/115.php
% http://www.dsprelated.com/showarticle/116.php
%
% ------------------------------------------------------
function [hx] = formfilter(n_stages,branch,h0,h1)
p = branch;
% Seed vector
hx = 0;
hx(1) = 1;
switch n_stages
case 1
% If it is a single stage filter
% the branches are just the base filters
if mod(branch,2) ~= 0
hx = h0;
else
hx = h1;
end
case 2
% For a 2 stage filter bank, the
% equivalent filters are simply
% the convolution between the corresponding
% base filters, and one of them is upsampled
% The use of upsample2 is needed to avoid having
% a large quantitiy of zeros at the end of the vector
% which certainly difficults its usage to build a
% convolution matrix.
switch branch
case 1
hx = conv(h0,upsample2(h0,2));
case 2
hx = conv(h0,upsample2(h1,2));
case 3
hx = conv(h1,upsample2(h0,2));
case 4
hx = conv(h1,upsample2(h1,2));
otherwise
beep;
fprintf('\nFor a 2 stage filter bank there can not be a fifth branch');
end
otherwise
% For a n>2 stages filter bank, a more ellaborated
% process must be made. A series of upsamplings and convolutions
% are needed to get an equivalent vector.
for i=0:n_stages-2
q = floor(p /(2^(n_stages-1-i)));
if (q == 1)
hx = conv(hx,upsample2(h1,2^i));
else
hx = conv(hx,upsample2(h0,2^i));
end
p = mod(p,2^(n_stages-1-i));
end
% Depending on the parity of the branch number, the filter
% goes through a last convolution
t = mod(branch,2);
if(t == 1)
hx = conv(hx,upsample2(h0,2^(n_stages-1)));
else
hx = conv(hx,upsample2(h1,2^(n_stages-1)));
end
end
function [m, b] = ls_linear(x_observed, y_observed)
%
% Perform Least Squares curve-fitting techniques to solve the coefficients
% of the linear equation: y = m*x + b. Input and Output equation
% observations must be fed into the algorithm and a best fit equation will
% be calculated.
%
% Usage: [M,B] = ls_linear(observed_input, observed_output);
%
% observed_input is the x-vector observed data
% observed_output is the y-vector observed data
%
% Author: sparafucile17 06/04/03
%
if(length(x_observed) ~= length(y_observed))
error('ERROR: Both X and Y vectors must be the same size!');
end
%Calculate vector length once
vector_length = length(x_observed);
%
% Theta = [ x1 1 ]
% [ x2 1 ]
% [ x3 1 ]
% [ x4 1 ]
% [ ... ... ]
%
Theta = ones(vector_length, 2);
Theta(1:end, 1) = x_observed;
%
% Theta_0 = [ y1 ]
% [ y2 ]
% [ y3 ]
% [ y4 ]
% [ ... ]
%
Theta_0 = y_observed;
%
% Theta_LS = [ M ]
% [ B ]
%
Theta_LS = ((Theta' * Theta)^-1 ) * Theta' * Theta_0;
%Return the answer
m = Theta_LS(1);
b = Theta_LS(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
