Harmonics of a vibrating string.

Grain · 2005-08-29T07:21:02

How can I determine which harmonics will be present and in what amplitudes when a string is plucked at a given point along the string. For example let’s say I pluck the string in the exact middle I want to find out what will be amplitude of the fundamental and the subsequent harmonics?

Math and Physics Programming

Started by Grain August 25, 2005 12:43 AM

20 comments, last by Winograd 18 years, 7 months ago

Grain

500

Author

August 29, 2005 02:05 AM

Ok I put together a very basic waveguide just to get accustomed to the concept.

Let me know if I got the basic concept. From what I read I under stand the Delay line to be a sequential container containing the wave’s amplitude at multiple points along the string. Is this correct.

One thing I noticed right away was the shape the was tracing out wasn’t a sinusoid like here, but rather a parallelogram.

You can just copy and compile this, it’s a complete program.

#include <windows.h>#include <cmath>#include <gl/gl.h>/************************** * Function Declarations * **************************/LRESULT CALLBACK WndProc (HWND hWnd, UINT message,WPARAM wParam, LPARAM lParam);void EnableOpenGL (HWND hWnd, HDC *hDC, HGLRC *hRC);void DisableOpenGL (HWND hWnd, HDC hDC, HGLRC hRC);int factorial(int x){    int rtn = 1;    for(int i = 1; i <= x; i++)        rtn *= i;    return rtn;}void SwapFloat(float &A, float &B){     float T = A;     A = B;     B = T;}void ShiftFloatArrayLeft(float *ARR, int Elements){     int i = 0;     for(i = 0;i<Elements-1; i++)     {         SwapFloat( ARR,ARR[i+1]);     } }void ShiftFloatArrayRight(float *ARR, int Elements){     int i = Elements-1;     for(i = Elements-1; i >= 1; --i)     {         SwapFloat( ARR,ARR[i-1]);     } }/************************** * WinMain * **************************/int WINAPI WinMain (HINSTANCE hInstance,                    HINSTANCE hPrevInstance,                    LPSTR lpCmdLine,                    int iCmdShow){    WNDCLASS wc;    HWND hWnd;    HDC hDC;    HGLRC hRC;            MSG msg;    BOOL bQuit = FALSE;    float theta = 0.0f;    /* register window class */    wc.style = CS_OWNDC;    wc.lpfnWndProc = WndProc;    wc.cbClsExtra = 0;    wc.cbWndExtra = 0;    wc.hInstance = hInstance;    wc.hIcon = LoadIcon (NULL, IDI_APPLICATION);    wc.hCursor = LoadCursor (NULL, IDC_ARROW);    wc.hbrBackground = (HBRUSH) GetStockObject (BLACK_BRUSH);    wc.lpszMenuName = NULL;    wc.lpszClassName = "GLSample";    RegisterClass (&wc);    /* create main window */    hWnd = CreateWindow (      "GLSample", "OpenGL Sample",       WS_CAPTION | WS_POPUPWINDOW | WS_VISIBLE,      0, 0, 640, 480,      NULL, NULL, hInstance, NULL);    /* enable OpenGL for the window */    EnableOpenGL (hWnd, &hDC, &hRC);        /************** local vars *******************/        const float Pi = 3.14159f;    const float X_step = 0.01f;    //Granularity in the X direction    const int numHarmonics = 9;    //Number of Fourier terms for Harmonics    const int maxPoints = Pi/X_step + 1; //Resoultion of the wave     float *Delay_B, *Delay_A, *FIN;    Delay_B = new float[maxPoints];    Delay_A = new float[maxPoints];    FIN = new float[maxPoints];    float HRM[numHarmonics];        float endpoint = Pi;        float MaxAmp = 0.3f;    float PluckPosition = Pi/2;     /*    for(int i = 0; i < numHarmonics; i++)    // Fourier     {        HRM =  MaxAmp * sin(PluckPosition*Pi*(i))/factorial(i);          //MaxAmp = MaxAmp - (sin(PluckPosition*Pi*(i))/(factorial(i)));            }        float x = 0;        for(int I = 0;  I<maxPoints; ++I, x+=X_step)    {            Delay_B = 0;            for(int j = 0; j < numHarmonics; j++)            {                  Delay_B += sin(x*j)*HRM[j];                  Delay_A = Delay_B;                              }                 }*/    float x = 0;    Delay_B[0] = 0;    Delay_A[0] = 0;    float Slope_A = (MaxAmp)/(PluckPosition),          Slope_B = (MaxAmp)/(PluckPosition-endpoint) ;    // ********* Fill The Delay Lines With initial string shape, pluck Triangle **************    for(int I = 1;  I < maxPoints - 1; ++I, x+=X_step)     {        if(x < PluckPosition)        {            Delay_B = x*Slope_A;            Delay_A = x*Slope_A;        }        else if(x > PluckPosition)        {            Delay_B = MaxAmp + (x-PluckPosition)*Slope_B;            Delay_A = MaxAmp + (x-PluckPosition)*Slope_B;        }        else        {            Delay_B = MaxAmp;            Delay_A = MaxAmp;        }    }    Delay_B[maxPoints - 1] = 0;    Delay_A[maxPoints - 1] = 0;                       /* program main loop */    while (!bQuit)    {        /* check for messages */        if (PeekMessage (&msg, NULL, 0, 0, PM_REMOVE))        {            /* handle or dispatch messages */            if (msg.message == WM_QUIT)            {                bQuit = TRUE;            }            else            {                TranslateMessage (&msg);                DispatchMessage (&msg);            }        }        else        {            //******************   Step WaveGuide     ***********************            for(int i = 0; i < 10; i++)            {                                           float Last_A  = Delay_A[maxPoints-1];   //Save Last Element from A                float First_B = Delay_B[0];             //Save First Element from B                ShiftFloatArrayRight(Delay_A,maxPoints);//Move waves right along delay line A                                   ShiftFloatArrayLeft(Delay_B,maxPoints); //Move waves left along delay line B                   Delay_B[maxPoints-1] = -1 * Last_A;     //Input Last Element from A into B and invert                             Delay_A[0]           = -1 * First_B;    //Input First Element from B into A and invert            }            for(int i = 0; i < maxPoints; i++)                    FIN = Delay_B+ Delay_A;                        /* OpenGL animation code*/            glClearColor (0.0f, 0.0f, 0.0f, 0.0f);            glClear (GL_COLOR_BUFFER_BIT);            glPushMatrix ();            glTranslatef(-.75, 0.0f, 0.0f);                 glScalef(0.5f,0.5f,0.5f);                  x = 0;            glColor3f(0,0,.7f); // Blue            glBegin (GL_LINE_STRIP);                 for(int i = 0; i < maxPoints; i++, x+=X_step)    // Draw Delay_B in Blue                {                               glVertex2f (x, Delay_B);                }            glEnd ();             x = 0;            glColor3f(0.7f,0,0); // Red            glBegin (GL_LINE_STRIP);                 for(int i = 0; i < maxPoints; i++, x+=X_step)   // Draw Delay_A in Red                {                                glVertex2f (x, Delay_A);                }            glEnd ();            glColor3f(1,1,1); // White            x = 0;            glBegin (GL_LINE_STRIP);                 for(int i = 0; i < maxPoints; i++, x+=X_step)  // Draw Combined A & B in white.                {                               glVertex2f (x, FIN);                }            glEnd ();                        glPopMatrix ();            SwapBuffers (hDC);        }    }    delete []Delay_B;    delete []Delay_A;    delete []FIN;    /* shutdown OpenGL */    DisableOpenGL (hWnd, hDC, hRC);    /* destroy the window explicitly */    DestroyWindow (hWnd);    return msg.wParam;}/******************** * Window Procedure * ********************/LRESULT CALLBACK WndProc (HWND hWnd, UINT message,                          WPARAM wParam, LPARAM lParam){    switch (message)    {    case WM_CREATE:        return 0;    case WM_CLOSE:        PostQuitMessage (0);        return 0;    case WM_DESTROY:        return 0;    case WM_KEYDOWN:        switch (wParam)        {        case VK_ESCAPE:            PostQuitMessage(0);            return 0;        }        return 0;    default:        return DefWindowProc (hWnd, message, wParam, lParam);    }}/******************* * Enable OpenGL * *******************/void EnableOpenGL (HWND hWnd, HDC *hDC, HGLRC *hRC){    PIXELFORMATDESCRIPTOR pfd;    int iFormat;    /* get the device context (DC) */    *hDC = GetDC (hWnd);    /* set the pixel format for the DC */    ZeroMemory (&pfd, sizeof (pfd));    pfd.nSize = sizeof (pfd);    pfd.nVersion = 1;    pfd.dwFlags = PFD_DRAW_TO_WINDOW |       PFD_SUPPORT_OPENGL | PFD_DOUBLEBUFFER;    pfd.iPixelType = PFD_TYPE_RGBA;    pfd.cColorBits = 24;    pfd.cDepthBits = 16;    pfd.iLayerType = PFD_MAIN_PLANE;    iFormat = ChoosePixelFormat (*hDC, &pfd);    SetPixelFormat (*hDC, iFormat, &pfd);    /* create and enable the render context (RC) */    *hRC = wglCreateContext( *hDC );    wglMakeCurrent( *hDC, *hRC );}/****************** * Disable OpenGL * ******************/void DisableOpenGL (HWND hWnd, HDC hDC, HGLRC hRC){    wglMakeCurrent (NULL, NULL);    wglDeleteContext (hRC);    ReleaseDC (hWnd, hDC);}

Winograd

380

August 29, 2005 07:21 AM

Quote:Original post by Grain
Ok I put together a very basic waveguide just to get accustomed to the concept.

Let me know if I got the basic concept. From what I read I under stand the Delay line to be a sequential container containing the wave’s amplitude at multiple points along the string. Is this correct.

Yes, this is correct. Delay line is just a buffer which is shifted by one unit at every "tick". Thus if you're sampling rate is 44100Hz then the unit is 1/44100s (as well as the "tick").

Quote:Original post by Grain
One thing I noticed right away was the shape the was tracing out wasn’t a sinusoid like here, but rather a parallelogram.

You're simulating a lossless string which does not lose energy at all. In reality when the wave travels, it loses energy due to internal "friction". This "friction" is frequency dependant, thus you would have to add a digital filter between every unit-delay in the delay line. Fortunately when we use linear time-invariant filters we can commute them (i.e group them from the delay line to the front of the delay line). The two delay lines can also be combined. This is done in Karplus-Strong algorithm. Because we combined the filters between unit-delays into one filter, the delay line no longer represents the true shape of the vibrating string. The result is that you have to do your measurement right after the filter. So you _can_ find the sound the string produces, but you _can't_ see the true shape of the string. Karplus-Strong Algorithm

In KP-Algorithm the plucking is done with white noise. This is just a choice made by the developers of the algorithm. White noise gives richer sound altough this does not realisticly correspond to plucking.

I found a good paper which explains well many things related to string synthesis with waveguides. Julius Smith explains shortly resistive losses in the string (the ones you were missing) and also frequency dependent damping. Read this and refer to JOS web-site for deeper information.

Quote:Original post by Grain
You can just copy and compile this, it’s a complete program.

*** Source Snippet Removed ***

Unfortunately I can't. I don't own a copy of Windows. Your program seems to implement the waveguide correctly. You could consider using a ring buffer to implement the delay line much more efficiently (altough this is not necessary).

Harmonics of a vibrating string.

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Harmonics of a vibrating string.

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Reticulating splines