Ik zoek een algoritme om een ruisbandje te filteren.

[Bert F]

Ik heb wat gelezen op je site en zie dat je toch al wat watertjes doorzwomen heb daarom denk ik ook dat volgende ook wel moet lukken:

Code: Selecteer alles

fs=20000;

 % sampling rate

F=[3000 4000 6000 8000];  % band limits

A=[0 1 0];

% band type: 0='stop', 1='pass'

dev=[0.0001 10^(0.1/20)-1 0.0001]; % ripple/attenuation spec

[M,Wn,beta,typ]= kaiserord(F,A,dev,fs);  % window parameters

b=fir1(M,Wn,typ,kaiser(M+1,beta),'noscale'); % filter design

dit is een stukje matlab code waarmee je snel een filter operaties kan doorvoeren.

*)specificeer je sample frequentie

*)specificeer wat je limits zijn van je bandpass filter

*)dit is een bandpass filter een lowpass is [1 0 0] een high pass [0 0 1]

*)specifeer je atenuaties in je gabarit

*)de functie kaiserord maakt met behulp van deze eigenschappen een aantal andere parameters die je dan op zijn beurt weer invoerd bij fir1

*)nu kan je een y = filter(b,a,X) waarbi X je getallen zijn je b net bepaalt hebt en a gewoon op 1 kan stellen

y geeft nu de gefilterde getallen terug.

Matlab is een duur pakket en enkel voor dit moet je dat niet aankopen maar octave is een gelijkaardig pakket, ik heb voor jouw nagekeken en volgens mij werkt bovenstaande stukje code ook daarin zie hier voor een referenties van de functies in oactav:

http://octave.sourceforge.net/functions_by_alpha.php

rest nu nog je getalen in te lezen deze in je vector x te steken een filter te kiezen door de getaltjes te veranderen en het resultaat terug in je tekst file op te slaan.

Het spijt me, maar ik heb gisteren net mijn computer geformateerd en mijn matlab cd ligt nog op mijn kot anders kon ik je snel wel eens een voorbeeldje geven post mss je getallen alvast in een gewoon tekst bestand.

Groeten.

[Wimapon]

317070

Hoi, heel erg bedankt voor dit werk.

Ik snap niet helemaal wat je bedoelt.

Het eerste deel snap ik y(n) =( 1 * x(n-20)) + ( 0 * x(n-19)) etc etc....

maar na de lege regel .... moet ik ineens een antwoord er bij optellen dat ik nog niet berekend heb..... namelijk y(n-20) en y(n-19) etc...

+ (-0.001782875 * y(n-20)) + ( 0.0387309183 * y(n-19)) etc etc etc...

Ik hoop dat ik je uitleg verkeerd begrijp..

[Wimapon]

Hoi Bert,

dankjewel

Mijn probleem is dat ik alles zelf wil programmeren... want dan heb ik alles in eigen hand.

tenzij ik vanuit een Quick Basic programma iets van Matlab of Octave kan aanroepen...

Het punt is dat ik een hele serie programma's achter elkaar laat werken om tot een sterrenkaart te komen ..

Ik kan dus niet 100 000 maal matlab opstarten....( grin)

Ik weet het, ik kies niet de gemakkelijkste weg.

Trouwens het lukt me ondanks je bericht niet om een bestand in deze post te krijgen..

ik kan gemakkelijk een rij getallen of een .jpg op mijn eigen site zetten.

Als ik nu vanuit dit forum een link kan plaatsen zodat het plaatje hier te zien is... zou ik heeeel gelukkig zijn.

Ik snap wel dat het moet kunnen, maar ik kan het gewoon niet vinden....

voorlopig kan ik natuurlijk een meting op mijn site zetten, die zou dan het url krijgen: http://home.kpn.nl/apon001/vector1.dat bijvoorbeeld..

Ik denk dat je die dan zo binnen kunt halen.. het zijn dan 1024 * 8 getallen in ascii format

als je het wil, zal ik het direct doen

[317070]

317070

Hoi, heel erg bedankt voor dit werk.

Ik snap niet helemaal wat je bedoelt.

Het eerste deel snap ik y(n) =( 1 * x(n-20)) + ( 0 * x(n-19)) etc etc....

maar na de lege regel .... moet ik ineens een antwoord er bij optellen dat ik nog niet berekend heb..... namelijk y(n-20) en y(n-19) etc...

+ (-0.001782875 * y(n-20)) + ( 0.0387309183 * y(n-19)) etc etc etc...

Ik hoop dat ik je uitleg verkeerd begrijp..

Wel, in het begin bestaat y uit allemaal 0'en. Daarna bereken je y[20], hiervoor heb je y[0] tot y[19] nodig. Daarna bereken je y[21], hiervoor heb je y[1] tot y[20] nodig, (je ziet: y[20] heb je al berekend). Daarna bereken je y[22], hiervoor heb je y[2] tot y[21] nodig, enz...

[Wimapon]

hoi 317070

Dus ik moet gewoon al die 41 regels bij elkaar optellen , en zo het hele array doorlopen..

Nou, dat moet lukken.

Ik ga het proberen.

Hartstikke bedankt

[Wimapon]

317070

HOi, ik het de besseltruuk even geprogrammeerd.

Het spectrum ziet er geweldig uit, maar het uitgangs signaal is niet goed meer.

De eerste 100 samples laten precies de gewenste frequentie zien, maar de rest van de 8000 samples staan zo goed als op 0.......

Ik wil heel graag het plaatje laten zien. maar hoe moet ik dat hier nu doen!!!!!

ik denk dat je het plaatje hier kunt zien--- > http://home.kpn.nl/apon001/bessel.jpg

Links boven is het spectrum van het ingangssignaal ( is 1 meting van 1 antenne)

links onder staat het signaal zelf... de eerste 300 samples

Rechtsboven staat het spectrum van het uitgangssignaal

rechtsonder staat dit signaal zelf..... ( ik het het 2.4 volt opgetild om in het de grafiek te krijgen)

Je ziet alleen links een paar mooie sinussen.... de rest van de samples.. ( tot 8000) staan op 0

[Bert F]

Ik denk dat je die dan zo binnen kunt halen.. het zijn dan 1024 * 8 getallen in ascii format

Dat zou al helpen. uiteindelijk moet het wel mogelijk zijn om zo'n klein progje aan te roepen vanuit jouw programma, dus uiteindelijk te komen tot één programma.

e ziet alleen links een paar mooie sinussen.... de rest van de samples.. ( tot 8000) staan op 0

Kan je, je code daarvan er even bijzetten?

[Wimapon]

Hoi Bert F,

Ik heb een meting van 8 * 1024 getallen op mijn site gezet.

Deze meting is ongeveer in 0.4 seconde gedaan.

http://home.kpn.nl/apon001/verctor1.txt"

[Wimapon]

mocht het downloaden niet lukken dan staat hij ook in:

http://home.kpn.nl/apon001/vector1.htm

even de html dingen eraf slopen en dan heb je 8 * 1024 getallen

[Wimapon]

de basic code waarmee ik filter staat hieronder:

/code

' dit prog probeert een bessel band filter te berekenen

' input = kanaala.dat

' output = kanaalb.dat





CLS

DEFINT I-J



aantal = 1024 * 8



ixax = 9000

DIM x(ixax) AS DOUBLE

DIM y(ixax) AS DOUBLE

' leest eerst kanaala in in het eerste array



OPEN "i", #1, "kanaala.dat"

FOR i = 1 TO 1024 * 8

INPUT #1, x(i)

NEXT i

CLOSE #1



' doe nu een bessel filtering volgens 317070 van het forum

FOR n = 21 TO aantal - 21

y(n) = (1 * x(n - 20))

y(n) = y(n) + (0 * x(n - 19))

y(n) = y(n) + (-10 * x(n - 18))

y(n) = y(n) + (0 * x(n - 17))

y(n) = y(n) + (45 * x(n - 16))

y(n) = y(n) + (0 * x(n - 15))

y(n) = y(n) + (-120 * x(n - 14))

y(n) = y(n) + (0 * x(n - 13))

y(n) = y(n) + (210 * x(n - 12))

y(n) = y(n) + (0 * x(n - 11))

y(n) = y(n) + (-252 * x(n - 10))

y(n) = y(n) + (0 * x(n - 9))

y(n) = y(n) + (210 * x(n - 8))

y(n) = y(n) + (0 * x(n - 7))

y(n) = y(n) + (-120 * x(n - 6))

y(n) = y(n) + (0 * x(n - 5))

y(n) = y(n) + (45 * x(n - 4))

y(n) = y(n) + (0 * x(n - 3))

y(n) = y(n) + (-10 * x(n - 2))

y(n) = y(n) + (0 * x(n - 1))

y(n) = y(n) + (1 * x(n - 0))

y(n) = y(n) + (-.0017820875# * y(n - 20))

y(n) = y(n) + (.0387309183# * y(n - 19))

y(n) = y(n) + (-.4111305294# * y(n - 18))

y(n) = y(n) + (2.8308581223# * y(n - 17))

y(n) = y(n) + (-14.1680177848# * y(n - 16))

y(n) = y(n) + (54.7519171722# * y(n - 15))

y(n) = y(n) + (-169.4387712114# * y(n - 14))

y(n) = y(n) + (429.8313691642# * y(n - 13))

y(n) = y(n) + (-907.5633411356# * y(n - 12))

y(n) = y(n) + (1610.3756932883# * y(n - 11))

y(n) = y(n) + (-2414.0692545962# * y(n - 10))

y(n) = y(n) + (3062.2622671977# * y(n - 9))

y(n) = y(n) + (-3280.7595198059# * y(n - 8))

y(n) = y(n) + (2951.7958480361# * y(n - 7))

y(n) = y(n) + (-2207.9535492666# * y(n - 6))

y(n) = y(n) + (1351.3740145742# * y(n - 5))

y(n) = y(n) + (-660.5471327188# * y(n - 4))

y(n) = y(n) + (248.3330369343# * y(n - 3))

y(n) = y(n) + (-67.49323853# * y(n - 2))

y(n) = y(n) + (11.812001138# * y(n - 1))

NEXT n

' schrijf het resultaat nu naar de output-file

OPEN "o", #2, "kanaalb.dat"

FOR i = 1 TO 1024 * 8

PRINT #2, y(i) / 500000 'anders gaat ie te ver offscale

NEXT i

CLOSE #2



END

/endcode

[Bert F]

ik heb je data verwerkt met volgende gegevens:

Code: Selecteer alles

tijd_meting=0.4;

n_metingen=1024;

tijdmin=tijd_meting/n_metingen;

fs=2/tijdmin;

Code: Selecteer alles

F=[0 400 800 1000];	 % band limits

A=[0 1 0];

  % band type: 0='stop', 1='pass'

dev=[0.0001 10^(0.1/20)-1 0.0001]; % ripple/attenuation spec

[M,Wn,beta,typ]= kaiserord(F,A,dev,fs);		% window parameters

b=fir1(M,Wn,typ,kaiser(M+1,beta),'noscale');

hieruit kan je afleiden wat ik gebruikt heb voor fs de band waarin ik filter 400 tot 800 doorlaten de rimpel ...

Volgende grafiek geeft een idee van het signaal ervoor:

origineel 603 keer bekeken

erna geeft dit:

gefilterd 603 keer bekeken

zie in bijlage voor de waarde van elke burst van 1024.

De code geeft weer in detail wat ik heb uitgespookt met je data. Indien je graag voor elk frequentie bandje de waarden kent van de coëfficiënten kan ik je die gegeven dan kan je van al die bandjes apart je filter functie (dus die convolutie programmeren)

het automatiseren van het vinden van deze coëfficiënten denk ik dat nogal moeilijk is, groeten.

[Wimapon]

Hoi Bert,

Je krijgt zo te zien hetzelfde eruit als ik.

de enige sinus die in het filterbandje zit, staat helemaal links in de grafiek van de output.

Misschien is het goed om een frequency spectrogram van de output te laten zien.. Ik ben heel benieuwd of jij er hetzelfde uitkrijgt als ik

ik vind het vreemd dat het hele signaal eigenlijk weg is, en dat je alleen maar 3 of 4 sinussen overhoudt .... de rest van de waarden vallen

buiten de doorlaat en zijn dus voor mij niet meer interessant....

ik weet niet goed wat ik hiermee aan moet.

ik ben wel geinterresseerd in de coefficienten die ik kan gebruiken om de volgende frequentiebandjes ook te kunnen filteren.

800- 1200 Hz , 2000 - 2400 Hz en 2800 - 3200 Hz

[Gringo]

Haai, ik heb laatst een FIR en een IIR moeten programmeren voor een opdracht van school in C.

//Fir.c FIR filter. Include coefficient file with length N

#include "LP5300_51.cof" //coefficient file

#include "dsk6713_aic23.h"

//codec-dsk support file

Uint32 fs=DSK6713_AIC23_FREQ_24KHZ; //set sampling rate

int yn = 0;

//initialize filter's output

short dly[N]; //delay samples

interrupt void c_int11() //ISR

{

short i;



dly[0]=input_sample(); //input newest sample

yn = 0; //initialize filter's output

for (i = 0; i< N; i++)

yn += (h * dly); //y(n) += h(i)* x(n-i)

for (i = N-1; i > 0; i--) //starting @ end of buffer

dly = dly[i-1]; //update delays with data move

output_sample(yn >> 15); //scale output filter sample

return;

}

void main()

{

comm_intr(); //init DSK, codec, McBSP

while(1); //infinite loop

}

Het belangrijkste stuk zijn de twee for lussen.

De eerste is de convulsie lus:

for (i = 0; i< N; i++)

yn += (h * dly);

De formule voor convulsie staat hier: http://nl.wikipedia.org/wiki/FIR. (sorry, tex kost me teveel moeite)

De tweede lus is de delay lus.



for (i = N-1; i > 0; i--)

dly = dly[i-1];

Hier wordt elke waarde in het delay een plekje verschoven. Wiskundig gezien y(n) wordt y(n-1), y(n-1) wordt y(n-2) enzovoort

Ik heb een coefficienten file lp5300_51.cof gemaakt waarin h gedefinieerd wordt met daarin de coefficienten en N met het orde nummer van het filter. In dit geval dus een lowpass filter met een cut off frequency van 5300Hz van de 51ste orde.

Via matlab is het vrij simpel om de coefficienten te creeeren. Helaas is mijn desktop met daarop matlab gecrashed dus ik kan je even niet helpen. Als je zelf de beschikking hebt over matlab, type dan "fdatool" op de command prompt van matlab en je komt in de filter design tool van matlab. Hieruit kun je makkelijk de coefficienten exporteren.

[Gringo]

Ik hoop dat je hier wat mee kan. Ik vond zelf de FIR makkelijker dan de IIR. Ik heb ook nog wel mijn code die ik gebruikt heb voor de IIR. Grootste verschil is dat de FIR niet recursief is en de IIR wel. In mijn geval bestond de IIR dus uit twee sets van 4 lussen. Een voor het niet recursieve deel en eentje voor het recursieve deel.

Ik ben geen ster op het gebied van programmeren of filters, maar misschien dat dit je wat verder helpt. Als je nog vragen hebt hoor ik het wel.

[317070]

ik vind het vreemd dat het hele signaal eigenlijk weg is, en dat je alleen maar 3 of 4 sinussen overhoudt .... de rest van de waarden vallen

buiten de doorlaat en zijn dus voor mij niet meer interessant....

ik weet niet goed wat ik hiermee aan moet.

Wel, je moet ze niet verwaarlozen. Wat er gebeurt is dat je een deel van de energie in je signaal weggefilterd is, je resultaat bevat dus minder energie dan je oorspronkelijke signaal. Het gedeelte wat je wil hebben staat nu voor je, alleen is het uiteraard niet zo sterk als je volledige signaal. Ik zou dit niet te snel als 'zuivere ruis' bestempelen.

Die 3 golfjes in het begin zijn een deel van het overgangsverschijnsel. Hoe het komt dat ze hier zo sterk zijn weet ik ook niet.

ik ben wel geinterresseerd in de coefficienten die ik kan gebruiken om de volgende frequentiebandjes ook te kunnen filteren.

800- 1200 Hz , 2000 - 2400 Hz en 2800 - 3200 Hz

Kun je even je samplefrequentie meegeven? Ik heb de laatste keer op 6000Hz gegokt, maar ik zit er precies wel een beetje naast. Als je aan 6000Hz samplet, kun je bovendien (meestal) niets nuttigs zeggen over frequenties boven 3000Hz.