1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Welcome to SE.SP. You apprarently already have performed soome work on your problem. You might have a look at "50 Hz notch filters", or "50 Hz band-stop filters", perhaps with hints here electronics.stackexchange.com/questions/36005/… $\endgroup$ – Laurent Duval Apr 6 '16 at 6:11
  • $\begingroup$ If this is 50 Hz noise induced by AC power line (aka "hum"), removing 50 Hz is typically not sufficient. While the fundamental frequency of the hum is 50 Hz (in Europe) the most energy tends to be in the harmonics (100 Hz, 150 Hz, etc.) $\endgroup$ – Hilmar Apr 6 '16 at 18:24
  • $\begingroup$ How to add pins input and output ? $\endgroup$ – ted Feb 8 '18 at 15:20
2
$\begingroup$

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.

enter image description here


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")
$\endgroup$
0
$\begingroup$

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]];
        };
}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.