using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using System.IO;
using System.Xml;
using ZedGraph;
using System.Collections;
using dnAnalytics.LinearAlgebra;
namespace Auger
{
public partial class Form1 : Form
{
public step_n final_step;
public step_0 STEP_0;//deklaracja zmiennej do wczytania danych
public chart CHART;//deklaracja zmiennej do tworzenia wykresów
public Form1()
{
InitializeComponent(); //inicjalizacja komponentów API
Rectangle r = Screen.PrimaryScreen.WorkingArea;//prostokąt który przechowuje wymiary obszaru roboczego ekranu
this.StartPosition = FormStartPosition.Manual;//ręczne ustalanie pozycji startowej okna
this.Size = new Size(r.Width, r.Height);//rozmiar okna - pochodzi od Rectangle r
this.Location = new Point(r.X,r.Y);//pozycja okna - pochodzi od Rectangle r
zedGraphControl1.Location = new Point(this.numericUpDown1.Width + this.numericUpDown1.Location.X + 25, 0);//lokacja i rozmiar poszczególnych wykresów
zedGraphControl1.Size = new Size(this.Size.Width - zedGraphControl1.Location.X - 100, this.Size.Height / 2-20);
zedGraphControl2.Location = new Point(zedGraphControl1.Location.X, zedGraphControl1.Location.Y + zedGraphControl1.Size.Height);
zedGraphControl2.Size = zedGraphControl1.Size;
zedGraphControl5.Size = zedGraphControl1.Size;
zedGraphControl5.Location = zedGraphControl1.Location;
zedGraphControl5.Visible = false;
zedGraphControl3.Visible = false;
zedGraphControl3.Location = zedGraphControl2.Location;
zedGraphControl3.Size = zedGraphControl2.Size;
zedGraphControl4.Location = zedGraphControl1.Location;
zedGraphControl4.Size = zedGraphControl1.Size;
zedGraphControl4.Visible = false;
zedGraphControl6.Size = zedGraphControl2.Size;
zedGraphControl6.Location = zedGraphControl2.Location;
zedGraphControl6.Visible = false;
button4.Location = new Point(zedGraphControl1.Location.X + zedGraphControl1.Size.Width + 5, zedGraphControl1.Size.Height / 2);
button5.Location = new Point(zedGraphControl1.Location.X + zedGraphControl1.Size.Width + 5, (int)(zedGraphControl1.Size.Height * 1.5));
button6.Location = new Point(button4.Location.X,zedGraphControl1.Location.Y+zedGraphControl1.Size.Height - (int)button6.Size.Height/2);
button4.Click += new EventHandler(przelacz_gore);
button5.Click += new EventHandler(przelacz_dol);
button6.Click += new EventHandler(przelacz_dol);
button6.Click += new EventHandler(przelacz_gore);
button7.Location = new Point(button4.Location.X,button1.Location.Y);
}
private void charts_initialization()
{
double [] X={};
double [] Y={};
chart my_chart = new chart(X, Y, zedGraphControl1, Color.Red, "Wyniki z tłem");
my_chart.remove();
my_chart.print();
my_chart = new chart(X, Y, zedGraphControl2, Color.Blue, "Po usunięciu tła");
my_chart.remove();
my_chart.print();
}
private void button1_Click(object sender, EventArgs e)//funkcja obsługująca kliknięcie przycisku 1
{
try//obsługa sytuacji wyjątkowych
{
STEP_0 = new step_0();//definicja zminnej do wczytywania
CHART = new chart(STEP_0.E, STEP_0.S, zedGraphControl1, Color.Red, "Wyniki z tłem");//podstawowe ustawienia wykresów
CHART.remove();//usunięcie poprzedniego wykresu, jeżeli był
CHART.print();//stworzenie nowego
}
catch (Exception ex)//jeśli wczytywanie danych zakończyło się błędem wyrzuci błąd
{
MessageBox.Show("Corrupted file format");
}
}
private void button2_Click(object sender, EventArgs e)//przycisk usun tlo
{
if (STEP_0 != null)//jeżeli istenieje wykres
remove(STEP_0.E, STEP_0.S, zedGraphControl2, "Po usunięciu tła");//remove - funkcja realizująca cały algorytm
else
{
MessageBox.Show("Input function doesn't exist!");//jeżeli wykresu nie mam wyskakuje info o błędzie
}
}
private void Form1_Load(object sender, EventArgs e)
{
}
// private void button3_Click(object sender, EventArgs e)//przycisk dodaj szum i usuń tło
// {
// if (STEP_0!=null)//jeżeli jest wykres
// {
// Random r = new Random();//zmienna do losowania
// double max = STEP_0.S[1];//szukanie maksymalnej wartości. Na początku przyjmujemy że pierwsza wartość jest największa
//for (int i = 1; i < STEP_0.E.GetLength(0); i++)//dla całej tablicy STEP_0
// if (STEP_0.S[i] > max) max = STEP_0.S[i];//jeżeli aktualna wartość jest większa od max, mamy nowego maxa
//double percent = (double)numericUpDown2.Value;//wczytanie wartości z pola procent szumu
//double [] S_noice = new double[STEP_0.E.GetLength(0)];//tworzy nową tablicę z zaszumionym sygnałem
//for (int i = 1; i < STEP_0.E.GetLength(0); i++)
// S_noice[i] = STEP_0.S[i] + (2 * r.NextDouble() - 1) * percent * max / 100;//dla każdej próbki dodaje szum z zakresu -procent * max do +procent * max
//remove(STEP_0.E, S_noice, zedGraphControl3, "Zaszumione po usunięciu tła");//usuwa tło zaszumionego sygnału i wyświetla w wykresie 3
//CHART.change_settings(STEP_0.E, S_noice, zedGraphControl4, Color.Pink, "Zaszumiony");//wyświetlenie przebiegu zaszumionego w wykresie 4
//CHART.remove();//usuwa stary wykres
//CHART.print();//rysuje
//}
// else
// MessageBox.Show("Input function doesn't exist!");//jeżeli funkcji nie ma wyświetla message
//}
private void przelacz_gore(object sender, EventArgs e)//kliknięcie w wykres 1
{
bool[] tab = new bool[3];
tab[0] = zedGraphControl1.Visible;
tab[1] = zedGraphControl4.Visible;
tab[2] = zedGraphControl5.Visible;
bool buf;
buf= tab[0];
tab[0] = tab[2];
tab[2] = tab[1];
tab[1] = buf;
zedGraphControl1.Visible = tab[0];
zedGraphControl4.Visible = tab[1];
zedGraphControl5.Visible = tab[2];
}
private void przelacz_dol(object sender, EventArgs e)
{
bool[] tab = new bool[3];
tab[0] = zedGraphControl2.Visible;
tab[1] = zedGraphControl3.Visible;
tab[2] = zedGraphControl6.Visible;
bool buf;
buf = tab[0];
tab[0] = tab[2];
tab[2] = tab[1];
tab[1] = buf;
zedGraphControl2.Visible = tab[0];
zedGraphControl3.Visible = tab[1];
zedGraphControl6.Visible = tab[2];
}
private void remove(double[] E, double[] S, ZedGraphControl control, string title)//usunięcie tła z przebiegu
{
step_1 STEP_1 = new step_1(E, S);//krok pierwszy wyniki do STEP_1
step_2 STEP_2 = new step_2(STEP_0.E, STEP_1.S_1);//krok drugi
step_3 STEP_3 = new step_3(STEP_0.E, STEP_2.S_2, STEP_1.S_1);//krok trzeci
ArrayList steps = new ArrayList();//kolekcja wyników kolejnych kroków
steps.Add(new step_n(STEP_0.E, STEP_3.S_3, (double)numericUpDown1.Value));//dodanie wyników króku czwartego
while (((step_n)steps[steps.Count - 1]).dalej)//dopóki sygnał nie zostanie odfiltrowany wykonuje następne kroki
steps.Add(new step_n(STEP_0.E, ((step_n)steps[steps.Count - 1]).S_n, (double)numericUpDown1.Value));
CHART.change_settings(E, ((step_n)steps[steps.Count - 1]).S_n, control, Color.Blue, title);//wyśiwetlenie odfiltrowanego sygnału
CHART.remove();
CHART.print();
SNR(STEP_1.S_1, ((step_n)steps[steps.Count - 1]).S_n);
double[] diff = new double[((step_n)steps[steps.Count - 1]).S_n.Length];
for (int i = 1; i < ((step_n)steps[steps.Count - 1]).S_n.Length - 1; i++)
{
double xp = (double)((step_n)steps[steps.Count - 1]).S_n[i + 1];
double xm = (double)((step_n)steps[steps.Count - 1]).S_n[i - 1];
diff[i-1] = ((xp-xm)/(E[i+1]-E[i-1]));
diff[0] = diff[1];
diff[diff.Length - 1] = diff[diff.Length - 2];
// MessageBox.Show(((xp - xm) / (E[i + 1] - E[i - 1])).ToString());
}
CHART.change_settings(E, diff, this.zedGraphControl6, Color.Red, "Odszumione po zrozniczkowaniu");
CHART.remove();
CHART.print();
final_step = (step_n)steps[steps.Count - 1];
}
private void SNR(double[] signal, double[] signal_filetred)//obliczenie SNR
{
if ((signal == null) || (signal_filetred == null)) return;//jeżeli tablice są puste nie wykonuje
if (signal.GetLength(0) != signal_filetred.GetLength(0)) return;//jeżeli długość sygnału odfiltrowanego i sygnału zaszumionego jest różna nie wykonuje
double average_signal = 0;//średnia sygnału
double average_noice = 0;//średnia szumu
for (int i = 0; i < signal.GetLength(0); i++)//liczy sumę kwadratów wartości szumu i sygnału
{
average_signal += signal_filetred[i] * signal_filetred[i];
average_noice += (signal[i] - signal_filetred[i]) * (signal[i] - signal_filetred[i]);
}
average_signal = Math.Sqrt(average_signal);//piewriastek
average_noice = Math.Sqrt(average_noice);
average_noice /= signal_filetred.GetLength(0);//średnia
average_signal /= signal_filetred.GetLength(0);
//textBox2.Text = (average_signal / average_noice).ToString();//wyświetla stosunek sygnału
//textBox1.Text = (20 * Math.Log10((average_signal / average_noice))).ToString();//wyświetla stosunek sygnału w dB
}
private void button7_Click(object sender, EventArgs e)
{
try
{
button1_Click(sender, e);
//button3_Click(sender, e);
}
catch (Exception ex)//jeśli wczytywanie danych zakończyło się błędem wyrzuci błąd
{
//MessageBox.Show("Corrupted file format");
}
}
private void button5_Click(object sender, EventArgs e)
{
}
private void button8_Click(object sender, EventArgs e)
{
step_n final;
try
{
final = final_step;
smooth(final,(int)Math.Pow(2,(double)((int)numericUpDown4.Value)),((int)numericUpDown3.Value));
MessageBox.Show("koniec");
}
catch (Exception ex)
{
MessageBox.Show("There is no data available to proceed");
return;
}
}
private void smooth(step_n final, int scope,int degree)
{
try
{
double[] finish = new double[final.S_n.Length];
for (int ii = 0; ii < final.S_n.Length; ii++)
{
if (ii < scope || ii > final.S_n.Length - scope)
{
finish[ii] = final.S_n[ii];
continue;
}
DenseMatrix m = new DenseMatrix(2 * scope + 1, degree + 1);
for (int i = 0; i < m.Rows-1; i++)
for (int j = 0; j < m.Columns; j++)
{
double pod = STEP_0.E[ii - scope + i];
if (i == scope) pod = -pod;
m[i, j] = Math.Pow(pod, j);
}
double[,] function = new double[2 * scope + 1, 1];
for (int s = 0; s < 2 * scope ; s++)
function[s, 0] = final.S_n[ii - scope + s];
DenseMatrix first, second, third,pomocnicza;
pomocnicza = (DenseMatrix)(m.Transpose().Multiply((Matrix)m));
first = pomocnicza;
try
{
first = (DenseMatrix)(m.Transpose().Multiply((Matrix)m)).Inverse();
second = (DenseMatrix)((Matrix)m.Transpose().Multiply((Matrix)new DenseMatrix(function)));
third = (DenseMatrix)first.Multiply(second);
}
catch (Exception ex)
{
first = (DenseMatrix)(m.Transpose().Multiply((Matrix)m));//.Inverse();
second = (DenseMatrix)((Matrix)m.Transpose().Multiply((Matrix)new DenseMatrix(function)));
third = (DenseMatrix)first.Multiply(second);
}
double [] coef = new double[2*scope+1];
for(int w=0;w<2*scope+1;w++)
{
double sum = 0;
for (int ww = 0; ww < degree + 1; ww++)
{
sum += first[0, ww] * Math.Pow(function[w, 0], ww);
}
coef[w] = sum;
}
double g = 0;
for (int fun = 0; fun < function.Length; fun++)
{
g += coef[fun] * function[fun,0];
}
finish[ii] = g;
}
CHART.change_settings(STEP_0.E, finish, this.zedGraphControl5, Color.Red, "Wygladzone");
CHART.remove();
CHART.print();
}
catch (Exception ex)
{
MessageBox.Show(ex.Message + " " + ex.Source + " " + ex.StackTrace);
}
}
private void numericUpDown3_ValueChanged(object sender, EventArgs e)
{
if ((int)numericUpDown3.Value % 2 == 1)
numericUpDown3.Value = numericUpDown3.Value - 1;
}
private void zedGraphControl1_Load(object sender, EventArgs e)
{
}
private void textBox1_TextChanged(object sender, EventArgs e)
{
}
private void textBox2_TextChanged(object sender, EventArgs e)
{
}
//private void groupBox4_Enter(object sender, EventArgs e)
//{
//}
//private void groupBox3_Enter(object sender, EventArgs e)
//{
//}
}
public class chart//klasa wykresu
{
private double[] X_axis;//tablica x-ów
private double[] Y_axis;//tablica y-ków
private ZedGraphControl control;//wykres
private Color color;//kolor
private String figure_title;//tytuł wykresu
public chart(double[] X_axis, double[] Y_axis, ZedGraphControl control, Color color, String figure_title)//konstruktor
{
if (X_axis.GetLength(0) != Y_axis.GetLength(0))
{
MessageBox.Show("Unequal axis dimensions");
return;
}
this.X_axis = X_axis;
this.Y_axis = Y_axis;
this.control = control;
this.color = color;
this.figure_title = figure_title;
}
public void change_settings(double[] X_axis, double[] Y_axis, ZedGraphControl control, Color color, String figure_title)
{
if (X_axis.GetLength(0) != Y_axis.GetLength(0))
{
MessageBox.Show("Unequal axis dimensions");
return;
}
this.X_axis = X_axis;
this.Y_axis = Y_axis;
this.control = control;
this.color = color;
this.figure_title = figure_title;
}
public void change_X_axis(double[] X_axis)
{
if (X_axis.GetLength(0) != Y_axis.GetLength(0))
{
MessageBox.Show("Unequal axis dimensions");
return;
}
this.X_axis = X_axis;//ustawia tablicę x-ów
}
public void change_Y_axis(double[] Y_axis)
{
if (X_axis.GetLength(0) != Y_axis.GetLength(0))
{
MessageBox.Show("Unequal axis dimensions");
return;
}
this.Y_axis = Y_axis;//ustawia tablicę y-ków
}
public void change_control(ZedGraphControl control)
{
this.control = control;//zmienia obszar wykresu
}
public void change_line_color(Color color)
{
this.color = color;//zmienia kolor linii
}
public void change_chart_title(String figure_title)
{
this.figure_title = figure_title;//zmienia tytuł wykresy
}
public void print()//wyświetla wylres
{
GraphPane myPane = control.GraphPane;
myPane.Title.Text = figure_title;//tytuł
myPane.XAxis.Title.Text = "E/eV";//opis osi
myPane.YAxis.Title.Text = "d(E*N(E))/dE";
myPane.XAxis.IsVisible = true;//widzialność osi
myPane.YAxis.IsVisible = true;
PointPairList list = new PointPairList();//para punktów
for (int ii = 1; ii < X_axis.GetLength(0); ii++)
{
list.Add(X_axis[ii], Y_axis[ii]);//dodanie punktu do wykresu
}
LineItem myCurve = myPane.AddCurve("", list, color, SymbolType.Circle);//nowa linia wykresu
myCurve.Line.Width = 0.1f;//
myCurve.Symbol.Size = 0.3f;
control.IsShowHScrollBar = true;
control.IsShowVScrollBar = true;
control.IsAutoScrollRange = true;
control.GraphPane = myPane;
control.AxisChange();//wyświetelnie i dopasowanie osi
}
public void remove()
{
GraphPane myPane = control.GraphPane;
myPane.CurveList.Clear();//usuwa linie z wykresu
}
}
public class step_0//krok 0
{
public double[] E;//tablica z energią
public double[] S;//tablica z próbkami
public step_0()//konstruktor
{
OpenFileDialog open = new OpenFileDialog();//okinenko z wyborem pliku
open.Filter = "pliki Xml (*.xml)|*.XML";//wczytuj tylko pliki xml
open.Multiselect = false;//nie można zaznaczyć wielu plików
if (open.ShowDialog() != DialogResult.OK) return;//pokazuje okno i jeżeli jest ok dalej
XmlReaderSettings settings = new XmlReaderSettings();//klasa odczytu xml
settings.ValidationType = ValidationType.None;
XmlReader r = XmlReader.Create(open.FileName, settings);
int ile_serii = 0, ile_pomiarow = 0;
while (r.Read())//odczyt pliku xml
{
if (r.Name == "pomiary")
{
Int32.TryParse(r.GetAttribute("ile_serii"), out ile_serii);
Int32.TryParse(r.GetAttribute("ile_pomiarow"), out ile_pomiarow);
break;
}
}
E = new double[ile_pomiarow + 1];
S = new double[ile_pomiarow + 1];
int i = 0;
while (r.Read())
{
if (r.Name == "pomiar")
{
Int32.TryParse(r.GetAttribute("nr"), out i);
if (i > E.GetLength(0) - 1)
{
break;
}
E[i] = Double.Parse(r.GetAttribute("x"));
S[i] = Double.Parse(r.GetAttribute("y"));
}
}
}
}
public class step_1//krok 1
{
public double[] S_1;//tablica z wynikiem działania
public step_1(double[] E, double[] S)//konstruktor
{
if(E.GetLength(0)!= S.GetLength(0))return;//czy E i S mają równą długość
S_1 = new double[E.GetLength(0)];
double min = S[1];//szukanie wartości minimalnej na początek ustalamy że pierszy element S jest namniejszy
for (int i = 1; i < E.GetLength(0); i++)
if (S[i] < min) min = S[i];//jeżeli aktualny element jest mniejszy od min to min = aktualny element
for (int i = 1; i < E.GetLength(0); i++)
S_1[i] = S[i] - min;//dla każdego elementu tablic S odejmujemy od niego min i wpisujemy do S_1
}
}
public class step_2
{
public double[] S_2;//tablica wyników
public step_2(double[] E, double[] S_1)//konstruktor
{
if (E.GetLength(0) != S_1.GetLength(0)) return;
S_2 = new double[E.GetLength(0)];
int k = 1;
for (; ; k++)
if (S_1[k] == 0) break;//Szukam pierwszego elementu dla którego S_1 jest równy zero k - indeks tego elementu
if (k == 1) S_1.CopyTo(S_2, 0);//jeżeli znaleziony element jest pierwszy
else if (k == 2)//jeżeli znaleziony element jest drugi
{
S_1.CopyTo(S_2, 0);
S_2[1] = 0;
}
else//jeżeli jest dalej
{
double alfa = oblicz(S_1, E, 2, k);//oblicza wsp alfa
for(int i = 2;i<k;i++)
{
double a = oblicz(S_1, E, i, k);//oblicz alfa
if (alfa < a) alfa = a;//znajduje wartść największą
}
if (alfa < 0) alfa = 0;//jeżelki alfa jest mniejsza od zero alfa równa się zero
S_1.CopyTo(S_2, 0);
for (int i = 1; i < k; i++)
S_2[i] = S_1[i]-S_1[1]*Math.Pow(((E[k]-E[i])/(E[k]-E[1])),alfa);//dla każdego elementu od 1 do k wykonaie wzoru
}
}
double oblicz(double[] S_1,double [] E,int i, int k)
{
return (Math.Log((S_1[i] / S_1[1]), Math.E)) / (Math.Log(((E[k] - E[i]) / (E[k] - E[1])), Math.E));//wzór z algorytmu
}
}
public class step_3
{
public double[] S_3;//wynik kroku 3
public step_3(double[] E, double[] S_2, double [] S_1)//konstruktor
{
if (E.GetLength(0) != S_2.GetLength(0)) return;
S_3 = new double[E.GetLength(0)];
int k = E.GetLength(0)-1;
for (; ; k--)
if (S_1[k] == 0) break;//szukanie ostatniego elementu dla którego S_1 jest równy 0
if (k == E.GetLength(0) - 1) S_3 = S_2;//jeżeli znaleziony element jest ostatni
else if (k == E.GetLength(0) - 2)//jeżeli przedostani
{
S_2.CopyTo(S_3, 0);
this.S_3[E.GetLength(0) - 1] = 0;
}
else//jeżeli dalej
{
S_2.CopyTo(S_3, 0);
double alfa = znajdz_max(S_1, E, k/*, j*/);//maksymalny współczynnik alfa
if (alfa < 0) alfa = 0;//jeżeli alfa <0
for (int j = k + 1; j < E.GetLength(0); j++)
{
S_3[j]=S_2[j]-S_2[E.GetLength(0)-1]*Math.Pow(((E[j]-E[k])/(E[E.GetLength(0)-1]-E[k])),alfa);//wzór z algorytmu
}
}
}
public double znajdz_max(double[] S_1, double[] E, int k/*, int j*/)//znajduje maksymalny wsp alfa
{
double max = oblicz(S_1, E, k + 1, k/*, j*/);//oblicz wsp alfa
for (int i = k + 1; i < E.GetLength(0) - 1; i++)
{
double mala = oblicz(S_1, E, i, k/*, j*/);//dla każdego elementu od k do końca
if (mala > max) max = mala;//szukanie maksymalnej wartości alfa
}
return max;
}
public double oblicz(double[] S_1, double[] E, int i, int k/*, int j*/)
{
return (Math.Log((S_1[i] / S_1[E.GetLength(0) - 1]), Math.E)) / (Math.Log(((E[i] - E[k]) / (E[E.GetLength(0) - 1] - E[k])), Math.E));//wzór z algorytmu
}
}
public class step_n
{
public double[] S_n;//wynik z kroku n
public int n,ktory=0;
public bool dalej = false;
public step_n(double [] E,double [] S_n_minus1,double delta)//konsytruktor
{
int l=1, m=1;
ArrayList zeros = new ArrayList();//tablica o dowolnym rozmiarze i nieokroślonym typie przechowująca elementy S_n_minus_1 które są równe 0
for (int i = 1; i < E.GetLength(0); i++)
if (S_n_minus1[i] == 0) zeros.Add(i);//dodawnie do tablicy lementu, który równy jest zero
for (int i=0;i<zeros.Count-1;i++)//przeglądanie elementów tablicy S_n_minus_1 i sprawdzenie czy spełniają warunki kontynuowania określone w algorytmie
{
if (((E[(int)zeros[i + 1]] - E[(int)zeros[i]])>delta) && ((int)zeros[i + 1] > ((int)zeros[i] + 1/*indeks następnego zera większy o 2 od popredniego*/)))
{
l = (int)zeros[i];//piersze zero
m = (int)zeros[i + 1];//następne zero
dalej = true;
break;
}
else dalej = false;//jeżeli nie znaleziono punktów które mogą być przetwarzane dalej = false co oznacza że pętla zostaje przerwana
}
S_n = new double[E.GetLength(0)];
S_n_minus1.CopyTo(S_n, 0);
if (dalej)
{
double A = S_n_minus1[l + 1] / ((E[l + 1] - E[l]) * (E[m] - E[l + 1]));//wsp A
ktory = l + 1;
for (int i = l + 1; i < m; i++)
{
double min = S_n_minus1[i] / ((E[i] - E[l]) * (E[m] - E[i]));
if (min < A)
{
// ktory = i;
A = min;//minimalny wsp A od l+1 do k
}
}
int ile = 0;
for (int i = l + 1; i < m; i++)
{
S_n[i] = S_n_minus1[i] - A * (E[i] - E[l]) * (E[m] - E[i]);//przetworenie tablicy ze wzoru
if (S_n[i] < 1e-16)
{
S_n[i] = 0;//rozwiązanie problemu niedokładności obliczeń w komputerze
ile++;
}
}
}
}
}
}
Copyright © 2008-2010 EPrace oraz autorzy prac.