This is Pvc_Object.m in view mode; [Download] [Up]
#import <sys/message.h>
#import <mach_error.h>
#import <soundkit/Sound.h>
#import <string.h>
#import "Pvc_Object.h"
void bitreverse( float *buf,int N );
@implementation Pvc_Object : Object
{
/*
int N,
N2,
Nw = 0,
Nw2,
decim = 0,
interp = 0,
in,
on,
valid = -1,
eof = 0;
float freqmlt = 0.,
srate = 44100.,
pi,
twopi,
synt = 0.,
*Hwin,
*Wanal,
*Wsyn,
*input,
*buf,
*channel,
*output;
BOOL obank = 0,
aflag = 0,
sflag = 0;
NXZone *ozone;
*/
}
+ create {
id newInstance;
newInstance = [ self new ];
return newInstance;
}
- init {
[super init];
[self setDefaults];
inputFilename = malloc(1024);
status = IDLE;
return self;
}
- setInputFilename: (char *) arg
{
int err;
if (inputFilename && !(strcmp(arg,inputFilename)))
return self;
if (insfh != NULL)
SNDFree(insfh);
if (strlen(arg) > 1023)
return nil;
strcpy(inputFilename,arg);
err = SNDReadSoundfile(inputFilename, &insfh);
if (err != SND_ERR_NONE)
return nil;
return self;
}
- setOutputFilename: (char *) arg
{
outputFilename = arg;
return self;
}
static BOOL isPowerOfTwo(int n)
/* Query whether n is a pure power of 2 */
{
while (n > 1) {
if (n % 2) break;
n >>= 1;
}
return (n == 1);
}
- checkArgs:(char *)msg
/* If failure, sets msg to error message and returns nil.
Msg should be at least 1024 long */
{
if (Nw == 0)
Nw = N;
if (interp == 0)
interp = decim;
obank = freqmlt != 0.;
N2 = N>>1;
Nw2 = Nw>>1;
if (interp > Nw) {
sprintf(msg,"I must be less than or equal to window size.");
return nil;
}
if (decim > Nw) {
sprintf(msg,"D must be less than or equal to window size.");
return nil;
}
if (decim > Nw) {
sprintf(msg,"D must be less than or equal to window size.");
return nil;
}
if (!isPowerOfTwo(N)) {
sprintf(msg,"FFT size (N) must be a power of 2.");
}
return self;
}
- allocateDataspace {
int width,
size;
if ( ozone != NULL )
NXDestroyZone(ozone);
ozone = NXCreateZone( 10 * Nw * sizeof(float), vm_page_size, 0 );
/* analysis window */
Wanal = (float *) zspace( ozone, Nw, sizeof(float) );
/* synthesis window */
Wsyn = (float *) zspace( ozone, Nw, sizeof(float) );
/* input buffer */
input = (float *) zspace( ozone, Nw, sizeof(float) );
/* hamming window */
Hwin = (float *) zspace( ozone, Nw, sizeof(float) );
/* FFT buffer */
buf = (float *) zspace( ozone, N, sizeof(float) );
/* analysis channels */
channel = (float *) zspace( ozone, N+2, sizeof(float) );
/* output buffer */
output = (float *) zspace( ozone, Nw, sizeof(float) );
SNDGetDataPointer( insfh, (char **)&inData, &insize, &width );
if (outsfh != NULL)
SNDFree(outsfh);
SNDAlloc( &outsfh,
( (int) ((N+insize) * width * ( (float) interp / decim ))),
SND_FORMAT_LINEAR_16, insfh->samplingRate,
insfh->channelCount, 0);
SNDGetDataPointer(outsfh, (char **)&outData, &size, &width );
srate = insfh->samplingRate;
return self;
}
- (SNDSoundStruct *)outsfh
{
return outsfh;
}
- (SNDSoundStruct *)insfh
{
return insfh;
}
- setDefaults {
N = 1024;
Nw = 0;
valid = -1;
decim = 128;
interp = 0;
eof = 0;
freqmlt = 0.;
srate = 44100.;
pi = 4 * atan(1.);
twopi = 8 * atan(1.);
synt = 0.;
obank = 0;
inPoint = 0;
outPoint = 0;
return self;
}
- (float)howFar {
return ((float)outPoint) / ((outsfh->dataSize)/sizeof(short));
}
- N: (int) aN {
N = aN;
return self;
}
- Nw: (int) aNw {
Nw = aNw;
return self;
}
- decim: (int) aDecim {
decim = aDecim;
return self;
}
- interp: (int) aInterp {
interp = aInterp;
return self;
}
- freqmlt: (float) aFreqmlt {
freqmlt = aFreqmlt;
return self;
}
- srate: (float) aSrate {
srate = aSrate;
return self;
}
- synt: (float) aSynt {
synt = aSynt;
return self;
}
- aflag: (BOOL) aAflag {
aflag = aAflag;
return self;
}
- sflag: (BOOL) aSflag {
sflag = aSflag;
return self;
}
- writeOutputSound {
int err;
err = SNDWriteSoundfile( outputFilename, outsfh );
if (err != SND_ERR_NONE)
return nil;
return self;
}
- kill
{
[self resume];
eof = 1;
return self;
}
- resume
/* Gets called from the other thread */
{
kern_return_t ec;
msg_header_t msg = {0, /* msg_unused */
TRUE, /* msg_simple */
sizeof(msg_header_t),/* msg_size */
MSG_TYPE_NORMAL, /* msg_type */
0}; /* Fills in remaining fields */
if (eof != 2 && appToObjPort != PORT_NULL)
return nil;
msg.msg_local_port = PORT_NULL;
msg.msg_remote_port = appToObjPort;
ec = msg_send(&msg, SEND_TIMEOUT, 0);
if (ec == SEND_TIMED_OUT)
; /* message queue is full, don't need to send another */
return self;
}
- pause
{
eof = 2;
return self;
}
- (int)status
{
return status;
}
- _waitForMessage
/* This is the main loop for a separate-threaded performance. */
{
struct {
msg_header_t header;
char data[MSG_SIZE_MAX];
} msg;
status = PAUSED;
(void)port_allocate(task_self(), &appToObjPort);
msg.header.msg_size = MSG_SIZE_MAX;
msg.header.msg_local_port = appToObjPort;
(void)msg_receive(&msg.header, MSG_OPTION_NONE,0);
(void)port_deallocate(task_self(),appToObjPort);
appToObjPort = PORT_NULL;
status = RUNNING;
return self;
}
/* phase vocoder */
- runPvc
/* Returns self if successfully writes file, else nil. */
{
status = RUNNING;
/* initialize input and output time values (in samples) */
in = -Nw;
if ( decim )
on = (in * interp) / decim;
else
on = in;
/* main loop--perform phase vocoder analysis-resynthesis */
while ( eof != 1) {
if (eof == 2) {
[self _waitForMessage];
if (eof == 1)
break;
}
/* increment times */
in += decim;
on += interp;
/* analysis: input D samples; window, fold and rotate input
samples into FFT buf; take FFT; and convert to
amplitude-frequency (phase vocoder) form */
eof = [self shiftin];
[self fold];
[self rfft:YES];
[self convert];
/* at this point channel[2*i] contains amplitude data and
channel[2*i+1] contains frequency data (in Hz) for phase
vocoder channels i = 0, 1, ... N/2; the center frequency
associated with each channel is i*f, where f is the
fundamental frequency of analysis R/N; any desired spectral
modifications can be made at this point: pitch modifications
are generally well suited to oscillator bank resynthesis,
while time modifications are generally well (and more
efficiently) suited to overlap-add resynthesis */
/* oscillator bank resynthesis */
if ( obank ) {
[self oscbank];
/* oscbank( channel, N2, srate, Nw, interp, freqmlt, output ); */
[self shiftout];
}
/* overlap-add resynthesis */
else {
[self unconvert];
[self rfft:NO];
[self overlapadd];
[self shiftout];
}
}
if (![self writeOutputSound]) {
status = IDLE;
return nil;
}
status = IDLE;
return self;
}
/* overlap-add FFT analysis/resynthesis without phase unwrapping */
- runFFT
{
/* argument checking */
status = RUNNING;
/* initialize input and output time values (in samples) */
in = -Nw;
if ( decim )
on = (in * interp) / decim;
else
on = in;
/* main loop--perform phase vocoder analysis-resynthesis */
while ( !eof ) {
/* increment times */
in += decim;
on += interp;
/* analysis: input D samples; window, fold and rotate input
samples into FFT buf; take FFT; and convert to
amplitude-frequency (phase vocoder) form */
eof = [self shiftin];
[self fold];
[self rfft:YES];
[self convertFFT];
/* at this point channel[2*i] contains amplitude data and
channel[2*i+1] contains frequency data (in Hz) for phase
vocoder channels i = 0, 1, ... N/2; the center frequency
associated with each channel is i*f, where f is the
fundamental frequency of analysis R/N; any desired spectral
modifications can be made at this point: pitch modifications
are generally well suited to oscillator bank resynthesis,
while time modifications are generally well (and more
efficiently) suited to overlap-add resynthesis */
/* overlap-add resynthesis */
[self unconvertFFT];
[self rfft:NO];
[self overlapadd];
[self shiftout];
}
status = IDLE;
return self;
}
/* shift next D samples into righthand end of array A of
length N, padding with zeros after last sample (A is
assumed to be initially 0); return 0 when more input
remains, otherwise return 1 after N-2*D zeros have been
padded onto end of input */
- (int) shiftin
{
register int i;
if ( valid < 0 ) /* first time only */
valid = Nw;
for ( i = 0 ; i < Nw - decim ; i++ )
*(input+i) = *(input + i + decim);
if ( valid == Nw ) {
for ( i = (Nw - decim); i < Nw; i++) {
if ( inPoint >= insize ) {
*(input+i) = (float) (*(inData + inPoint) / 32768.);
valid = i - (Nw-decim);
break;
}
*(input+i) = (float) (*(inData + inPoint) / 32768.);
++inPoint;
}
}
if ( valid < Nw ) { /* pad with zeros after EOF */
for (i=valid; i < Nw; i++)
*(input+i) = 0.;
valid -= decim;
}
return (valid <= 0);
}
/* if output time n >= 0, output first I samples in
array A of length N, then shift A left by I samples,
padding with zeros after last sample */
- shiftout
{
int i;
if ( on >= 0 ) {
for ( i=0; i < interp; i++ )
*(outData + outPoint + i) = (short) (32767. * *(output+i));
outPoint += interp;
}
for ( i = 0; i < Nw - interp; i++ )
*(output+i) = *(output + i + interp);
for ( i = (Nw - interp); i < Nw; i++ )
*(output+i) = 0.;
return self;
}
/* S is a spectrum in rfft format, i.e., it contains N real values
arranged as real followed by imaginary values, except for first
two values, which are real parts of 0 and Nyquist frequencies;
convert first changes these into N/2+1 PAIRS of magnitude and
phase values to be stored in output array C; the phases are then
unwrapped and successive phase differences are used to compute
estimates of the instantaneous frequencies for each phase vocoder
analysis channel; decimation rate D and sampling rate R are used
to render these frequency values directly in Hz. */
- convert
{
static int first = 1;
static float *lastphase,
fundamental,
factor;
float phase,
phasediff;
int real,
imag,
amp,
freq;
float a,
b;
int i;
/* first pass: allocate zeroed space for previous phase
values for each channel and compute constants */
if ( first ) {
first = 0;
lastphase = (float *) space( N2+1, sizeof(float) );
bzero(lastphase,(N2+1)*sizeof(float));
fundamental = srate / (float) (N2<<1);
factor = srate / (decim * twopi);
}
/* unravel rfft-format spectrum: note that N2+1 pairs of
values are produced */
for ( i = 0; i <= N2; i++ ) {
imag = freq = ( real = amp = i<<1 ) + 1;
a = ( i == N2 ? *(buf+1) : *(buf+real) );
b = ( i == 0 || i == N2 ? 0. : *(buf+imag) );
/* compute magnitude value from real and imaginary parts */
*(channel+amp) = hypot( a, b );
/* compute phase value from real and imaginary parts and take
difference between this and previous value for each channel */
if ( *(channel+amp) == 0. )
phasediff = 0.;
else {
phasediff = ( phase = -atan2( b, a ) ) - *(lastphase+i);
*(lastphase+i) = phase;
/* unwrap phase differences */
while ( phasediff > pi )
phasediff -= twopi;
while ( phasediff < -pi )
phasediff += twopi;
}
/* convert each phase difference to Hz */
*(channel+freq) = (phasediff * factor) + (i * fundamental);
}
return self;
}
- convertFFT
{
int amp,
phase;
float a,b;
register int i;
for ( i = 0; i <= N2; i++ ) {
phase = (amp = i<<1) + 1;
a = ( i == N2 ? *(buf+1) : *(buf+amp) );
b = ( i == 0 || i == N2 ? 0. : *(buf+phase) );
/* compute magnitude value from real and imaginary parts */
*(channel+amp) = hypot( a, b );
*(channel+phase) = -atan2( b, a );
}
return self;
}
/* unconvert essentially undoes what convert does, i.e., it
turns N2+1 PAIRS of amplitude and frequency values in
C into N2 PAIR of complex spectrum data (in rfft format)
in output array S; sampling rate R and interpolation factor
I are used to recompute phase values from frequencies */
- unconvert
{
static int first = 1;
static float *lastphase,
fundamental,
factor;
int i,
real,
imag,
amp,
freq;
float phase;
/* first pass: allocate memory and compute constants */
if ( first ) {
first = 0;
lastphase = (float *) space( N2+1, sizeof(float) );
fundamental = srate / (float) (N2<<1);
factor = twopi * interp / srate;
}
/* subtract out frequencies associated with each channel,
compute phases in terms of radians per I samples, and
convert to complex form */
for ( i = 0; i <= N2; i++ ) {
imag = freq = ( real = amp = i<<1 ) + 1;
if ( i == N2 )
real = 1;
*(lastphase+i) += *(channel+freq) - i * fundamental;
phase = *(lastphase+i) * factor;
*(buf+real) = *(channel+amp) * cos( phase );
if ( i != N2 )
*(buf+imag) = - *(channel+amp) * sin( phase );
}
return self;
}
- unconvertFFT
{
int amp,
phase;
register int i;
for ( i = 0; i <= N2; i++ ) {
phase = (amp = i<<1) + 1;
if ( i == N2 )
*(buf+1) = *(channel+amp) * cos( *(channel+phase) );
else
*(buf+amp) = *(channel+amp) * cos( *(channel+phase) );
if ( i != N2 )
*(buf+phase) = - *(channel+amp) * sin( *(channel+phase) );
}
return self;
}
/* make balanced pair of analysis (A) and synthesis (S) windows;
window lengths are Nw, FFT length is N, synthesis interpolation
factor is I, and osc is true (1) if oscillator bank resynthesis
is specified */
- makewindows
{
int i;
float sum;
/* basic Hamming windows */
for ( i = 0; i < Nw; i++ )
*(Hwin+i) = *(Wanal+i) = *(Wsyn+i) = 0.54 - (0.46 * cos(
twopi * i / (Nw - 1) ));
/* when Nw > N, also apply interpolating (sinc) windows to
ensure that window are 0 at increments of N (the FFT length)
away from the center of the analysis window and of I away
from the center of the synthesis window */
if ( Nw > N ) {
float x;
/* take care to create symmetrical windows */
x = -(Nw - 1) / 2.;
for ( i = 0; i < Nw; i++, x += 1. )
if ( x != 0. ) {
*(Wanal+i) *= N * sin( (pi * x) / N ) / (pi*x);
if ( interp )
*(Wsyn+i) *= interp * sin( (pi*x) / interp ) / (pi*x);
}
}
/* normalize windows for unity gain across unmodified
analysis-synthesis procedure */
for ( sum = i = 0; i < Nw; i++ )
sum += *(Wanal+i);
for ( i = 0; i < Nw; i++ ) {
float afac = 2./sum;
float sfac = Nw > N ? 1./afac : afac;
*(Wanal+i) *= afac;
*(Wsyn+i) *= sfac;
}
if ( Nw <= N && interp ) {
for ( sum = i = 0; i < Nw; i += interp )
sum += *(Wsyn+i) * *(Wsyn+i);
for ( sum = 1./sum, i = 0; i < Nw; i++ )
*(Wsyn+i) *= sum;
}
return self;
}
/* input I is a folded spectrum of length N; output O and
synthesis window W are of length Nw--overlap-add windowed,
unrotated, unfolded input data into output O */
- overlapadd
{
int i,
n;
n = on;
while ( n < 0 )
n += N;
n %= N;
for ( i = 0 ; i < Nw ; i++ ) {
*(output+i) += *(buf+n) * *(Wsyn+i);
if ( ++n == N )
n = 0;
}
return self;
}
/* multiply current input I by window W (both of length Nw);
using modulus arithmetic, fold and rotate windowed input
into output array O of (FFT) length N according to current
input time n */
- fold
{
int i,
n;
n = in;
for ( i = 0; i < N; i++ )
*(buf+i) = 0.;
while ( n < 0 )
n += N;
n %= N;
for ( i = 0; i < Nw; i++ ) {
*(buf+n) += *(input+i) * *(Wanal+i);
if ( ++n == N )
n = 0;
}
return self;
}
/* If forward is true, rfft replaces 2*N real data points in buf with
N complex values representing the positive frequency half of their
Fourier spectrum, with *(buf+1) replaced with the real part of the Nyquist
frequency value. If forward is false, rfft expects buf to contain a
positive frequency spectrum arranged as before, and replaces it with
2*N real values. N MUST be a power of 2. */
- rfft: (BOOL) forward
{
float c1,c2,
h1r,h1i,
h2r,h2i,
wr,wi,
wpr,wpi,
temp,
theta;
float br,bi;
int i,
i1,i2,i3,i4,
N2p1;
theta = pi/N2;
wr = 1.;
wi = 0.;
c1 = 0.5;
if ( forward ) {
c2 = -0.5;
[self cfft:forward];
br = *buf;
bi = *(buf+1);
}
else {
c2 = 0.5;
theta = -theta;
br = *(buf+1);
bi = 0.;
*(buf+1) = 0.;
}
wpr = -2. * pow( sin( 0.5*theta ), 2. );
wpi = sin(theta);
N2p1 = (N2 <<1) + 1;
for ( i = 0; i <= N2 >> 1; i++ ) {
i1 = i<<1;
i2 = i1 + 1;
i3 = N2p1 - i2;
i4 = i3 + 1;
if ( i == 0 ) {
h1r = c1 * ( *(buf+i1) + br );
h1i = c1 * ( *(buf+i2) - bi );
h2r = -c2 * ( *(buf+i2) + bi );
h2i = c2 * ( *(buf+i1) - br );
*(buf+i1) = h1r + (wr * h2r) - (wi * h2i);
*(buf+i2) = h1i + (wr * h2i) + (wi * h2r);
br = h1r - (wr * h2r) + (wi * h2i);
bi = -h1i + (wr * h2i) + (wi * h2r);
}
else {
h1r = c1 * ( *(buf+i1) + *(buf+i3) );
h1i = c1 * ( *(buf+i2) - *(buf+i4) );
h2r = -c2 * ( *(buf+i2) + *(buf+i4) );
h2i = c2 * ( *(buf+i1) - *(buf+i3) );
*(buf+i1) = h1r + wr*h2r - wi*h2i;
*(buf+i2) = h1i + wr*h2i + wi*h2r;
*(buf+i3) = h1r - wr*h2r + wi*h2i;
*(buf+i4) = -h1i + wr*h2i + wi*h2r;
}
wr = ((temp = wr) * wpr) - (wi * wpi) + wr;
wi = (wi * wpr) + (temp * wpi) + wi;
}
if ( forward )
*(buf+1) = br;
else
[self cfft:forward];
return self;
}
/* cfft replaces float array x containing NC complex values
(2*NC float values alternating real, imagininary, etc.)
by its Fourier transform if forward is true, or by its
inverse Fourier transform if forward is false, using a
recursive Fast Fourier transform method due to Danielson
and Lanczos. NC MUST be a power of 2. */
- cfft: (BOOL) forward
{
float wr,wi,
wpr,wpi,
theta,
scale;
int mmax,
ND,
m,
i,j,
delta;
ND = N2<<1;
bitreverse( buf, ND );
for ( mmax = 2; mmax < ND; mmax = delta ) {
delta = mmax<<1;
theta = twopi / ( forward? mmax : -mmax );
wpr = -2. * pow( sin( 0.5*theta ), 2. );
wpi = sin(theta);
wr = 1.;
wi = 0.;
for ( m = 0; m < mmax; m += 2 ) {
register float rtemp, itemp;
for ( i = m; i < ND; i += delta ) {
j = i + mmax;
rtemp = ( wr * *(buf+j) ) - ( wi * *(buf+j+1) );
itemp = ( wr * *(buf+j+1) ) + ( wi * *(buf+j) );
*(buf+j) = *(buf+i) - rtemp;
*(buf+j+1) = *(buf+i+1) - itemp;
*(buf+i) += rtemp;
*(buf+i+1) += itemp;
}
wr = ((rtemp = wr) * wpr) - (wi * wpi) + wr;
wi = (wi * wpr) + (rtemp * wpi) + wi;
}
}
/* scale output */
scale = forward ? 1./ND : 2.;
{
register float *bi=buf, *be=buf + ND;
while ( bi < be )
*bi++ *= scale;
}
return self;
}
/* bitreverse places float array x containing N/2 complex values
into bit-reversed order */
void bitreverse( float *buf,int N )
{
float rtemp,itemp;
int i,j,
m;
for ( i = j = 0; i < N; i += 2, j += m ) {
if ( j > i ) {
rtemp = *(buf+j); /* complex exchange */
itemp = *(buf+j+1);
*(buf+j) = *(buf+i);
*(buf+j+1) = *(buf+i+1);
*(buf+i) = rtemp;
*(buf+i+1) = itemp;
}
for ( m = N>>1; m >= 2 && j >= m; m >>= 1 )
j -= m;
}
}
/* oscillator bank resynthesizer for phase vocoder analyzer
uses sum of N+1 cosinusoidal table lookup oscillators to
compute I (interpolation factor) samples of output O
from N+1 amplitude and frequency value-pairs in C;
frequencies are scaled by P */
- oscbank
{
static int NP,
tablesize = 8192,
first = 1;
static float Iinv,
*lastamp,
*lastfreq,
*index,
*table,
Pinc,
ffac;
int amp,
freq,
n,
chan;
/* first pass: allocate memory to hold previous values
of amplitude and frequency for each channel, the table
index for each oscillator, and the table itself; also
compute constants */
if ( first ) {
float twopioL = twopi / tablesize, tabscale;
first = 0;
lastamp = (float *) space( N2+1, sizeof(float) );
lastfreq = (float *) space( N2+1, sizeof(float) );
index = (float *) space( N2+1, sizeof(float) );
table = (float *) space( tablesize+1, sizeof(float) );
tabscale = ( Nw >= N2 ? N2 : 8 * N2 );
for ( n = 0; n < tablesize+1; n++ )
*(table+n) = tabscale*cos( twopioL*n );
Iinv = 1. / interp;
Pinc = freqmlt * ( (float) tablesize ) / srate;
ffac = freqmlt * pi / ((float) N2);
if ( freqmlt > 1. )
NP = (int) ((float) N2 / freqmlt);
else
NP = N2;
}
/* for each channel, compute I samples using linear
interpolation on the amplitude and frequency
control values */
for ( chan=0; chan < NP; chan++ ) {
register float a,
ainc,
f,
finc,
address;
freq = ( amp = ( chan << 1 ) ) + 1;
if ( *(channel+amp) < synt ) /* skip the little ones */
continue;
*(channel+freq) *= Pinc;
finc = ( *(channel+freq) - ( f = *(lastfreq+chan) ) ) *Iinv;
ainc = ( *(channel+amp) - ( a = *(lastamp+chan) ) ) *Iinv;
address = *(index+chan);
/* accumulate the I samples from each oscillator into
output array O (initially assumed to be zero);
f is frequency in Hz scaled by oscillator increment
factor and pitch (Pinc); a is amplitude; */
for ( n=0; n < interp; n++ ) {
*(output+n) += a * *(table + ( (int) address ));
address += f;
while ( address >= (float) tablesize )
address -= (float) tablesize;
while ( address < 0. )
address += (float) tablesize;
a += ainc;
f += finc;
}
/* save current values for next iteration */
*(lastfreq+chan) = *(channel+freq);
*(lastamp+chan) = *(channel+amp);
*(index+chan) = address;
}
return self;
}
@end
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.