Eliminare 50hz din microfon cu MVC-6

[Vizitator win-paul]

Buna dimineata,

 

Am realizat un program pe comp. in MVC-6 (microsoft visual c-6) care compune o sinusoida de frecventa pf, amplitudine pamp, si faza fphi, il compune cu un semnal WAV. Problema care ma deranjeaza pe mine provine dela faptul ca inregistrarea WAV se va face cu microfonul. La inregistrare apare purtatoarea frecventei de retea, cu faza cu care acesta se suprapune peste semnal, si constat ca purtatoarea nu este mentinuta constanta in timp – pe intervale mai lungi, la inregistrarea cu microfonul compului, pe langa faza care trebuie mereu s-o caut. Eu vreau sa obtin semnalul curat din microfon, cu zgomote mici acolo, dar cu magnitudine nu atat de mare la 50hz.

 

Care ar fi sansele/metoda ca sa pot obtine un semnal curat cu purtatoarea de 50hz eliminata utilizand doar programarea software? Mai exista vreo sansa cu un c.i. extern (565).

 

Stiu ca ar trebui un filtru opreste banda, dar pentru o singura frecventa poate exista ceva mai usor (o legatua inversa de reactie care intra in preamp-a placii de sunet) insa in elt. nu sunt cel mai priceput. Folosesc o placa de sunet onboard, 44k1, 16bit.

 

fisierul WAV arata cam asa:

Editat Aprilie 4, 2016 de win-paul

[Vizitator win-paul]

Vreau sa stiti ca la inregistrare s-a folosit o placa de sunet Realtek onboard, si ca am constatat ca pentru conversie se consuma o cantitate mica de timp, pierzandu-se astfel sincronizarea de faza cu semnalul de 50hz pur din retea, pe care o pot calcula. Din aceasta cauza programul meu nu poate face compunerea aceasta in mod realtime.

 

Nu este nimeni amabil si care sa poate face o schema cu c.i. β E 565?

[cirip]

Salut,Incerc sa te ajut. Posteaza te rog wav-ul. Nu te astepta la rezultate de azi pe maine, ca nu e chiar asa de simplu.Da, se poate face soft si in timp real. Nu ai nevoie de niciun fel de hardware sau circuit suplimentar.Sunt mai multe tehnici posibile de filtrare, dar mai intai ar trebui analizat spectrul. Uitandu-ma la poza, suspectez ca nu ai numai 50Hz, ci si ceva armonice. Asa pe viteza imi vin in minte 2 metode: un notch cu biquad sau o filtrare adaptiva.Ar mai fi cateva chestii:Nomenclatura - Aia de la retea nu e "purtatoare", ca nu poarta nimic. Ala e brum de 50Hz. Purtatoare e aia modulata, de frecventa inalta, care "poarta" muzica sau vorbe pana la difuzorul radioului.Bube - Microfonul ala are ceva bube. Nu am intalnit niciodata un microfon intern care sa capteze asa un brum. Ceva nu se leaga. E chiar microfonul intern sau e unul extern cuplat la mufa? Daca e extern, o fi buna ecranarea?Posteaza wav-ul si vedem.

Editat Aprilie 6, 2016 de cirip

[Vizitator win-paul]

Am inteles diferenta dintre “brum” si “purtatoare”; multumesc mult ca mi-ai atras atentia cum se diferentiaza.

 

Microfonul este unul capacitiv extern de calculator cu picior stativ.

 

Nu mai este necesar niciun circuit cu c.i. fiindca am schimbat calculatorul si vad ca zgomotul de 50hz a disparut complet. ramane insa intrebarea urmatoare: cum se implementeaza filtrarea in software, care sa fie cel mai putin sensibila la variatiile mici de faza (intre frecventa filtrata, si un interval centrat pe frecventa cu pricina)?

Daca are cineva o astfel de implementare rog sa posteze aici.

Multumesc mult cirip pt raspuns!

 

Nu cred ca este posibil sa se faca prin software filtrarea frecventei la modul cum ma gandeam eu.

[cirip]

Nu cred ca este posibil sa se faca prin software filtrarea frecventei la modul cum ma gandeam eu.

Posteaza te rog wav-ul cu brum. Chiar ar fi un exercitiu interesant.

[Vizitator win-paul]

Fisierul este mai mare, am postat numai o parte va rog nu va suparati pe mine.

numere-15.wav

Editat Aprilie 7, 2016 de win-paul

[cirip]

Salut Paul,M-am jucat un pic cu fisierul tau si, intr-adevar, filtrarea este foarte dificila. Motivul nu este acela ca nu se poate prin soft, ci pentru ca zgomotul este de banda foarte larga. Contine 50Hz, dar sunt o gramada de alte componente care se suprapun peste spectrul vocal. Chestia asta face filtrarea imposibila. Daca erau numai cei 50Hz era destul de usor. Probabil ca de asta ai avut si tu dificultati.Cateva poze ...Spectrul original (albastru) si cel filtrat (portocaliu). Se vede ca pe langa cei 50Hz, mai exista si zgomot de banda destul de larga.Semnalul in timp, original si filtratSemnalul in timp marit in zona care are numai zgomot, original si filtrat. Se vede ca brumul de 50Hz a fost mult atenuat, dar raman componentele din banda vocala.Si wavul filtrat ...numere_15_filtrat.zipChiar daca nu a iesit mare valoare, a fost un exercitiu interesant , iar problema ta e oricum rezolvata de ieri.Spor in continuare!Cirip

[Vizitator win-paul]

salut Cirip,

 

 

Filtrarea ta a iesit foarte bine. Eu reusisem sa fac ce mi-am propus dar vad ca partea de

program iarasi nu merge (cu operatii de scadere). Cu ocazia asta va rog, as vrea sa verific

programul meu pt care as avea nevoie de putina indrumare:

-se incarca fisierul “50hz.wav” cu faza zgomotului de fond

-se calculeaza parametrii: amplitudine zgomot, faza zgomot, frecventa de baza, si

dc-offset. Pentru asta se fol cmd “sepsilon0”.

-cu cmd “Save values (pf, phis, pamp, pdc)” se salveaza parametrii masurati

-se incarca alt fisier, cel cu spectrul vocal (dar care pastreaza faza phis)

-se aplica functia “Remove signal rc0” cu parametrii cititi din fisierul cu parametrii.

-automat se obtine fisierul “50hz1.00.WAV” calculate prin metoda mea :p (scadere).

 

Ca sa fie simplu am generat un fisier de 50hz, cu faza diferita de 0 pe care l-am

prelucrat ca mai sus, si acum postez rezultatul. Acum esti bun sa-mi spui o metoda ca sa

caut/calculez variabila 7.1189f ca sa fie cat mai exact. Adica rezultatul meu sa fie

ca “numere_15_filtrat”? As dori si cateva cuvinte despre procedura urmatoare...

 

Iata sectiunea de program cu probleme:

 

//variabile globale

pdc=0;

pamp=10454;//amplitudinea zgomotului generat

phis=2.41667;//faza zgomotului generat. [rad]

pf=0.87266;//frecventa retea [rad], se calc. din 44100/882.0f*PI/180.0f

 

 

////////////////////////////////////////////////////////////////////////

//

// remove 50hz carrier

// modifica semnalul pentru a scoate efectul unui semnal de 50hz.

// foloseste var. globale (pamp, phis, pdc) pe care le gasesc

// functiile (spamp, sphis, sdc)

////////////////////////////////////////////////////////////////////////

void rs(short* dwData, short* dwDataOut, DWORD dwLength){

float omegaf=pf/(float)dwLength*7.1189f;

DWORD dwi;

float fu;

 

// MessageBox(cvar.get_hwnd(), "inainte de bucla", "titlu rs", 0);

 

for (dwi=0;dwi<dwLength;dwi++){

fu=dwi*omegaf+phis;

while(fu>=2*PI)

fu-=2*PI;

 

dwDataOut[dwi]=short(dwData[dwi]-pbEnableAntiPhase*pamp*sqrt(2.0f)*sin(fu)-pdc*0);

};

 

//unda stationara

phis+=pf*7.1189f;

//reducere la cerc;;

while(phis>=2*PI)

phis-=2*PI;

}

 

 

 

Daca esti amabil imi poti trimite prin privat si procedura matematica prin care sa trasez eu

fig.1 ca sa nu deschid un alt topic te rog? MS :-)

50hz1.00.WAV

Editat Septembrie 21, 2017 de win-paul

[cirip]

Salut win-paul,

 

"Filtrarea ta a iesit foarte bine"

Multumesc de apreciere.

 

"-se calculeaza parametrii: amplitudine zgomot, faza zgomot, frecventa de baza, si dc-offset. Pentru asta se fol cmd “sepsilon0”.
Sunt curios cum calculezi parametrii zgomotului de fond. Ala e aleator, ca e zgomot, deci cum ii poti calcula faza, amplitudinea si alte marimi instantanee care-l caracterizeaza?
Habar n-am ce e "sepsilon0".

"o metoda ca sa caut/calculez variabila 7.1189f ca sa fie cat mai exact."
Cine e variabila asta, ce vrea ea si la ce se foloseste? Cum contribuie ea la bunastarea si progresul programului tau? Vad pe-acolo //unda stationara. //reducere la cerc, dar habar n-am despre ce e vorba.

 

"Iata sectiunea de program cu probleme:"
Nu o pot urmari. Descrie intai cu ecuatii algebrice ce algoritm matematic vrei sa implementezi si mai ales verifica daca intr-adevar face ce vrei. Procesarea de semnale e un pic diferita de programarea logica.

"trimite prin privat si procedura matematica prin care sa trasez eu fig.1"
Pai nu e niciun secret. Este Transformata Fourier.

Editat Septembrie 22, 2017 de cirip

[Vizitator win-paul]

multumesc ca mai sunteti aici sa cititi, chiar dupa atata timp de inactivitate din partea mea.

 

In program folosesc datele din fisiere wav deci faza zgomotului este deja stabilita in ele: trebuie doar sa stii s-o cauti.

In fisierul "numere-15.wav" nu mai stiu cat este precis, dar in 50hz1.00.WAV este phis = -40grade adica phis=320grade.

 

Eu nu am putut sa fac sa mearga procedura de capturare in timp-real, scopul programului fiind si de procesare avand un necesar de timp

de procesare putin mai indelungat decat timpul capturat. Astfel am zis ca este mai sigur sa lucrez cu fisiere: faza amplitudinii ramane acolo

imprimata, ca si restul parametrilor zgomotului. Desigur asta inseamna sa pot izola niste sample-uri din care sa-i pot deriva. Eu zic ca-mi ajung

882samples de liniste esantionati la 44100hz.

 

Cam multa vorba...

Dupa ce reusesc sa fac programul sa sa detecteze macar faza, aplitudinea, si dc-offset-ul unui zgomotului de fond,

o sa postez aici codul poate va ajuta pe cineva.

 

la revedere, Paul

[Vizitator win-paul]

cirip, on 21 sept 2017 said:

"-se calculeaza parametrii: amplitudine zgomot, faza zgomot, frecventa de baza, si dc-offset. Pentru asta se fol cmd “sepsilon0”.
Sunt curios cum calculezi parametrii zgomotului de fond. Ala e aleator, ca e zgomot, decicum ii poti calcula faza, amplitudinea si alte marimi instantanee care-l caracterizeaza?
Habar n-am ce e "sepsilon0".

*sepsilon0 este procedura din programul meu care face tocmai asta. Din pacate se pare ca nu o face foarte bine, dar algoritmul ar fi urmatorul:

*algoritm ptr eliminare 50hz de faza phi din microfon:

1.izolez o sinusoida (2PI) din zgomot, adica ma opresc sa analizez 882sample-uri din fisier.

2.Generez o sinusoida de 50hz in tot atatea sample-uri, adica dupa formula

dwDatasin(dwi) = pamp sin(omega f * dwi+ phi) = pamp sin(const * dwi + phi).

const pt ca frecv. retea e const. Am constatat ca ac. const. este

const = 44100/882*PI/180 = 0.87266.

3.Reactie de compunere cu restul sinusoidei (acelasi phi dupa 882sample-uri, nu?) va da
dwDataOut[dwi] = dwData[dwi]-dwDatasin[dwi] = 0
pt orice dwi := 0..882. 882 este durata unei sinusoide. Se obtine un vector dwDataOut nul (avand efectul semnalului de zgomot scos).
Avantajul folosirii ac. algoritm este ca odata stabilit defazajul phi si amplitudinea pamp, vom obtine semnalul curatat de zgomot, dar dupa cele 882frames in care se face analiza pt phi si pamp. Mai merita sa postez codul pt searchpamp() si searchphis()? Eu zic nu deocamdata ca nu functioneaza corect. Si aici revin la titlul topicului care asta cere.



cirip, on 21 sept 2017, 4:29 said:
"o metoda ca sa caut/calculez variabila 7.1189f ca sa fie cat mai exact."
Cine e variabila asta, ce vrea ea si la ce se foloseste? Cum contribuie ea la bunastarea si progresul programului tau? Vad pe-acolo //unda stationara. //reducere la cerc, dar habar n-am despre ce e vorba.
*a)este vorba despre o reducerea la cerc a unui numar (parca asa se spune). De exemplu sin (3 PI) se reduce la cerc ca sin (PI). Problema deriva din faptul ca MVC6 da eroare cand inghite sin(3 PI) si de aceea, repet in cazul lui MVC-6 trebuie ca unghiul-argument al functiilor trigonometrice simple sa fie "redus la cerc".

*b)7.1189f trebuia sa fie (2*PI), am scris eu gresit - vina mea. Este corecta formula de generat semnalul

pamp*sin(omega f t + phi)

, nu?

Problema este cum trec hertii din retea in formula cu sinus: ca unghi sau ca radiani? De asta si dureaza atat de mult pt mine sa finalizez secventa sepsilon0.

Acum...

omega f = 0.87266 / 2PI * dwi, unde 2PI eu il notasem gresit cu 0.7198f, dar chiar asa se poate obtine formula (*1)?

 

 

cirip, on 21 sept 2017, 4:29 said:

"trimite prin privat si procedura matematica prin care sa trasez eu fig.1"

Pai nu e niciun secret. Este Transformata Fourier.

*nu gaseam pe internet implementare ptr procedura matematica Fourier folosind MVC-6, si date de intrare short. N-am zis ca e secret, doar ca nu stiu... doar ca mi-ar prinde bine in program. MS

 

cam astea ar fi.

Editat Septembrie 23, 2017 de win-paul

[Vizitator win-paul]

*Algoritm pentru gasirea lui pamp (amplitudine zgomot de 50hz).

1.se stie pentru retea 2PIf = 0.82667 / 2 PI (si nu imartit la chiar 7.1189 ca mai inainte)

2.se calculeaza pentru bucla 2PI, marimea putere din semisinusoida pozitiva pleft dupa formula

if (dwData[dwi]>=0)

pleft+=dwData[dwi];

3.deasemenea pentru semisinusoida negativa, pright

else

pright+=dwData[dwi];

4.se numara cate sample-uri sau acumulat in semisinusoida negativa si in cea pozitiva

dwl++;//nr de sample-uri pozitive

dwr++;//nr de sample-uri negative

5.puterea totala in sinusoida este suma modulelor celor doua semisume (pright este negativa deci se ia cu semn schimbat)

pamp = (pleft[0]-pright[0])*sqrt(2).

Asa gasesc amplitudinea.

 

 

*Algoritm pentru gasirea lui phis (unghiul de faza al zgomotului de 50hz).

1.se stie ca datele contin un zgomot (dwData) de 50hz, dar poate si altele. Pentru a gasi pe acela de 50hz, nu trebuie decat sa compar mostrele care compun o bucla 2PI de 50hz cu o sinusoida curata de 50hz avand faza care se schimba (cred ca pot s-o numesc faza alunecatoare).

2.scriu formula zgomotului, in cazul nostru de 50hz.

float omegaf=0.87266f;

dwDatasin[dwi]=pamp sin(omegaf * dwi + phi);

3.Vectorul dwData il adun cu dwDatasin, pe elemente. Adica

dwDataOut[0]=dwData[0]+dwDatasin[0];

dwDataOut[1]=dwData[1]+dwDatasin[1];

dwDataOut[2]=dwData[2]+dwDatasin[2];

samd pana la

dwDataOut[360]=dwData[360]+dwDatasin[360];

4.caut minimul din matricea dwDataOut, si gasesc de exemplu ca el este elementul din

dwDataOut[40];

5Asta inseamna ca 40 reprezinta gradele defazajului sau

phis=40+180;//

6.Daca suplimentar vrem sa anulam brumul de 50hz, scriem

fu = dwi omegaf + 40;

short(dwData[dwi]+pbEnableAntiPhase*pamp*sqrt(2.0f)*sin(fu)-pdc*0);

, totul pentru dwi dela 0 la 882.

7.reimprostatez pe dwData[dwi] si tratam secvential in blocuri de cate 882 pana terminam tot fisierul.

 

Acesta e algoritmul pt gasirea defazajului. Am scris 0.87266f deci e pentru 50hz.

 

Sunt bune algoritmurile acestea? Mie mi se par in regula si le-am implementat. Chestia este ca nu merge si imi rezulta fisiere defazate care parca ar fi modulate in frecventa.

[cirip]

Salut win-paul,

 

Este destul de greu sa urmaresc ce scrii, fie din cauza numeroaselor erori, fie din cauza abordarii principiale incorecte. Exemple:

 

"se stie pentru retea 2PIf = 0.82667 / 2 PI (si nu imartit la chiar 7.1189 ca mai inainte)"

2*PI*f = 2*3.141592*50 = 314.1592. Nu stiu de unde apar sau ce reprezinta alte valori.

 

"se calculeaza pentru bucla 2PI, marimea putere"

Ca sa discutam de putere, trebuie sa apara o ridicare la patrat. Altfel, ceea ce faci seamana cu o integrare a sinusului pe portiuni. E necesara si asta daca vrei sa calculezi valoarea medie (care la sinus e oricum 0), dar procedura ta este incompleta.

Mai mult, rezultatul procedurii tale este dependent de rata de esantionare, ceea ce nu este corect. Puterea unei sinusoide de o anumita amplitudine debitata intr-o rezistenta de 1 Ohm trebuie sa fie aceeasi indiferent de rata de esantionare folosita, presupunand ca satisfaci teorema esantionarii (Nyquist-Shannon).

 

"puterea totala in sinusoida[...] Asa gasesc amplitudinea"

Pai cum trecusi de la putere la amplitudine? Si ce te faci, ca semnalul tau nu contine numai sinusoida de 50Hz. Uita-te in in postul din Aprilie 2016 in care am postat spectrul si semnalul in timp. Acolo ai suprapuse o multime de componente. Cand faci suma aia pe o perioada, incluzi multe altele.

 

"pamp*sin(omega f t + phi)"

Cum adica "omega f" ? Pai daca inmultesti omega cu f, apare frecventa la patrat. Asta e intentia?

 

"Mie mi se par in regula si le-am implementat. Chestia este ca nu merge"

Prin urmare nu sunt in regula.

 

Deocamdata ma opresc aici. Hai sa o luam pe rand, ptr ca sunt prea multe de discutat. Parca as incerca sa beau apa de la furtunul de pompieri

 

Cirip

Editat Septembrie 24, 2017 de cirip

[Vizitator win-paul]

Unde eu ziceam "puterea sinusoidei" era mai corect "integrare pe portiuni". Da!

"se stie pentru retea 2PIf = 0.82667 / 2 PI (si nu imartit la chiar 7.1189 ca mai inainte)"

2*PI*f = 2*3.141592*50 = 314.1592. Nu stiu de unde apar sau ce reprezinta alte valori.

 

"se calculeaza pentru bucla 2PI, marimea putere"

"puterea totala in sinusoida[...] Asa gasesc amplitudinea""pamp*sin(omega f t + phi)"

Cum adica "omega f" ? Pai daca inmultesti omega cu f, apare frecventa la patrat. Asta e intentia?

 

 

Este corect algoritmul pt integ. pe portiuni a amplitudinii zgomotului?

 

"se stie pentru retea 2PIf = 0.82667 / 2 PI (si nu imartit la chiar 7.1189 ca mai inainte)"

2*PI*f = 2*3.141592*50 = 314.1592. Nu stiu de unde apar sau ce reprezinta alte valori.

trebuie sa lucrez in radiani. Dar formula nu este 2*PI*f*t+phi?

 

Deasemenea, in algoritmul pt aflarea fazei, pasii se repeta _doar_ pana la

dwDataOut[359]=dwData[359]+dwDatasin[359];

 

"pamp*sin(omega_f * t + phi)" (*n.r. transcriere modificata)

Cum adica "omega f" ? Pai daca inmultesti omega cu f, apare frecventa la patrat. Asta e intentia?

 

Nu asta era intentia. omegaf este un param care include parametrul de frecventa 50 si se exprima in radiani. omegafXPI/180Xt=0.87266xt.

Deasemenea el include si parametrul frecventa de esantionare f=44100.

 

Nu tin neaparat sa lucrez in putere pentru ca necesita si formule pt decibeli si presiune acustica, pe care nu stiu acum sa le aplic in program. Formula pentru putere acustica este P=10 log(p2/p1) [dB], unde p1 este o presiune acustica de referinta.

Trebuie folosita si aceasta?

 

"amplitudine a unei sinusoide" per total este zero insa caut (integrez pe portiuni) o semisinusoida de 50hz, pe care o inmultesc cu 2.

 

 

La revedere, MS cirip

Paul

Editat Septembrie 24, 2017 de win-paul

[Vizitator win-paul]

I-am dat de cap , aveam o greseala in procedura de gasire a maximului dintr-un vector de valori – findmin() si findmax().

 

Pt cirip:

Pt urmatoarea val. a frecventei din retea se genereaza o sinusoida de 50hz

pf=2*50*PI/180.0f*90/100.0f;

float omegaf=pf/(float)dwLength*4.0f;

dwLength=882;

si sin (omegaf x dwi + phi)

cu dwi pe post de variabila timp 0 < dwi < dwLength

semnalul se filtreaza ca in fisier.

 

Problema a aparut ca urmare a generarii sinusoidei, de ex de 50hz. Frecventa ei trebuie adusa in radiani (deci sin (2 PI 50) apoi redus la cerc) si, zic eu, se transforma gradele centizecimale in cele 90 din sfertul cercului trigonometric. Mai mult nu stiu. Un topic interesant ar fi "cum sa generez o sinusoida de frecventa data".

 

Prefer filtrarea aceasta pt ca nu apar armonici ale lui 50hz.

Imi pare rau nu pot uploada fisierul, nu merge flash-ul o sa rog webmasterul poate ma ajuta

Paul

numere-15-2.orig+filtrat.00.WAV