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')