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);
% Filename: Beat_Frequency.m
%
% [Richard Lyons, Feb. 2013]
clear, clc
Fs = 8192; % Sample rate of dig. samples
N = 8192; % Number of time samples
n = 0:N-1;
Wave_1 = sin(2*pi*210*n/Fs); % First tone, 210 Hz
Wave_2 = sin(2*pi*200*n/Fs); % Second tone, 200 Hz
% Plot the two tones
figure(1)
plot(n/Fs, Wave_1, '-b')
ylabel('200 Hz'); xlabel('Time (sec.)');
hold on
plot(n/Fs, Wave_2, '-r')
axis([0, 0.05, -1.2, 1.5]);
ylabel('Input tones'); xlabel('Time (sec.)');
title('red = 200 Hz tone, blue = 210 Hz tone');
grid on, zoom on
hold off
Product = Wave_1.*Wave_2;
Sum = Wave_1 + Wave_2;
% Plot the tones' product and sum
figure(2)
subplot(2,1,1)
plot(n/Fs, Product, '-b'),
ylabel('Product'); xlabel('Time (sec.)');
grid on, zoom on
hold on
Red_Curve = 0.5*cos(2*pi*10*n/Fs) + 0.5; % Used for plotting only
plot(n/Fs, Red_Curve, '-r')
axis([0, 0.3, -1.25, 1.5]);
hold off
grid on, zoom on
subplot(2,1,2)
plot(n/Fs, Sum, '-b')
hold on
Red_Curve = 2*cos(2*pi*5*n/Fs); % Used for plotting only
plot(n/Fs, Red_Curve, '-r')
axis([0, 0.3, -2.4, 3]);
hold off
ylabel('Sum'); xlabel('Time (sec.)');
grid on, zoom on
% Play all the signals
sound(Wave_1, Fs)
pause(1.2)
sound(Wave_2, Fs)
pause(1.2)
sound(Product, Fs)
pause(1.2)
sound(Sum, Fs)
% Spec analysis of the "Sum" signal
Spec = fft(Sum);
Spec_Mag = abs(Spec);
Freq = n*Fs/N; % Freq vector in Hz
figure (3) % Plot positive-freq spec amg
plot(Freq(1:N/16), Spec_Mag(1:N/16))
title('Spec Mag of "Sum" signal')
ylabel('Mag.'), xlabel('Hz')
grid on, zoom on
//Determination of spectrum of a signal
//With maximum normalized frequency f = 0.4
//using Blackmann window
clear all;
close;
clc;
N = 11;
cfreq = [0.4 0];
[wft,wfm,fr]=wfir('lp',N,cfreq,'re',0);
wft; // Time domain filter coefficients
wfm; // Frequency domain filter values
fr; // Frequency sample points
for n = 1:N
h_balckmann(n)=0.42-0.5*cos(2*%pi*n/(N-1))+0.08*cos(4*%pi*n/(N-1));
wft_blmn(n) = wft(n)*h_balckmann(n);
end
wfm_blmn = frmag(wft_blmn,length(fr));
WFM_blmn_dB =20*log10(wfm_blmn);
plot2d(fr,WFM_blmn_dB)
xtitle('Frequency Response of Blackmann window Filtered output N = 11','Frequency in cycles per samples f','Energy density in dB')
clf;
L = 100;
for (L_factor = 1:1:6)
%L_factor = 3;
N_factor = 1;
L = round(L*L_factor);
n = [0:1:L-1];
x = sin(pi/1000*(n.*n));
subplot(3,1,1);
stem(n,x);
title ('x[n]');
subplot(3,1,2);
y = fft(x,L*N_factor);
plot(abs(y));
title ('DFT');
subplot(3,1,3);
xr = ifft(y);
stem(n,xr(1:L));
title ('IDFT');
pause;
end
//Caption: Determination of Frequency Resolution of the
//[1]. Barlett [2]. Welch [3].Blackman Tukey Power Spectrum Estimation
//Methods
clear;
clc;
close;
Q = 10; //Quality factor
N = 1000; //Length of the sample sequence
//Bartlett Method
F_Bartlett = Q/(1.11*N);
disp(F_Bartlett,'Frequency Resolution of Bartlett Power Spectrum Estimation')
//Welch Method
F_Welch = Q/(1.39*N);
disp(F_Welch,'Frequency Resolution of Welch Power Spectrum Estimation')
//Blackmann-Tukey Method
F_Blackmann_Tukey = Q/(2.34*N);
disp(F_Blackmann_Tukey,'Frequency Resolution of Blackmann Tukey Power Spectrum Estimation')
//Result
//Frequency Resolution of Bartlett Power Spectrum Estimation
// 0.0090090
//Frequency Resolution of Welch Power Spectrum Estimation
// 0.0071942
//Frequency Resolution of Blackmann Tukey Power Spectrum Estimation
// 0.0042735
function [b,a] = rc_filter(R, C, Fs, filter_type)
% Returns equivalent IIR coefficients for an analog RC filter
%
% Usage: [B,A] = RC_FILTER(r, c, fs, type);
%
% R is the resistance value (in ohms)
% C is the capacitance value (in farrads)
% FS is the digital sample rate (in Hz)
% type is a character string defining filter type
% Choices are: 'high' or 'low'
%
% Highpass filters have the analog structure:
%
%
% | |
% Vi o--------| |----------+---------o Vo
% | | |
% C |
% ---
% | | R
% | |
% ---
% |
% |
% o---------------------+---------o
% GND
%
%
% Lowpass filters have the analog structure:
%
%
% |-----|
% Vi o--------| |------+---------o Vo
% |-----| |
% R |
% ----- C
% -----
% |
% |
% o---------------------+---------o
% GND
%
% This function uses a pre-calculated equation for both of these circuits
% that only requires the resistance and capacitance value to get a true
% digital filter equivalent to a basic analog filter.
%
% The math behind these equations is based off the basic bilinear transform
% technique that can be found in many DSP textbooks. The reference paper
% for this function was "Conversion of Analog to Digital Transfer
% Functions" by C. Sidney Burrus, page 6.
%
% This is also the equivalent of a 1st order butterworth with a cuttoff
% frequency of Fc = 1/(2*pi*R*C);
%
% Author: sparafucile17 07/01/02
%
% Verify that cutoff of the analog filter is below Nyquist
if( (1/(2*pi*R*C)) > (Fs/2) )
error('This analog filter cannot be realized with this sample rate');
end
% Default to allpass if invalid type is selected
b = [ 1 0 ];
a = [ 1 0 ];
% Constants
RC = R * C;
T = 1 / Fs;
% Analog Cutoff Fc
w = 1 / (RC);
% Prewarped coefficient for Bilinear transform
A = 1 / (tan((w*T) / 2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following equations were derived from
%
% s
% T(s) = -------
% s + 1
%
%
% using Bilinear transform of
%
% 1 ( 1 - z^-1 )
% s --> ----------- * ------------
% tan(w*T/2) ( 1 + z^-1 )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(strcmp(filter_type,'high'))
b(1) = (A) / (1 + A);
b(2) = -b(1);
a(2) = (1 - A) / (1 + A);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following equations were derived from
%
% 1
% T(s) = -------
% s + 1
%
%
% using Bilinear transform of
%
% 1 ( 1 - z^-1 )
% s --> ----------- * ------------
% tan(w*T/2) ( 1 + z^-1 )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(strcmp(filter_type,'low'))
b(1) = (1) / (1 + A);
b(2) = b(1);
a(2) = (1 - A) / (1 + A);
end
function [y] = signal_unmute(x, index, duration, Fs)
% Return the input vector with a unmute that occurs at a specific
% index and with an exponential ramp-up to reduce "pop" sounds.
% The output will start off at -100dB gain and end at 0dB gain.
%
% Usage: y = SIGNAL_UNMUTE(x, index, duration, Fs);
%
% X is your one-dimensional input array
% INDEX is where in the input signal you want the unmute to begin
% DURATION is how long (in seconds) to exponentially ramp up the
% input signal. 100ms is recommended.
% FS is the sample rate
%
% Example:
% You want to unmute your signal at 0.5sec with a duration of 100ms
% and with a 48k Fs. (signal is unmuted completely 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 unmute. ' ...
'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; %What do we consider "muted" to be in decibels
decayrate = -dB_atten / duration;
coeff = 1.0 - 10^(-decayrate/(20.0*Fs));
% Build the Gain array
gain = [zeros(index, 1); ones((length(x)-index), 1);];
b = [ coeff 0];
a = [ 1 (-1*(1-coeff))];
gain = filter(b,a,gain);
% Apply Mute (gain) to the input signal
y = gain .* x;
% Flip the vector (if required)
if(flip_back == true);
y = y';
end
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
function thr = sureshrink(CD,T)
%function used to calculate threshold using sureshrink method
CD = CD(:)';
n = length(CD);
sx2 = sort(abs(CD)-T).^2; % sorting in descending order
b = cumsum(sx2); %cumulative sum
risks = (n-(2*(1:n))+b)/n;
[risk,best] = min(risks);
thr = sqrt(sx2(best));
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;
