# OH - 2018 - Subsetsum implementation to reduce large datasets
# Using  third order intermodulation products

# based on "A general algorithm to calcuate third order intermodulation product locations for any number of tones" (2015)
# by: Chris Arnott

library(stats)
library(pracma)


plot.fourier <- function(fourier.series, f.0, ts) {
  w <- 2*pi*f.0
  trajectory <- sapply(ts, function(t) fourier.series(t,w))
  plot(ts, trajectory, type="l", xlab="time", ylab="f(t)"); abline(h=0,lty=3)
}
plot.frequency.spectrum <- function(X.k, xlimits=c(0,length(X.k))) {
  plot.data  <- cbind(0:(length(X.k)-1), Mod(X.k))

  # TODO: why this scaling is necessary?
  plot.data[2:length(X.k),2] <- 2*plot.data[2:length(X.k),2] 
  
  plot(plot.data, t="h", lwd=2, main="", 
       xlab="Frequency (Hz)", ylab="Strength", 
       xlim=xlimits, ylim=c(0,max(Mod(plot.data[,2]))))
}


#f.sum   <- function(t,w) { 
#  dc.component + 
#  sum( component.strength * sin(component.freqs*w*t + component.delay)) 
#}
f.sum   <- function(t,w) { 
  sum(  sin(component.freqs*w*t )) 
}
 

readkey <- function()
{
    cat ("Press [enter] to continue")
    line <- readline()
}


# MAIN
n      <-  4
orig   <-  c(1,2,3,4)  
print("original set:",cat(sets))
targetsum <- 7
f1 <- 1212.5
fch <- 60

sets <- f1 + fch*orig
targetfreq <-  f1 + fch*target

 

# move components to frequency domain to reduce the list in "sets"
acq.freq <- 100                    # base data acquisition (sample) frequency (Hz)
time     <- 60                      # measuring time interval (seconds)
ts       <- seq(0,time-1/acq.freq,1/acq.freq) # vector of sampling time-points (s) 
f.0      <- 1/time
component.freqs <- sets                        # frequency of signal components (Hz)

# how this looks in a plot
w <- 2*pi*f.0
trajectory <- sapply(ts, function(t) f.sum(t,w))

# Intermodulation - third order product
trajectory<-trajectory**3

X.k <- fft(trajectory)                   # find all harmonics with fft()
plot.frequency.spectrum(Mod(X.k), xlimits=c(0,acq.freq*time/2))


#plot(trajectory,type="l")

# frequency location of  sum distortion products
#find peaks

spectrum <- Mod(X.k)[0:5000]
#plot(spectrum,type='l')
p<-findpeaks(spectrum, npeaks=50, threshold=4, sortstr=TRUE)
points(p[, 2], p[, 1], pch=20, col="maroon")


peaks <- p[,2][order(p[,2])]
peaks
sets