My problems starts out with a variable length of samples. Usually, it is 1024 or higher powers of 2. The DFT of this "signal" is taken and only the amplitude spectrum is retained and the phase information is ignored. Now, it is my job to define a phase for this amplitude to get back a complex fourier transform such that it's inverse transform yields a signal which is within a predefined maximum and minimum. If we were talking about audio(wav) data, the range for each sample would be from -1 to 1 or -0.5 to 0.5 in some cases.
Here is what I have tried so far- As a prototype algorithm, I tried working with 8 samples at a time. I generated 8 inequalities for each of the samples and tried working with it.
s = [3 1 4 1 5 9 2 6];
#s is the original 8 samples
#s = rand(1,8).-0.5
#s = [1 2 3 4 5 6 7 8]
f = fft(s);
#f is the dft
a = abs(f);
#a is the amplitude spectrum
an = arg(f);
#p is the phase information to be ignored
p = an;
#vat1 is the value of the signal at time 1
vat1 = 0.25*(a(2)*cos(p(2)) + a(3)*cos(p(3)) + a(4)*cos(p(4))) + (1/8)*(a(1) + a(5)*cos(p(5)))
#vat2 is the value of the signal at time 2
vat2 = 0.25*(a(2)*(0.7071*cos(p(2)) - 0.7071*sin(p(2))) - a(3)*sin(p(3)) - a(4)*(0.7071*cos(p(4)) + 0.7071*sin(p(4)))) + 0.125*(a(1) - a(5)*cos(p(5)))
#vat3 is the value of the signal at time 3
vat3 = 0.25*(-a(2)*sin(p(2)) - a(3)*cos(p(3)) + a(4)*(sin(p(4)))) + 0.125*(a(1) + a(5)*cos(p(5)))
#vat4 is the value of the signal at time 4
vat4 = 0.25*(a(2)*(-0.7071*cos(p(2)) - 0.7071*sin(p(2))) + a(3)*sin(p(3)) + a(4)*(0.7071*cos(p(4)) - 0.7071*sin(p(4)))) + 0.125*(a(1) - a(5)*cos(p(5)))
#vat5 is the value of the signal at time 5
vat5 = 0.25*(a(2)*(-cos(p(2))) + a(3)*cos(p(3)) - a(4)*cos(p(4))) + (1/8)*(a(1) + a(5)*cos(p(5)))
#vat6 is the value of the signal at time 6
vat6 = 0.25*(a(2)*(-0.7071*cos(p(2)) + 0.7071*sin(p((2)))) - a(3)*(sin(p(3))) + a(4)*(0.7071*cos(p(4)) + 0.7071*sin(p(4)))) + 0.125*(a(1) - a(5)*cos(p(5)))
#vat7 is the value of the signal at time 7
vat7 = 0.25*(a(2)*(sin(p(2))) + a(3)*(-cos(p(3))) + a(4)*(-sin(p(4)))) + 0.125*(a(1) + a(5)*cos(p(5)))
#vat 8 is the value of the signal at time 8
vat8 = 0.25*(a(2)*(0.7071*cos(p(2)) + 0.7071*sin(p(2))) + a(3)*(sin(p(3))) + a(4)*(-0.7071*cos(p(4)) + 0.7071*sin(p(4)))) + 0.125*(a(1) - a(5)*cos(p(5)))
The "vat" variables are all expected to be within those predefined bounds mentioned earlier. I tried solving them as a system of trigonometric inequalities, but that yielded good results very rarely. This was because 8 samples provided amplitudes such that simply giving random phases were sufficient to keep the signal within the required range. Extending this to 1024 samples would be great, but that is too tedious to solve.
My current program to perform this for 8 samples is as follows
aud = audioread('president.wav');
aud = aud(:,1);
aud = aud';
wc = 49848;
#Audio file has been converted to a single row array
a = [];#a is the amplitude spectrum
f = [];#f is the fourier transform
p = [];#p is the phase spectrum
t = [];#t is the temporary array
s = [];#s is the reconstructed signal array
l = length(aud);
l/8
fflush(stdout);
r = rem(l,8);
aud = aud(1:l-r);
l = length(aud);
for i = 1:(l/8)
t = [t;(aud((i*8)-7:(i*8)))];
end
display("Here");
fflush(stdout);
for i = 1:rows(t)
i
fflush(stdout)
f = [f ; fft(t(i,:))];
end
p = arg(f);
a = abs(f);
na = a;
sf = [];
sig = zeros(1,length(aud));
p1 = 0;
p2 = 0;
p3 = 0;
p4 = 0;
p5 = 0;
p6 = 0;
p7 = 0;
p8 = 0;
ai = [];
p24a = [];
la = [];#array for l values
ha = [];#array for h values
p3s = [];#The cos p3 values
aac15 = [];
aac15n = [];
lla = [];
ula = [];
l = min(aud);
h = max(aud);
p24bca = [];
p3la = [];
p3ha = [];
p4la = [];
p4ha = [];
for i = 1:rows(a)
i
fflush(stdout);
ac = a(i,:);#Current a values
pc = p(i,:);#Current p values. Can be ignored
h = sum(ac)/8;#h is the highest possible value that can be generated with the given amplitude spectrum
l = -h;#l is the lowest possible value
la = [la;l];#Just an array for all l values
ha = [ha;h];#Just an array for all h values
ac15 = ac(1)-(ac(5)*cos(p(i)));#ac15 is defined like this.
ac15n = ac(1)+(ac(5)*cos(p(i)));
aac15 = [aac15 ; ac15];
aac15n = [aac15n ; ac15n];
ll = (ac15n - (8*l))/(2*ac(3));
ul = (ac15n - (8*h))/(2*ac(3));
lla = [lla ; ll];
ula = [ula ; ul];
p3h = (ac15n - 8*l)/(2*ac(3));
p3l = (ac15n - 8*h)/(2*ac(3));
#np3l = (ac15n - 8*l)/(2);
#np3h = (ac15n - 8*h)/(2);
p3la = [p3la ; p3l];
p3ha = [p3ha ; p3h];
#p3s = [p3s;(2*ac15n - 8*(l+h))/(4*ac(3))];
#p3 = pc(3);
if p3l<-1 || p3h>1
p3l = -1;
p3h = 1;
wc--;
end
p3 = acos(p3l + (p3h - p3l)*rand())
p2l = real(((2*l*(sqrt(2)-1)) - 0.5*((ac15/sqrt(2)) + (ac15n/2)) - (2*ac(3)*cos(p3)))/ac(2));
p2h = real(((2*h*(sqrt(2)-1)) - 0.5*((ac15/sqrt(2)) + (ac15n/2)) - (2*ac(3)*cos(p3)))/ac(4));
if p2l<-1 || p2h>1
p2l = -1;
p2h = 1;
wc--;
end
p2 = acos(p2l + (p2h-p2l)*rand(1,1));
p4l = real(((4*l) - 0.5*(ac15n) - ac(2)*cos(p2) - ac(3)*cos(p3))/ac(4));
p4h = real(((4*h) - 0.5*(ac15n) - ac(2)*cos(p2) - ac(3)*cos(p3))/ac(4));
p4la = [p4la ; p4l];
p4ha = [p4ha ; p4h];
if p4l<-1 || p4h>1
p4l = -1;
p4h = 1;
wc--;
end
p4 = acos(p4l + (p4h-p4l)*rand(1,1));
#p4 = pc(4);
#p2 = pc(2);
#p2 = pc(2);
#p4 = pc(4);
# p2 = p24/2;
# p4 = p24/2;
#p4 = p2;#Just for now. Change it later
p6 = -p4;
p7 = -p3;
p8 = -p2;
p5 = p(i);
ip = [0 p2 p3 p4 p5 p6 p7 p8];
ip = real(ip);
ip(isinf(ip)) = 0;
ip(isnan(ip)) = 0;
ip;
fflush(stdout)
ai = [ai ; ip];
sf(i,1) = complex(ac(1)*cos(pc(1)),ac(1)*sin(pc(1)));
sf(i,2) = complex(ac(2)*cos(ip(2)),ac(2)*sin(ip(2)));
sf(i,3) = complex(ac(3)*cos(ip(3)),ac(3)*sin(ip(3)));
sf(i,4) = complex(ac(4)*cos(ip(4)),ac(4)*sin(ip(4)));
sf(i,5) = complex(ac(5)*cos(ip(5)),ac(5)*sin(ip(5)));
sf(i,6) = complex(ac(6)*cos(ip(6)),ac(6)*sin(ip(6)));
sf(i,7) = complex(ac(7)*cos(ip(7)),ac(7)*sin(ip(7)));
sf(i,8) = complex(ac(8)*cos(ip(8)),ac(8)*sin(ip(8)));
end
sig = [];
for i = 1:rows(sf)
i
fflush(stdout);
s(i,:) = ifft(sf(i,:));
sig = [sig s(i,:)];
end
So basically, the question boils down to this- what phase will keep the signal bound within an upper limit and a lower limit?
