To view the presentation, please click on the below link:
http://www.slideshare.net/minh65/electrocardiogram-ecg-or-ekg
Study biomedical signals and images, Matlab, and LabView code
Saturday, August 27, 2016
Sunday, August 21, 2016
To view the presentation, please click the below link:
http://www.slideshare.net/minh65/simulation-results-of-induction-heating-coil
http://www.slideshare.net/minh65/simulation-results-of-induction-heating-coil
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');
% 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')
% 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
% 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)
% 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)
Subscribe to:
Posts (Atom)










