Friday, August 19, 2016

matlab code to study upsampling(Sampling rate) of audio file

%% matlab code to study upsampling(Sampling rate) of audio file
% sounds like a problem of asynchronous sample rate conversion. 
% To convert from one sample rate to another: 
% use sinc interpolation to compute the continuous time representation 
% of the signal then resample at our new sample rate.
% resample the signal to have sample times that are fixed.
% Read in the sound
% returns the sample rate(FS) in Hertz
% N = number sample 
% the number of bits per sample (BITS) used to encode the data in the file
% Minh Anh Nguyen
% minhanhnguyen@q.com

% close
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables.
clc;% Clear the command window.
imtool close all; % Close all imtool figures.
workspace; % Make sure the workspace panel is showing.
%% read and plot .wav file
figure;

myvoice1 = audioread('good2.wav');
% Reads in the sound file, into a big array called y.
[y1, fs1] = audioread('good2.wav');
size(y1);
left1=y1(:,1);
right1=y1(:,2);
N = length(right1);
% Normalize y; that is, scale all values to its maximum. 
y = right1/max(abs(right1));
t1 = (0:1:N-1)/fs1;
plot(t1,right1 );
str1=sprintf('Sound with sampling rate = %d Hz and number sample = %d', fs1, N);
title(str1);
xlabel('time (sec)'); ylabel('relative signal strength'); grid on;
axis tight
grid on;






%% Compute the power spectral density, a measurement of the energy at various frequencies
% DFT to describe the signal in the frequency
NFFT = 2 ^ nextpow2(N);
Y = fft(y1, NFFT) / N;
f = (fs1 / 2 * linspace(0, 1, NFFT / 2+1))'; % Vector containing frequencies in Hz
amp = ( 2 * abs(Y(1: NFFT / 2+1))); % Vector containing corresponding amplitudes
figure;
plot (f, amp); xlim([0,1600]);  ylim([0,0.5e-4]);
title ('plot single-sided amplitude spectrume of signal')
xlabel ('frequency (Hz)')
ylabel ('|y(f)|')
grid on;

% display markers at specific data points on DFT graph 
indexmin1 = find(min(amp) == amp);
xmin1 = f(indexmin1);
ymin1 = amp(indexmin1);
indexmax1 = find(max(amp) == amp);
xmax1 = f(indexmax1);
ymax1 = amp(indexmax1);
strmax1 = [' Maximum = ',num2str(xmax1),' Hz','    ', num2str(ymax1),' dB'];
text(xmax1,ymax1,strmax1,'HorizontalAlignment','Left');





fsorg=44100;
fs = 8000;
%%  resample is your function. To downsample signal from 44100 Hz to 8000 Hz:
%(the "1" and "2" arguments define the resampling ratio: 8000/44100 = 1/2)
%To upsample back to 44100 Hz: x2 = resample(y,2,1);
figure;
y_1 = resample(y1,fs,fsorg);
N1_1 = length(y_1);

t1_1 = (0:1:length(y_1)-1)/fs;
plot(t1_1, y_1);
str1=sprintf('Sound with sampling rate = %d Hz and number sample = %d', fs, N1_1);
title(str1);
xlabel('time (sec)');
ylabel('signal strength')
axis tight
grid;



%% Compute the power spectral density, a measurement of the energy at various frequencies
% DFT to describe the signal in the frequency
NFFTd = 2 ^ nextpow2(N1_1);
Yd = fft(y_1, NFFTd) / N1_1;
fd = (fs / 2 * linspace(0, 1, NFFTd / 2+1))'; % Vector containing frequencies in Hz
amp_down = ( 2 * abs(Yd(1: NFFTd / 2+1))); % Vector containing corresponding amplitudes
figure;
plot (fd, amp_down); xlim([0,1600]);  ylim([0,0.5e-4]);
title ('\t plot single-sided amplitude spectrume of downsampling signal')
xlabel ('frequency (Hz)')
ylabel ('|y(f)|')
grid on;

% display markers at specific data points on DFT graph 
indexmin1 = find(min(amp_down) == amp_down);
xmin1 = fd(indexmin1);
ymin1 = amp_down(indexmin1);
indexmax1 = find(max(amp_down) == amp_down);
xmax1 = fd(indexmax1);
ymax1 = amp_down(indexmax1);
strmax1 = [' Maximum = ',num2str(xmax1),' Hz','    ', num2str(ymax1),' dB'];
text(xmax1,ymax1,strmax1,'HorizontalAlignment','Left');



Matlab code to compute the corresponding absorption coefficients

 %% oxyhemoglobin and deoxyhemoglobin extinction.m
 % Author:Minh Anh Nguyen
 % minhanhnguyen@q.com
 % Matlab code to compute the corresponding absorption coefficients and plot 
 % the three absorption spectra on the same graph. 
 % Identify the low-absorption near-IR window that provide deep
 % penetration.
 % data for molar extinction coefficients of oxy-and deoxyhemoglobin and 
 % absorption coefficient of pure water as a function of wavelength are
 % copied directly from this website: http://omlc.org/spectra/hemoglobin/summary.html
 % Use physiologically representative values for both oxygen saturation SO2
 % and total concentration of hemoglobin CHb.
 % These values for the molar extinction coefficient e in [cm-1/(moles/liter)] were compiled by Scott Prahl using data from

    %W. B. Gratzer, Med. Res. Council Labs, Holly Hill, London
    %N. Kollias, Wellman Laboratories, Harvard Medical School, Boston 

%To convert this data to absorbance A, multiply by the molar concentration and the pathlength. For example, if x is the number of grams per liter and a 1 cm cuvette is being used, then the absorbance is given by

        %(e) [(1/cm)/(moles/liter)] (x) [g/liter] (1) [cm]
  %A =  ---------------------------------------------------
                     %     64,500 [g/mole]

%using 64,500 as the gram molecular weight of hemoglobin.

%To convert this data to absorption coefficient in (cm-1), multiply by the molar concentration and 2.303,

   % µa = (2.303) e (x g/liter)/(64,500 g Hb/mole) 
%where x is the number of grams per liter. A typical value of x for whole blood is x=150 g Hb/liter. 


close all; clear; clc;
load ('oxy_deoxy.mat');
figure;
semilogy(oxy_deoxy(:,1), oxy_deoxy(:,2),  'b', 'linewidth',1.5);
hold on
semilogy(oxy_deoxy(:,1), oxy_deoxy(:,3),  'r', 'linewidth',1.5);
title('The molar oxyhemoglobin and deoxyhemoglobin extinction','FontSize',14);
set(gca,'FontSize',14);
ylabel('molar extinction coefficient (cm^{-1}(moles/l)^{-1})','fontsize',14);
xlabel('wavelength (nm)','fontsize',14);
axis([250 1000 0 10^6]);
set(gcf,'PaperUnits','inches','PaperPosition',[0 0 6 5]);
grid on;
legend('y = HbO2 ','y = Hb')



%% purewater

load ('purewater.mat');
figure;
semilogy(water(:,1), water(:,2),  'b', 'linewidth',1.5);
title('The pure water specific absorption coefficient ','FontSize',14);
set(gca,'FontSize',14);
ylabel(' specific absorption coefficient (um^{-1})','fontsize',14);
xlabel('wavelength (nm)','fontsize',14);
axis([300 1000 0 1*10^6]);
set(gcf,'PaperUnits','inches','PaperPosition',[0 0 6 5]);
grid on;
legend('y = pure water')





%% all three graph
figure;
semilogy(oxy_deoxy(:,1), oxy_deoxy(:,3)*0.0054,  'r', 'linewidth',1.5);
hold on
semilogy(oxy_deoxy(:,1), oxy_deoxy(:,2),  'b', 'linewidth',1.5);
hold on
semilogy(water(:,1), water(:,2),  'y', 'linewidth',1.5);
title('The molar hemoglobin and water extinction','FontSize',14);
set(gca,'FontSize',14);
ylabel('molar extinction coefficient (cm^{-1}(moles/l)^{-1})','fontsize',14);
xlabel('wavelength (nm)','fontsize',14);
%axis([250 1000 0 0.6*10^6]);
axis([250 1000 0 1*10^6]);
set(gcf,'PaperUnits','inches','PaperPosition',[0 0 6 5]);
grid on;
legend('y = HbO2 ','y = Hb', 'y = pure water')



%% Absorption
figure;
semilogy(oxy_deoxy(:,1), oxy_deoxy(:,2)*0.0054,  'b', 'linewidth',1.5);
hold on
semilogy(oxy_deoxy(:,1), oxy_deoxy(:,3)*0.0054,  'r', 'linewidth',1.5);
hold on
semilogy(water(:,1), water(:,2),  'y', 'linewidth',1.5);
title('The hemoglobin and water','FontSize',14);
set(gca,'FontSize',14);
ylabel('effective absorption coefficient (um^{-1})','fontsize',14);
xlabel('wavelength (nm)','fontsize',14);
axis([250 1000 0 1*10^6]);
set(gcf,'PaperUnits','inches','PaperPosition',[0 0 6 5]);
grid on;
legend('y = HbO2 ','y = Hb', 'y = pure water')



Thursday, August 18, 2016

Matlab code to calculate and plot pulses on a lossless transmission line at13.5Mhz

% this matlab code is calculate and plot pulses on a lossless transmission line
% Frequency is 13.5Mhz
%minhanhnguyen@q.com

% Clear the screen
clc;
clear;
% close all windows
close all;

% use these Parameters to calculate and plot Pulse Propagation on a Transline
Ro=50; %Ro = line impedance
Rl=100; %Rl = load resistance
Rg=25; %Rg = generator resistance
Vo=10; %Vo = voltage source
T=.08; %T = pulse width
v=1; %v = wave velocity
d=1; %d = line length
dt=.02; %dt = time step,
n=170; %n = number of time steps

% Voltage Divider at Input
vDivider = (Ro/(Ro+Rg));

% Reflection Coefficient at Load
gamLoad = (Rl-Ro)/(Rl+Ro);

% Reflection Coefficient at Generator
gamGenerator = (Rg-Ro)/(Rg+Ro);

% Transit Time
tTime = d/v;

% Position along the line
x=[0:d/500:d];

for j = 0:n;

% Time
t = j*dt;

% Leading edge of pulse
xl = v*t;

% Trailing edge of pulse
xt = v*(t-T);

% Voltage along the line at time t
V = Vo*(vDivider*(x<xl).*(x>xt));
V = V + Vo*(vDivider*gamLoad)*(x>(2*d-xl)).*(x<(2*d-xt));
V = V + Vo*(vDivider*gamLoad*gamGenerator)*(x<(xl-2*d)).*(x>(xt-2*d));
V = V + Vo*(vDivider*(gamLoad^2)*gamGenerator)*(x>(4*d-xl)).*(x<(4*d-xt));

% Plotting the result
plot(x,V);axis([0 d -Vo Vo]);grid;
title('Pulse voltage vs. distance');
xlabel('Distance Along The Line');
ylabel('Pulse Voltage');
pause(.5)
end



Matlab code to calculate the voltage on the transmission line


% This matlab code is setup for calculating the voltage on the 
% transmission line. The code will then plot the results. 
% frequency is 13.5Mhz


% Clear the screen
clc;
clear all;
% close all windows
close all;


%% input values 
RL = 20;        % resistive part of load
CL = 3.3e-9;    % capacitive part of load
len = 80;       % transmission line length
freq = 13.5e+6;  % frequency
omega = 2*pi*freq;
Z0 = 50;        % characteristic impedance of line
v = 2e+8;       % velocity in RG 58 cable
ZS = 50;        % source impedance

%% set up the node measurements with arbitrary values

z_node = [0:4:len];      % equivalent position of nodes on transmission line

V_node = [0 1 2 3.1 4 5 6 7 8.2 9 10 11.4 12 13 14 15 16 17];
V_node = [V_node 18.1 19 20.5];

%%set up phase measurements with arbitrary values
phase_node = [0 1 2 3.1 4 5 6 7 8.2 9 10 11.4 12 13 14 15 16 17];
phase_node = [phase_node  18.1 19 20.5];


%%calculating the transmission line voltage 

z = [0:1:len];                  % z position on line
ZL = RL - j./omega./CL;         % load impedance
V_line = ZL .* z .* exp(-j.*omega.*z./v) ./len;


%% Find the magnitude and angle of the voltage phasor
V_line_mag = abs(V_line);
V_line_deg = 180*angle(V_line)./pi;


%% Plot results
figure(1);
subplot(2,1,1), plot(z,V_line_mag, z_node, V_node, 'o');
title('Unmodified line vs. V magnitude');
ylabel('V magnitude');
legend('trans. line analysis','circuit analysis', 0);
subplot(2,1,2), plot(z,V_line_deg, z_node, phase_node, 'o');
title('Unmodified line vs. phase');
ylabel('phase(degrees)');
xlabel('z position on line (meters)');
legend('trans. line analysis','circuit analysis', 0);


MatLab function to plot Amplitude Envelopes of .wav file

% This MatLab function takes a "wav" sound file as an input
% perform the average and peak envelopes of an input signal,
% and then plot all three plots on a graph.
% usage: amplitude_envelopes ('*.wav')
% Author: Minh Anh Nguyen
% email: minhanhnguyen@q.com


function [amplitude_envelopes] = amplitude_envelopes ( file );
close all;% Close all figures (except those of imtool.)
 clf % Clears the graphic screen.


% Reads in the sound file, into a big array called y.
y1=audioread(file);
[y1, fs1]=audioread(file);
size(y1);
left1=y1(:,1);
right1=y1(:,2);

% Normalize y; that is, scale all values to its maximum. 
y = right1/max(abs(right1));

% If you want to hear the sound after you read it in 
% Sound just plays a file to the speaker.

sound(y, 44100);

% Initialize peak and average arrays to be the same as the signal. We’re 
% keeping three big arrays, all of the same length.

average_env = y;
max_env = y;

% Set the window lengths directly, in number of samples.
% the average should be smaller than the peak, since it will tend to 
% flatten out if it gets too big. 
% average and peak window size are two numbers, which can be played with 
% to vary what the pictures will look like. 
% Note: use two different window sizes so that we can have 
% different kinds of resolution for peak and average charts.

average_window_size = 512;
peak_window_size = 10000;

% from the input signal, taking an average of the previous, 
% average window size number of samples, and store that in the current 
% place in the average array.
% Use a loop that starts at the end of the first window and goes to the 
% end of the sound file, indicated by length(average_env). 
% use "sum" command to take some range of a vector.

for k = average_window_size:length(average_env)
running_sum = sum(abs(y(k-average_window_size+1:k)));
average_env(k) = running_sum /average_window_size ;
end

% from the input signal, taking the maximum value of the previous 
% peak_window_size number of samples.

for k = peak_window_size:length(max_env)
max_env(k)= max(abs(y(k-peak_window_size+1:k)));

end

% Plot the three "signals": the original in the color cyan, the 
% peak in magenta, and the average in blue. The colors are the 
% letters in single quotes.

plot(y, 'c')
hold on
plot (max_env, 'm')
hold on
plot(average_env, 'b')
title('Monochord: Signal, Average Signal Envelope, Peak Signal Envelope')
xlabel('sample number')
ylabel('amplitude (-1 to 1)')
legend('original signal', 'peak envelope', 'running average envelope', 0)