From 4a2cc8ff7154fcbb4d6e68102e32b7720dec2acb Mon Sep 17 00:00:00 2001
From: xenakios <xenakios@gmail.com>
Date: Sun, 3 Dec 2017 18:39:47 +0200
Subject: [PATCH] Initial work to separate spectrum processing implementations
 into stand alone functions from the ProcessedStretch class, to allow
 visualizations etc without having to instantiate the full ProcessedStretch
 object

---
 Source/PS_Source/ProcessedStretch.cpp |   9 --
 Source/PS_Source/ProcessedStretch.h   | 148 ++++++++++++++++++++++++++
 2 files changed, 148 insertions(+), 9 deletions(-)

diff --git a/Source/PS_Source/ProcessedStretch.cpp b/Source/PS_Source/ProcessedStretch.cpp
index 13291ba..3a26779 100644
--- a/Source/PS_Source/ProcessedStretch.cpp
+++ b/Source/PS_Source/ProcessedStretch.cpp
@@ -174,15 +174,6 @@ void ProcessedStretch::process_spectrum(REALTYPE *freq)
 //void ProcessedStretch::process_output(REALTYPE *smps,int nsmps){
 //};
 
-
-REALTYPE profile(REALTYPE fi, REALTYPE bwi){
-	REALTYPE x=fi/bwi;
-	x*=x;
-	if (x>14.71280603) return 0.0;
-	return exp(-x);///bwi;
-
-};
-
 void ProcessedStretch::do_harmonics(REALTYPE *freq1,REALTYPE *freq2){
 	REALTYPE freq=pars.harmonics.freq;
 	REALTYPE bandwidth=pars.harmonics.bandwidth;
diff --git a/Source/PS_Source/ProcessedStretch.h b/Source/PS_Source/ProcessedStretch.h
index 056b194..3a12640 100644
--- a/Source/PS_Source/ProcessedStretch.h
+++ b/Source/PS_Source/ProcessedStretch.h
@@ -180,6 +180,154 @@ struct ProcessParameters
 	}
 };
 
+inline REALTYPE profile(REALTYPE fi, REALTYPE bwi) {
+	REALTYPE x = fi / bwi;
+	x *= x;
+	if (x>14.71280603) return 0.0;
+	return exp(-x);///bwi;
+
+};
+
+inline void spectrum_do_harmonics(ProcessParameters& pars, std::vector<REALTYPE>& tmpfreq1, int nfreq, double samplerate, REALTYPE *freq1, REALTYPE *freq2) {
+	REALTYPE freq = pars.harmonics.freq;
+	REALTYPE bandwidth = pars.harmonics.bandwidth;
+	int nharmonics = pars.harmonics.nharmonics;
+
+	if (freq<10.0) freq = 10.0;
+
+	REALTYPE *amp = tmpfreq1.data();
+	for (int i = 0; i<nfreq; i++) amp[i] = 0.0;
+
+	for (int nh = 1; nh <= nharmonics; nh++) {//for each harmonic
+		REALTYPE bw_Hz;//bandwidth of the current harmonic measured in Hz
+		REALTYPE bwi;
+		REALTYPE fi;
+		REALTYPE f = nh * freq;
+
+		if (f >= samplerate / 2) break;
+
+		bw_Hz = (pow(2.0f, bandwidth / 1200.0f) - 1.0f)*f;
+		bwi = bw_Hz / (2.0f*samplerate);
+		fi = f / samplerate;
+
+		REALTYPE sum = 0.0f;
+		REALTYPE max = 0.0f;
+		for (int i = 1; i<nfreq; i++) {//todo: optimize here
+			REALTYPE hprofile;
+			hprofile = profile((i / (REALTYPE)nfreq*0.5f) - fi, bwi);
+			amp[i] += hprofile;
+			if (max<hprofile) max = hprofile;
+			sum += hprofile;
+		};
+	};
+
+	REALTYPE max = 0.0;
+	for (int i = 1; i<nfreq; i++) {
+		if (amp[i]>max) max = amp[i];
+	};
+	if (max<1e-8f) max = 1e-8f;
+
+	for (int i = 1; i<nfreq; i++) {
+		//REALTYPE c,s;
+		REALTYPE a = amp[i] / max;
+		if (!pars.harmonics.gauss) a = (a<0.368f ? 0.0f : 1.0f);
+		freq2[i] = freq1[i] * a;
+	};
+
+};
+
+inline void spectrum_add(int nfreq, REALTYPE *freq2, REALTYPE *freq1, REALTYPE a) {
+	for (int i = 0; i<nfreq; i++) freq2[i] += freq1[i] * a;
+};
+
+inline void spectrum_zero(int nfreq,REALTYPE *freq1) {
+	for (int i = 0; i<nfreq; i++) freq1[i] = 0.0;
+};
+
+inline void spectrum_do_freq_shift(ProcessParameters& pars, int nfreq, double samplerate, REALTYPE *freq1, REALTYPE *freq2) {
+	spectrum_zero(nfreq, freq2);
+	int ifreq = (int)(pars.freq_shift.Hz / (samplerate*0.5)*nfreq);
+	for (int i = 0; i<nfreq; i++) {
+		int i2 = ifreq + i;
+		if ((i2>0) && (i2<nfreq)) freq2[i2] = freq1[i];
+	};
+};
+
+inline void spectrum_do_pitch_shift(ProcessParameters& pars, int nfreq, REALTYPE *freq1, REALTYPE *freq2, REALTYPE _rap) {
+	spectrum_zero(nfreq,freq2);
+	if (_rap<1.0) {//down
+		for (int i = 0; i<nfreq; i++) {
+			int i2 = (int)(i*_rap);
+			if (i2 >= nfreq) break;
+			freq2[i2] += freq1[i];
+		};
+	};
+	if (_rap >= 1.0) {//up
+		_rap = 1.0f / _rap;
+		for (int i = 0; i<nfreq; i++) {
+			freq2[i] = freq1[(int)(i*_rap)];
+		};
+
+	};
+};
+
+inline void spectrum_do_octave(ProcessParameters& pars, int nfreq, double samplerate, 
+	std::vector<REALTYPE>& sumfreq, 
+	std::vector<REALTYPE>& tmpfreq1,
+	REALTYPE *freq1, REALTYPE *freq2) {
+	spectrum_zero(nfreq,sumfreq.data());
+	if (pars.octave.om2>1e-3) {
+		spectrum_do_pitch_shift(pars,nfreq, freq1, tmpfreq1.data(), 0.25);
+		spectrum_add(nfreq, sumfreq.data(), tmpfreq1.data(), pars.octave.om2);
+	};
+	if (pars.octave.om1>1e-3) {
+		spectrum_do_pitch_shift(pars,nfreq, freq1, tmpfreq1.data(), 0.5);
+		spectrum_add(nfreq,sumfreq.data(), tmpfreq1.data(), pars.octave.om1);
+	};
+	if (pars.octave.o0>1e-3) {
+		spectrum_add(nfreq,sumfreq.data(), freq1, pars.octave.o0);
+	};
+	if (pars.octave.o1>1e-3) {
+		spectrum_do_pitch_shift(pars,nfreq, freq1, tmpfreq1.data(), 2.0);
+		spectrum_add(nfreq,sumfreq.data(), tmpfreq1.data(), pars.octave.o1);
+	};
+	if (pars.octave.o15>1e-3) {
+		spectrum_do_pitch_shift(pars,nfreq, freq1, tmpfreq1.data(), 3.0);
+		spectrum_add(nfreq,sumfreq.data(), tmpfreq1.data(), pars.octave.o15);
+	};
+	if (pars.octave.o2>1e-3) {
+		spectrum_do_pitch_shift(pars, nfreq, freq1, tmpfreq1.data(), 4.0);
+		spectrum_add(nfreq,sumfreq.data(), tmpfreq1.data(), pars.octave.o2);
+	};
+
+	REALTYPE sum = 0.01f + pars.octave.om2 + pars.octave.om1 + pars.octave.o0 + pars.octave.o1 + pars.octave.o15 + pars.octave.o2;
+	if (sum<0.5f) sum = 0.5f;
+	for (int i = 0; i<nfreq; i++) freq2[i] = sumfreq[i] / sum;
+};
+
+inline void spectrum_do_filter(ProcessParameters& pars, int nfreq, double samplerate, REALTYPE *freq1, REALTYPE *freq2) {
+	REALTYPE low = 0, high = 0;
+	if (pars.filter.low<pars.filter.high) {//sort the low/high freqs
+		low = pars.filter.low;
+		high = pars.filter.high;
+	}
+	else {
+		high = pars.filter.low;
+		low = pars.filter.high;
+	};
+	int ilow = (int)(low / samplerate * nfreq*2.0f);
+	int ihigh = (int)(high / samplerate * nfreq*2.0f);
+	REALTYPE dmp = 1.0;
+	REALTYPE dmprap = 1.0f - pow(pars.filter.hdamp*0.5f, 4.0f);
+	for (int i = 0; i<nfreq; i++) {
+		REALTYPE a = 0.0f;
+		if ((i >= ilow) && (i<ihigh)) a = 1.0f;
+		if (pars.filter.stop) a = 1.0f - a;
+		freq2[i] = freq1[i] * a*dmp;
+		dmp *= dmprap + 1e-8f;
+	};
+};
+
 class SpectrumProcess
 {
 public: