How to remove 50 Hz frequency noise from a PCM file?

Question

Are you kind to tell me the odds or the means by which should I grab a clean signal with the $~50 \mathrm{Hz}$ hum removed using just programming?

I managed to make a program in Visual-C that composes a sinusoid and overlaps it over a signal in a file. The file is recorded using a microphone. At recording time using a PC I get the frequency hum of the AC, with the phase with which it overlaps the signal. I see that the noise's frequency is not influenced and remains constant over time.

Answer 1

You should just use a notch filter.

Below shows two sinusoids added together (top) and the result of filtering the lower frequency one (bottom).

The filter is just: $$ H(z) = \frac{1 - 2 \cos(2\pi \omega_{b}) z^{-1} + z^{-2}} {1 - 2 \alpha \cos(2\pi \omega_{b}) z^{-1} + \alpha^2z^{-2}} $$ where $\omega_b$ is the bad frequency and $\alpha<1.0$ is related to the notch width / depth.

---

R Code Below

#29936

fs <- 8000
fbad <- 50
phibad <- 2*pi*runif(1,0,1)
fgood <- 101.98340234
phigood <- 2*pi*runif(1,0,1)

T <- 1000
t <- seq(1,T)

x_noisy <-  sin(2*pi*fgood/fs*t + phigood) + sin(2*pi*fbad/fs*t + phibad) 

num <- c(1, -2*cos(2*pi*fbad/fs), 1)
alpha <- 0.99
den <- c(1, -2*alpha*cos(2*pi*fbad/fs), alpha*alpha)

x_filtered <- filter(num, den, x_noisy)

par(mfrow = c(2,1))
plot(t,x_noisy, type="l")

plot(t,x_filtered, type="l")

Answer 2

Please note that the AC frequency is precisely 50.0hz. Sampling rate is 44100. C-code below.

How to use it:

You can compile it using any C compiler, and your program.

How to integrate it:

• Load dwData with the first samples from your file 882 samples at a time (or any with silence on it).

• Call spamp() with arguments. It will modify pamp=… short global value.

• Call sphis() with arguments. It will modify phis=… float global value.

• Call rs() with same arguments. It will return dwDataOut[] shorts with cleaned signal.

• Save or use dwDataOut samples.

• Advance the next 882 samples and repeat the algorithm.

Precautions:

Neither of the procedures has verification algorithms implemented. You have to implement it yourself.

---

C Code Only Below

//global variables
short pdc=0;//N/A
short pamp=2643;//amplitude [short]
float phis=71.00000;//degrees
float pf=1.57080;//radians

////////////////////////////////////////////////////////////////////////
//
//      remove 50hz humming
//  modifies the signal to remove the effect of a 50hz hum.
//      uses global vars. (pamp, phis, pdc) which are modified by 
//      functions (spamp, sphis, sdc)
////////////////////////////////////////////////////////////////////////
void rs(short* dwData, short* dwDataOut, DWORD dwLength){
    pf=2*50*PI/180.0f*90/100.0f;
float omegaf=pf/(float)dwLength*4.0f;
DWORD dwi;
float fu;

    for(dwi=0;dwi<dwLength;dwi++){
        fu=omegaf*dwi+phis*PI/180.0f;
        while(fu>=2*PI)
            fu-=2*PI;
        dwDataOut[dwi]=short(dwData[dwi]*bEnableAntiPhase-1/sqrt(2.0f)*pamp*sin(fu)-pdc*0);
    };

    phis+=pf*180/PI*4.0f;
    while(phis>=360.0f)
        phis-=360;
}

////////////////////////////////////////////////////////////////////////
//
//  search-for-amplitude
//      searches for maximum amplitude in a vector
//      of values. Uses findmin & findmax functions
//
////////////////////////////////////////////////////////////////////////
void spamp(short* dwData, DWORD dwLength, short* iValueOut){
short pleft=0;
short pright=0;
DWORD dwm=0;
DWORD dwi=0;
TCHAR sztemp[MAX_PATH];

    if (dwData[1]){
        findmax(dwData, dwLength, &dwm);
        findmin(dwData, dwLength, &dwi);
        iValueOut[0]=(dwData[dwm]-dwData[dwi])/2.0f;
    }else MessageBox(cvar.get_hwnd(), “semnal nul”, “titlu spamp”, 0);
}

////////////////////////////////////////////////////////////////////////
//
//  search-for-phase
//      compares any signal with a sinusoid of 50hz (if dwLength==882)
//      and finds out the phase so that by composing them altogether
//      you get the cleaned signal.
////////////////////////////////////////////////////////////////////////
void sphis(short* dwData, DWORD dwLength, float* fValueOut){
short* dwDataOut=(short*)calloc(1000, sizeof(short*));
short* dwValueLeft=(short*)calloc(1000, sizeof(short*));
short* dwValueRight=(short*)calloc(1000, sizeof(short*));
short* dwValueC=(short*)calloc(1000, sizeof(short*));
short pleft=0;
short pright=0;
DWORD iLeft=0;
DWORD iRight=0;
DWORD dwm;
float phisOld=phis;

    for(dwm=0;dwm<360;dwm++){
        phis=dwm;
        rs(dwData, dwDataOut, dwLength);

        cpower(dwDataOut, dwLength, &pleft, &pright);

        dwValueLeft[dwm]=pleft;
        dwValueRight[dwm]=pright;
        dwValueC[dwm]=(short)dwm;
    };

    phis=phisOld;

    findmin(dwValueLeft, dwm, &iLeft);
    findmax(dwValueRight, dwm, &iRight);

    fValueOut[0]=(dwValueC[iLeft]+dwValueC[iRight])/2.0f

    free(dwDataOut);
    free(dwValueLeft);
    free(dwValueRight);
    free(dwValueC);
}

////////////////////////////////////////////////////////////////////////
//
//  integrate by portions the amplitudes in a buffer (positive plus negative)
//
////////////////////////////////////////////////////////////////////////
void cpower(short* dwData, DWORD dwLength, short* pleft, short* pright){
DWORD dwi;
DWORD dwl=0;
DWORD dwr=0;
long cpl=0;
long cpr=0;

    for (dwi=0;dwi<dwLength;dwi++)
        if (dwData[dwi]>=0){
            cpl+=dwData[dwi];
            dwl++;
        }else{
            cpr+=(dwData[dwi]);
            dwr++;
        };

    pleft[0]=(short)(cpl/(float)(dwl+1*0));
    pright[0]=(short)(cpr/(float)(dwr+0*1));
}

void findmax(short* dwData, DWORD dwLength, DWORD* dwIndex){
DWORD dwi=0;
short dwValue=0-32768;//plec din minim

    if (dwLength==0){
        MessageBox(cvar.get_hwnd(), "dwLength==0", "titlu findmax", 0);
        return;
    };

    for(dwi=0;dwi<dwLength;dwi++)
        if (dwData[dwi]>=dwValue){
            dwValue=dwData[dwi];
            dwIndex[0]=dwi;
        }else{
            dwValue=dwData[dwIndex[0]];
        };

}
void findmin(short* dwData, DWORD dwLength, DWORD* dwIndex){
DWORD dwi;
short dwValue=32767;

    for(dwi=0;dwi<dwLength;dwi++)
        if (dwData[dwi]<dwValue){
            dwValue=dwData[dwi];
            dwIndex[0]=dwi;
        }else{
            dwValue=dwData[dwIndex[0]];
        };
}