-
Notifications
You must be signed in to change notification settings - Fork 0
/
note_decompose.m
100 lines (79 loc) · 1.76 KB
/
note_decompose.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
function note_decompose(filename, n)
close all;
[y,Fs]=audioread(filename);
sound(y,Fs);
l=length(y);
%compute fft
Y=fft(y);
P2 = abs(Y/l);
P1 = P2(1:l/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(l/2))/l;
P2_cplx = (Y/l);
P1_cplx = P2_cplx(1:l/2+1);
P1_cplx(2:end-1) = 2*P1_cplx(2:end-1);
%get fundamental frequency
[M,I]= max(P1);
f0=f(I);
%plot period in time domain
y_time=y(1:round(2*2*Fs/f0));
t = 0:1/Fs:(length(y_time)/Fs)-1/Fs;
subplot(2, 2, 1);
plot(t,y_time);
title('Time Domain')
xlabel('Seconds');
ylabel('Amplitude');
subplot(2, 2, 2);
plot(f,P1)
title('FFT')
xlabel('f (Hz)')
ylabel('Magnitude')
xlim([20 20000]) %audible range
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Decompose n Dominate sinusoids
[pks,locs] = findpeaks(P1);
pks_sort=pks;
locs_sort=locs;
for i=1:length(pks_sort)-1
for j=1:length(pks_sort)-1
if pks_sort(j)<pks_sort(j+1)
pksTemp=pks_sort(j);
pks_sort(j)=pks_sort(j+1);
pks_sort(j+1)=pksTemp;
locsTemp=locs_sort(j);
locs_sort(j)=locs_sort(j+1);
locs_sort(j+1)=locsTemp;
end
end
end
%addititive synthesis
tGen = 0:1/Fs:(length(y)/Fs)-1/Fs;
yGen=zeros(1,length(tGen));
if n>length(pks_sort)
n=length(pks_sort);
end
for i=1:n
phase=angle(P1_cplx(locs_sort(i)))
mag=pks_sort(i);
yGen=yGen-mag*sin(2*pi*f(locs_sort(i))*tGen-phase);
end
sound(yGen,Fs);
YGen=fft(yGen);
P2 = abs(YGen/l);
P1 = P2(1:l/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(l/2))/l;
%plot stuff
y_time=yGen(1:round(2*2*Fs/f0));
t = 0:1/Fs:(length(y_time)/Fs)-1/Fs;
subplot(2, 2, 3);
plot(t,y_time);
title('Time Domain')
xlabel('Seconds');
ylabel('Amplitude');
subplot(2, 2, 4);
plot(f,P1)
title('FFT')
xlabel('f (Hz)')
ylabel('Magnitude')
xlim([20 20000]) %audible range