rt300@4: //
rt300@4: //  frequencer.cpp
rt300@4: //  oscSenderExample
rt300@4: //
rt300@4: //  Created by Robert Tubb on 27/11/2012.
rt300@4: //
rt300@4: //
rt300@4: 
rt300@4: #include "frequencer.h"
rt300@4: Frequencer frequencer;
rt300@4: 
rt300@4: Frequencer::Frequencer(){
rt300@4:     
rt300@4:     N = 16;
rt300@4: 
rt300@4:     for(int i = 0; i<N;i++){
rt300@4:         timePoints.push_back(0.0);
rt300@4:         freqMag.push_back(0.0);
rt300@4:         freqPhase.push_back(0.0);
rt300@4:     }
rt300@4:     
rt300@4: }
rt300@4: /////////////////////////////////////////////////////
rt300@4: /////////////////////////////////////////////////////
rt300@9: // No time domain!
rt300@4: void Frequencer::timeEdit(int n, double value){
rt300@4: 
rt300@4:     // update internal
rt300@4:     timePoints[n] = value;
rt300@4:     
rt300@4:     // calculate transform
rt300@4:     dft();
rt300@4: 
rt300@4: }
rt300@4: 
rt300@4: void Frequencer::timeEdit(vector<double> allValues){
rt300@4:     
rt300@4:     // update internal
rt300@4:     vector<double>::iterator iter;
rt300@4:     vector<double>::iterator iter2;
rt300@4:     iter = allValues.begin();
rt300@4:     iter2 = timePoints.begin();
rt300@4:     for (iter = allValues.begin(); iter < allValues.end(); iter++,iter2++){
rt300@4:         *iter = *iter2;
rt300@4:     }
rt300@4:     
rt300@4:     
rt300@4:     // calculate transform
rt300@4:     dft();
rt300@4:     
rt300@4:     //TODO  output freq domain - to PD
rt300@4:     
rt300@4: }
rt300@4: 
rt300@4: /////////////////////////////////////////////////////
rt300@4: /////////////////////////////////////////////////////
rt300@4: 
rt300@4: 
rt300@4: vector<double> Frequencer::freqMagEdit(vector<int> indexes, vector<double> values){
rt300@4:     if(indexes.size() != values.size()){
rt300@4:         cout << "freqMagEdit error: indexes dont match values\n";
rt300@4:         return timePoints;
rt300@4:     }
rt300@4:     // put into local store,
rt300@4:     for(int i=0;i<indexes.size();i++){
rt300@4:         freqMag[indexes[i]] = values[i];
rt300@20:         if(i!=0)freqMag[N-indexes[i]] = freqMag[indexes[i]];
rt300@4:     }
rt300@4:    
rt300@4:     // calculate transform
rt300@4:     idft();
rt300@4: 
rt300@4:     return timePoints;
rt300@4: }
rt300@4: 
rt300@4: 
rt300@4: /////////////////////////////////////////////////////
rt300@4: /////////////////////////////////////////////////////
rt300@4: 
rt300@4: void Frequencer::freqPhaseEdit(int k, double phase){
rt300@4:     N = 16;
rt300@4:     // put into local store
rt300@4:     freqPhase[k] = phase;
rt300@4:     
rt300@4:     // making sure to conjugate
rt300@20:     if(k!=0)freqPhase[N-k] = - freqPhase[k];
rt300@4:     // calculate transform
rt300@4:     idft();
rt300@4:     // TODO output time domain
rt300@4:    
rt300@4: }
rt300@4: 
rt300@4: /////////////////////////////////////////////////////
rt300@4: /////////////////////////////////////////////////////
rt300@4: 
rt300@4: void Frequencer::dft(){
rt300@4:     
rt300@4:     // we have 16 time points. work out freqs and phases
rt300@4:     int n = 0;
rt300@4:     int k = 0;
rt300@4:     double xre = 0.;
rt300@4:     double xim = 0.;
rt300@4:     
rt300@4:     for(k=0; k < N; k++){
rt300@4:         // work out xn
rt300@4:         xre = 0.;
rt300@4:         xim = 0.;
rt300@4:         for(n=0;n < N;n++){
rt300@4:             xre = xre + timePoints[n]*cos( (2.*PI/(float)N)*n*k);
rt300@4:             xim = xim + timePoints[n]*sin( (2.*PI/(float)N)*n*k);
rt300@4:             
rt300@4:         }
rt300@4:         freqMag[k] = sqrt(pow(xre,2) + pow(xim,2.)); // MAGNITUDE
rt300@4:         freqPhase[k] = -atan2(xim,xre);
rt300@4:     }
rt300@4:     
rt300@4: }
rt300@4: 
rt300@4: 
rt300@4: 
rt300@4: /////////////////////////////////////////////////////
rt300@4: /////////////////////////////////////////////////////
rt300@4: 
rt300@4: void Frequencer::idft(){
rt300@4:     
rt300@4:     // we have some freqs and phases. work out real time points (conj)
rt300@4:     int n = 0;
rt300@4:     int k = 0;
rt300@4:     double xre = 0.;
rt300@4:     double xim = 0.;
rt300@4:     
rt300@4:     for(n=0; n < N; n++){
rt300@4:         // work out xn
rt300@4:         xre = 0.;
rt300@4:         xim = 0.;
rt300@4:         for(k=0;k < N;k++){
rt300@4:             xre = xre + freqMag[k]*cos( (2.*PI/(float)N)*n*k + freqPhase[k]);
rt300@4:             
rt300@4:         }
rt300@4:         timePoints[n] = xre/(float)N; // divided by N ???
rt300@4:     }
rt300@4:     
rt300@4: }
rt300@4: /////////////////////////////////////////////////////
rt300@4: /////////////////////////////////////////////////////