ISSN 0236-235X (P)
ISSN 2311-2735 (E)

Journal influence

Higher Attestation Commission (VAK) - К1 quartile
Russian Science Citation Index (RSCI)

Bookmark

Next issue

2
Publication date:
16 June 2024

Applying MATLAB in the design of digital filters for selecting Pc5 geomagnetic pulsations

Date of submission article: 10.11.2022
Date after edit article: 13.12.2022
Date of acceptance for publication: 15.12.2022
UDC: 537.86:621.396
The article was published in issue no. № 2, 2023 [ pp. 257-262 ]
Abstract:The paper considers a design procedure of an optimal nonrecursive bandpass digital filter with a finite impulse response (FIR-filter) using the MATLAB Filter Design tool and the method of best uniform (Chebyshev) approximation. The filter helps solving the problem of extracting Pc5 geomagnetic pulsations from a data set of geomagnetic field measurements. This type of pulsations was chosen due to the availability of 1-second data in a widely spaced network of ground-based geomagnetic observatories with standardized ground-based geophysical equipment. After proper processing, these data can be used in a detailed analysis of: the disturbance properties in the Earth's magnetic field in the range of long-period pulsations; the nature of the interaction of waves and particles in the magnetosphere. The results of this analysis can be used, for example, in calculating a space weather forecast, which makes this work relevant. The problem of selecting Pc5 ripples is solved by passing the original data set through a band-pass FIR filter with the required characteristics depending on the range of the Pc5 ripple period – 150÷600 seconds. Hereof it follows the limits of the bandwidth equal to 1.7÷6.7 mHz. The choice of a non-recursive band-pass FIR filter is due to the possibility of providing a linear phase-frequency characteristic that excludes phase distortions at the output of the FIR filter; stability by definition of this type of filters. The first condition also ensures that there are no requirements for the phase response of the FIR filter. The operability of the obtained digital filter is demonstrated on the example of processing a real set of measurement data of the geomagnetic field state obtained from the Lycksele geomagnetic observatory (Sweden, Geological Survey of Sweden, international IAGA code LYC), which is a part of the INTERMAGNET international network. The necessary information about this observatory is available on the Internet: https://www.intermagnet.org. Bandpass FIR filter design and calculations were carried out in MATLAB R2022b.
Аннотация:В статье рассмотрена процедура проектирования методом наилучшей равномерной (чебышевской) аппроксимации оптимального нерекурсивного полосового цифрового фильтра с конечно-импульсной характеристикой (КИХ-фильтра). Фильтр предназначен для решения задачи выделения геомагнитных пульсаций типа Pc5 из набора показателей измерений геомагнитного поля. Этот тип пульсаций выбран в связи с доступностью 1-секундных данных в широко разнесенной сети геомагнитных наземных обсерваторий, имеющих стандартизированную геофизическую наземную аппаратуру. Данные, обработанные соответствующим образом, можно использовать, например, при детальном анализе свойств возмущений магнитного поля Земли в диапазоне длиннопериодных пульсаций, а также характера взаимодействия волн и частиц в магнитосфере. Результаты анализа можно использовать, в частности, при расчете прогноза космической погоды, что обусловливает актуальность настоящей работы. Решение задачи выделения пульсаций Pc5 производится путем пропускания исходного набора данных через полосовой КИХ-фильтр, обладающий требуемыми характеристиками, зависящими от диапазона периода пульсаций Pc5, – 150–600 секунд. Отсюда следуют границы полосы пропускания, равные 1.7–6.7 мГц. Выбор нерекурсивного полосового КИХ-фильтра обусловлен возможностью обеспечения линейной фазочастотной характеристики, исключающей фазовые искажения на выходе КИХ-фильтра, а также устойчивостью по определению данного типа фильтров. Первое условие также гарантирует и отсутствие требований к фазочастотной характеристике КИХ-фильтра. Работоспособность полученного цифрового фильтра продемонстрирована на примере обработки реального набора данных измерений состояния геомагнитного поля, полученного с геомагнитной обсерватории Lycksele (Швеция), входящей в международную сеть INTERMAGNET. Проектирование полосового КИХ-фильтра и расчеты проводились в MATLAB R2022b.
Authors: Korobeynikov, A.G. (korobeynikov_a_g@mail.ru) - The National Research University of Information Technologies, Mechanics and Optics (Professor), St. Petersburg, Russia, Ph.D
Keywords: INTERMAGNET, PC5, geomagnetic ripple, bandpass filter, FIR-filter, Filter Design, matlab
Page views: 1283
PDF version article

Font size:       Font:

Система MATLAB используется при решении задач в различных предметных облас-тях [1–3]. Представленные в данной работе результаты относятся к области проектирования цифровых фильтров (ЦФ), параметры которых задаются исходя из свойств геомагнитных пульсаций (ГП), являющихся ультранизкочастотными колебаниями магнитного поля Земли (МПЗ) в диапазоне частот от тысячных долей герца до нескольких герц [4].

Актуальность исследования обусловлена ускоренным развитием научных направлений в области физики магнитосферы, где результаты анализа ГП служат для диагностики процессов в солнечном ветре и магнитосфере, прогноза динамики развития магнитосферных бурь, анализа свойств возмущений МПЗ в диапазоне длиннопериодных пульсаций, анализа харак-тера взаимодействия волн и частиц в магнитосфере [5, 6]. Полученные результаты можно использовать, например, при расчете прогноза космической погоды.

Аналогичные исследования начали проводиться довольно давно. Например, в [7] представлены основы цифровой фильтрации, методы которой применялись при обработке данных сейсмологии, гравиметрии и магнитометрии, а также описаны способы расчета и действие низкочастотных фильтров. В настоящее время разработка методов низкочастотной фильтрации геофизических данных, базирующихся на современных информационных технологиях, ведется достаточно интенсивно. Эти методы включают различные алгоритмы, например, вычисление ковариаций и линеаризации операторов прогноза и наблюдений, фильтр Калмана, ансамблевый алгоритм, метод частиц, байесовский подход [8–12]. Все это еще раз подчеркивает актуальность выбранной темы.

Исходные данные измерений состояния МПЗ в реальном масштабе времени с частотой дискретизации fд = 1 cек. взяты на сайте http://www.intermagnet.org [13].

Здесь имеются данные о более чем 150 геомагнитных обсерваториях, включенных в эту международную сеть. Для проверки полученного ЦФ использованы данные геомагнитной обсерватории Lycksele (Швеция) за 1 августа 2022 года.

В работе сделан акцент на выделение ГП типа Pc5, имеющих период колебаний в диапазоне Т = 150–600 сек. ГП типа Рс5 наиболее часто происходят во время возмущения магнитосферы. Что касается дневных ГП, то они обычно генерируются в моменты протекания буревых и суббуревых ночных процессов. Необходимо отметить, что при интенсификации ночной суббури амплитуда ГП Pc5 увеличивается.

За генерацию ГП Рс5 в процессе протекания суббури отвечает физический механизм интенсификации продольных токов и ионосферных электроструй. Необходимо отметить также отдельный класс ГП Рс5, за генерацию которых, вероятнее всего, отвечает внемагнитосферный механизм [14].

Вычислительные эксперименты были проведены с использованием системы MATLAB R2022b под управлением Windows 11. Их отличительной особенностью является фактическое отсутствие программного кода, то есть поставленная задача полностью решается при помощи специального инструментария. Это позволяет выполнять ее пользователям-геофизикам, не имеющим навыка написания исходного кода на MATLAB.

Постановка задачи

Проектирование оптимального нерекурсивного полосового ЦФ с конечно-импульсной характеристикой (КИХ) необходимо провести при помощи метода наилучшей равномерной (чебышевской) аппроксимации с учетом теоремы Котельникова [15]. Под оптимальным КИХ-фильтром понимается ЦФ минимально возможного порядка при заданных требованиях к амплитудно-частотной характеристике (АЧХ).

Зададим исходные данные к АЧХ проектируемого ЦФ. Так как метод синтеза КИХ-фильтра изначально предполагает линейность фазочастотной характеристики (ФЧХ), требования задаются только к АЧХ в основной полосе частот.

Исходя из заданного диапазона периода колебаний Т полоса пропускания имеет диапазон 1.7–6.7 мГц. Поэтому левая граничная частота полосы пропускания f-χ равна 1.7 мГц, а правая (fχ) – 6.7 мГц.

Значение левой граничной частоты полосы задержки f-k зададим равным 1.2 мГц, а правой (fk) – 7.0 мГц.

Значение максимально допустимого отклонения АЧХ от единицы δ1 в полосе пропускания зададим равным 5.8e-05. Тогда максимальное допустимое затухание amax будет равно 0.005 (дБ).

Значение минимально допустимого отклонения в левой (δ21) и правой (δ22) полосах задержки зададим равным 1e-5. Тогда максимально допустимое затухание в левой (a1min) и правой (a2min) полосах будет равно 100 (дБ).

Проектирование ЦФ

Запустив MATLAB, вызываем конструктор фильтров. Это можно сделать несколькими способами. Например, во вкладке APPS выбираем секцию signal processing and communications, а в ней – приложение Filter Designer. Такого же результата можно достичь, набрав в командной строке MATLAB команду filterDesig-ner.

В открывшемся меню необходимо выполнить следующие действия.

1. Выбрать полосовой тип фильтра. Для этого в Response Type нужно установить Bandpass.

2. Выбрать КИХ-фильтр. Для этого в Design Method установить FIR. В выпадающем списке надо оставить Equiripple.

3. В Filter Order выбрать Minimum Order. Величину Density Factor задать равной 20.

4. В Frequency Specifications задать следующие значения для параметров: Fs = 1, Fstop1 = = 0.0014, Fpass1 = 0.0017, Fpass2 = 0.007, Fstop2 = 0.0073.

5. В Magnitude Specifications, в окне Units, выбрать Linear и задать следующие значения параметров: Dstop1 = 1e-05, Dpass = 0.00058, Dstop2 = 1e-05.

6. Нажать Design Filter.

Вид инструмента Filter Designer, получае-мый в результате проделанных действий, пред-ставлен на рисунке 1, АЧХ спроектированного полосового цифрового КИХ-фильтра – на рисунке 2.

Для сохранения коэффициентов спроектированного цифрового КИХ-фильтра необходимо во вкладке File выбрать Generate MATLAB Code. Затем в выпадающем меню – одну из трех опций, например, Data Filtering Function (with System Objects). После этого бу-дет предложено выбрать имя функции. Назовем эту функцию, например, my_FIR_filter.m. Нажимаем на ввод. Выходим из конструктора фильтров.

Подготовка данных для проверки спроектированного цифрового КИХ-фильтра

Проверка спроектированного цифрового КИХ-фильтра была проведена на реальных данных, полученных 1 августа 2022 года на геомагнитной обсерватории Lycksele. С сайта http://www.intermagnet.org были скачаны соответствующие данные – файл lyc20220801psec.sec.

Выберем данные о X-компоненте МПЗ, то есть столбец под именем LYCX. Для этого во вкладке Home отыщем Import Data. В выпавшем меню выберем файл lyc20220801psec.sec. После загрузки этого файла появится окно с данными, где в разделе Range надо отметить D23:D8422. После этого в разделе Output Type выберем Numeric Matrix, затем в разделе Import Selection – Import Data. В рабочей области появится файл с данными об X-компоненте МПЗ. Для удобства этот файл можно переименовать. График данных можно получить разными способами. Например, выбрать вкладку PLOTS. Затем в рабочей области выделить my_massiv. После этого MATLAB предлагает различные типы графиков, например, график исходных данных (рис. 3).

На рисунке 3 видно, что данные достаточно чистые, то есть нет явных сбоев при их записи, а также пропусков. Кроме того, в его начале имеется достаточно большой выброс. Вероятнее всего, это связано с проявлением суббури.

Поскольку интерес представляют ГП, в исходных данных необходимо удалить постоянную составляющую (тренд) и провести нормализацию [16]. Для этого осуществляется переход в раздел Live editor, в нем выбирается вкладка Task, далее раздел Data preprocessing, а в нем Remove Trends. В выдающем меню подраздела Input data раздела Select data выбирается файл my_massiv, в выдающем меню подраздела Polynomial type раздела Specify trend parameters – Custom. В выпавшем меню выбирается 6. Тренд, представляющий себой полином 6-го порядка, выбирается на основе проведенных вычислительных экспериментов. Далее в разделе Display result галочкой отмечаются те данные, которые необходимо вывести на график.

Для нормализации данных снова в разделе Live editor выбирается вкладка Task, далее в разделе Data preprocessing – Normalize Data. В выдающем меню подраздела Input data раздела Select data выбирается файл detrendedData, далее в подразделе Normalization method раздела Specify method and parameters – Range. В появившемся меню выбирается диапазон [0, 1]. В разделе Display result галочкой отмечаются те данные, которые нужно вывести на график, например, Input data и Normalized data. На экране появляются соответствующие графики, а в рабочем пространстве – массив нормализованных исходных данных normalizedData.

Результаты применения спроектированного ЦФ КИХ-фильтра

Для фильтрации нормализованных данных необходимо выполнить команду my_massiv_filter_norm = my_FIR_filter (normalizedDa­ta). После этого в рабочем пространстве появ-ляется массив отфильтрованных данных – data my_massiv_filter_norm. Его можно вывести в графическом виде таким же способом, как и исходные данные. Графики нормализованных и фильтрованных данных начиная с 5-го часа и длиной в 1 час представлены на рисунках 4 и 5. Визуальный анализ этих графиков показывает практически полное отсутствие гармоник за пределами 150–600 сек., фазовый сдвиг ФЧХ (согласно рисунку 2), явное выделение ГП Pc5.

Выводы

Анализ полученных результатов позволяет сделать вывод о достаточно простой и эффективной процедуре применения инструментария MATLAB для предварительной обработки исходных данных. При этом профессиональное владе-ние системой MATLAB не требуется. На основе рассмотренной последовательности операций можно достаточно легко расширить множество типов выделяемых ГП с учетом их особенностей. Например, ГП Pc1 имеет период 0.2–5 сек., поэтому Pc1, согласно теореме Котельникова, не попадает в основную полосу частот. Для решения задачи необходимо применить процедуру интерполяции, например, при помощи выбора инструмента Signal Analyzer в разделе APPS.

Необходимо отметить, что спроектированный ЦФ имеет порядок, равный 15 737. Отсюда следует, что фильтрованные данные необ-ходимо брать после 15 737-го отсчета. Чтобы начинать работу с первого отсчета, можно, например, взять данные предыдущего дня. В случае обработки данных за несколько дней целесообразно использовать инструменты для работы с большими данными, также имеющиеся в MATLAB.

Полученные результаты могут являться основой для дальнейшей обработки и интерпретации данных, используемых, например, при прогнозировании космической погоды или анализе свойств возмущений МПЗ в диапазоне длиннопериодных пульсаций.

Список литературы

1. Коробейников А.Г. Применение искусственных нейронных сетей в системах автоматического управления магнитной левитацией // Программные продукты и системы. 2022. Т. 35. № 3. С. 452–457. doi: 10.15827/0236-235X. 139.452-457.

2. Korobeynikov A.G., Grishentsev A.Y., Velichko E.N., Korikov C.C., Aleksanin S.A., Fedosovskii M.E., Bondarenko I.B. Calculation of regularization parameter in the problem of blur removal in digital image. Optical Memory & Neural Networks, 2016, vol. 25, no. 3, pp. 184–191. doi: 10.3103/S1060992X16030036.

3. Дьяконов В.П. MATLAB и SIMULINK для радиоинженеров. М.: ДМК Пресс, 2011. 976 с.

4. Грунская Л.В. Электромагнетизм земной атмосферы. Владимир, 2019. 209 с. URL: https://dspace.www1. vlsu.ru/bitstream/123456789/8306/1/01934.pdf (дата обращения: 05.11.2023).

5. Леонович А.С., Мазур В.А. Линейная теория МГД-колебаний магнитосферы. М.: Физматлит, 2016. 480 с.

6. Clilverd M.A., Rodger C.J., Brundell J.B., Dalzell M., Martin I. et al. Long-lasting geomagnetically induced currents and harmonic distortion observed in New Zealand during the 7–8 September 2017 disturbed period. Space Weather, 2018, vol. 16, pp. 704–717. doi: 10.1029/2018SW001822.

7. Кулханек О. Введение в цифровую фильтрацию в геофизике; [пер. с англ.]. М.: Недра, 1981. 198 с.

8. Carrassi A., Bocquet M., Bertino L., Evensen G. Data assimilation in the geosciences: An overview of methods, issues and perspectives. WIREs Climate Change, 2018, vol. 9, no. 5, art. e535. doi: 10.1002/wcc.535.

9. Evensen G. Analysis of iterative ensemble smoother for solving inverse problems. Computational Geosciences, 2019, vol. 22, no. 3, pp. 885–908. doi: 10.1007/s10596-018-9731-y.

10. Potthast R., Walter A., Rhodin A. A localized adaptive particle filter within an operational NWP Framework. Monthly Weather Review, 2019, vol. 147, no. 1, pp. 345–362. doi: 10.1175/MWR-D-18-0028.1.

11. van Leeuwen P.J., Kunsch H.R., Nerger L., Potthast R., Reich S. Particle filters for high-dimensional geoscience applications: A review. Quarterly J. of the Royal Meteorological Society, 2019, vol. 145, no. 723, pp. 2335–2365. doi: 10.1002/qj.3551.

12. Климова Е.Г. Локальный ансамблевый алгоритм усвоения данных для нелинейных геофизических моделей // СибЖВМ. 2023. Т. 26. № 1. С. 27–42. doi: 10.15372/SJNM20230103.

13. Коробейников А.Г. Обработка и анализ данных с российского сегмента мировой сети магнитных обсерваторий ИНТЕРМАГНЕТ // Междунар. журнал гуманитарных и естественных наук. 2018. № 8. С. 91–98.

14. Potapov A.S., Polyushkina T.N., Pulyaev V.A. Observations of ULF waves in the solar corona and in the solar wind at the Earth's orbit. JASTP, 2013, vol. 102, pp. 235–242. doi: 10.1016/j.jastp.2013.06.001.

15. Макшанов А.В., Журавлев А.Е. Технологии интеллектуального анализа данных. СПб: Лань, 2019. 212 с.

16. Старовойтов В.В., Голуб Ю.И. Нормализация данных в машинном обучении // Информатика. 2021. T. 18. № 3. С. 83–96. doi: 10.37661/1816-0301-2021-18-3-83-96.

Reference List

1. Korobeynikov, A.G. (2022) ‘Applying artificial neural networks in automatic control systems for magnetic levitation’, Software & Systems, 35(3), pp. 452–457. doi: 10.15827/0236-235X.139.452-457 (in Russ.).

2. Korobeynikov, A.G., Grishentsev, A.Y., Velichko, E.N., Korikov, C.C., Aleksanin, S.A., Fedosovskii, M.E., Bonda-
renko, I.B. (2016) ‘Calculation of regularization parameter in the problem of blur removal in digital image’, Optical Memory & Neural Networks, 25(3), pp. 184–191. doi: 10.3103/S1060992X160300361.

3. Dyakonov, V.P. (2016) MATLAB and SIMULINK for Radio Engineers, Moscow, 976 p. (in Russ.).

4. Grunskaya, L.V. (2019) Electromagnetism of the Earth's Atmosphere, Vladimir, 209 p., available at: https://dspace.
www1.vlsu.ru/bitstream/123456789/8306/1/01934.pdf
(accessed November 05, 2023) (in Russ.).

5. Leonovich, A.S., Mazur, V.A. (2016) Linear Theory of MHD-Magnetosphere Oscillations, Moscow, 480 p. (in Russ.).

6. Clilverd, M.A., Rodger, C.J., Brundell, J.B., Dalzell, M., Martin, I. et al. (2018) ‘Long-lasting geomagnetically induced currents and harmonic distortion observed in New Zealand during the 7–8 September 2017 disturbed period’, Space Weather, 16, pp. 704–717. doi: 10.1029/2018SW001822.

7. Kulhánek, O. (1976) ‘Introduction to digital filtering in geophysics’, in Developments in Solid Earth Geophysics, 8, pp. 1–168. doi: 10.1016/c2009-0-07397-8. (Russ. ed.: Moscow (1981), 198 p.).

8. Carrassi, A., Bocquet, M., Bertino, L., Evensen, G. (2018) ‘Data assimilation in the geosciences: An overview of methods, issues and perspectives’, WIREs Climate Change, 9(5), art. e535. doi: 10.1002/wcc.535.

9. Evensen, G. (2019) ‘Analysis of iterative ensemble smoother for solving inverse problems’, Computational Geosciences, 22(3), pp. 885–908. doi: 10.1007/s10596-018-9731-y.

10. Potthast, R., Walter, A., Rhodin, A. (2019) ‘A localized adaptive particle filter within an operational NWP Framework’, Monthly Weather Review, 147(1), pp. 345–362. doi: 10.1175/MWR-D-18-0028.1.

11. van Leeuwen, P.J., Kunsch, H.R., Nerger, L., Potthast, R., Reich, S. (2019) ‘Particle filters for high-dimensional geoscience applications: A review’, Quarterly J. of the Royal Meteorological Society, 145(723), pp. 2335–2365. doi: 10.1002/qj.3551.

12. Klimova, E.G. (2023 ) ‘A local ensemble data assimilation algorithm for nonlinear geophysical models’, Numerical Analysis and Applications, 26(1), pp. 27–42. doi: 10.15372/SJNM20230103 (in Russ.).

13. Korobeynikov, A.G. (2018 ) ‘Processing and analysis of data from the Russian segment of the world network of magnetic observatory INTERMAGNET’, Int. J. of Humanities and Natural Sci., (8), pp. 91–98 (in Russ.).

14. Potapov, A.S., Polyushkina, T.N., Pulyaev, V.A. (2013) ‘Observations of ULF waves in the solar corona and in the solar wind at the Earth's orbit’, JASTP, 102, pp. 235–242. doi: 10.1016/j.jastp.2013.06.001.

15. Makshanov, A.V., Zhuravlev, A.E. (2019) Data Mining Technologies, St. Petersburg, 212 p. (in Russ.).

16. Starovoytov, V.V., Golub, Yu.I. (2021) ‘Data normalization in machine learning’, Informatics, 18(3), pp. 83–96. doi: 10.37661/1816-0301-2021-18-3-83-96 (in Russ.).


Permanent link:
http://swsys.ru/index.php?page=article&id=4999&lang=en
Print version
The article was published in issue no. № 2, 2023 [ pp. 257-262 ]

Perhaps, you might be interested in the following articles of similar topics: