ブログ

割とコンピュータよりの情報をお届けします。

TakeMe

Python numpyでFFTを実行する.周波数軸はどうしたら

numpyでFFTを実行する場合,周波数軸はnumpy.linspaceで設定するのが便利である.
よーく見るとずれてしまって困った話.

まず周波数解析の例を示す.周波数軸を作る部分は
np.linspace(0, サンプリング周波数, データ長さ)
で大体あっているのだが,微妙にずれる.
このずれの分が気になったのでよーく考えると以下の例の場合には
freq = np.linspace(0, fsmp - 1/(t_len + t[1] - t[0]), num=len(t));
ということになる.

ピったしの周波数軸が作れる.

import numpy as np
import matplotlib.pyplot as plt

# 例の信号を作成
t = np.linspace(0.001, 4, 4000);
z = 0.1 + 0.2 * np.sin(t * 10 * 2 * np.pi) + 0.2 * np.sin(t * 33 * 2 * np.pi);

# サンプリング周波数
fsmp = 1 / (t[1] - t[0]);
# 解析時間
t_len = t.max() - t.min();

Fz = np.fft.fft(z) / z.shape[0] * 2; # 折り返すのでパワーが2分の1になっている.
Fz[0] = Fz[0] / 2; # 平均成分は折り返さない.
Fz_abs = np.abs(Fz);

freq = np.linspace(0, fsmp, num=len(t));

fname = 'text.csv';
np.savetxt(fname, np.array([freq, np.real(Fz), np.imag(Fz)]).transpose(), fmt='%.18g', delimiter=',');

plt.figure(1)
plt.subplot(311)
plt.plot(t, z)

plt.subplot(312)
plt.plot(freq, Fz_abs)
plt.xlim([0, fsmp/2]);

#136/4
plt.subplot(313)
plt.plot(range(len(t)), Fz_abs)
#plt.xlim([0, 500]); # 後半(高周波側)は前半の折り返し

plt.show();

≫ 続きを読む

2019/01/23 コンピュータ   TakeMe
タグ:Python

Pythonでは複素数の演算が行えることを知った

最近の言語は複素数が使えることが多いが,Pythonも複素数演算が扱える.

例えば以下のような演算ができる.
以下の例では虚部と実部を分けて指定する方法と極座標表示から得る方法の例.

import cmath

# 虚部と実部を分けて指定
a1 = complex(1, 0);
a2 = complex(2, 0);

# 偏角を使った指定
mag1 = 1.0;
arg1 = 0.0;
a1 = complex(mag1 * cmath.cos(arg1), mag1 * cmath.sin(arg1))

mag2 = 2.0;
arg2 = 0.0;
a2 = complex(mag2 * cmath.cos(arg2), mag2 * cmath.sin(arg2))


trial = complex(100, 0);

a21 = a2 - a1;
b = a21 / trial

print(b)

いったん複素数になってしまえば掛け算,割り算は可能.
(例の偏角はラジアンなので注意)

≫ 続きを読む

2019/01/23 コンピュータ   TakeMe
タグ:Python

Python のプログラムはIDLEでもデバッグできる

PythonのGUIとしてWinPythonにはIDLEが添付している.
このIDLEはブレイクポイントなどの設定をどのように行うのか不明でデバッグには使えないと思っていたが,どうもShellの側のメニューバーにDebugという項目があり,Debuggerを選択すればブレイクポイントも設定できるようになるらしい.

IDLEを起動するとPython ShellというWindowが現れる.このとき,ファイルを開く方法は何ら難しいことはないと思う.
ファイルを開くとエディタの中にファイルが開く.
よくIDEを使った開発になれている方であればF5キーでプログラムが実行されるもの推測しやすいと思う.

ただ,ブレイクポイントはどのように設定するのか?

じつは,エディタの側ではなくPython Shellの側のメニューバーのDebug->Debuggerを選択してデバッグをオンにしておかなければならない.

これをオンにしておくとエディタで行を右クリックした時に,Set Breakpoint / Clear Breakpointという選択ができるようになっている.
毎回Visual Studio Codeをインストールしていた.もっと早く教えてよ.

参考
PythonのIDLEでデバッグを行う方法(GAMMSOFT: https://gammasoft.jp/python/debugging-python-idle/)

≫ 続きを読む

2019/01/18 コンピュータ   TakeMe
タグ:Python

SSD(Q-360)でプチフリ?LPM HIPM非対応

2011年製の中古パソコンのハードディスクをSSD Q-360に置き換えて半年.HDDに比べて起動時間が短縮.起動時間など気にならなくなっていた.

年明けになってプチフリとみられる現象が発生するようになった.

理由はよくわかっていないが,設定が原因ではないかと思うようになった.

調べてみるとQ-360はHIPMには対応していないようだ.

そこで,「PCの動作が突然止まってしまう!SSDの謎の不具合「プチフリ」の原因と対策は?」「Windows10のスリープ不安定問題 その1」を参考に設定をいじってみた.
これで設定がうまくいっているのかわからない....

レジストリエディタを使わないと設定できない項目だったのか

 

20190118 追記: ちなみにどうも私の場合にはプチフリーズの原因は違ったらしい.(遅いネットワーク経由でWindows Updateの準備が始まっていたらしい)

≫ 続きを読む

2019/01/12 コンピュータ   TakeMe

WinPython64-3.7.1.0ZeroでSavitzyky-Golayフィルタ

WinPython64-3.7.1.0Zeroを準備して参考ページ「SciPy で Savitzky-Golay フィルタ」を実行できるようにするまで記録する.

1. WinPythonをダウンロードする.

2. ダウンロードしたファイルを実行して適当なフォルダにインストールする.

3. WinPython Command Prompt.exeを実行して,
pip install scipy
を実行する.すると,numpyもインストールされる.
pip install matplotlib
も実行しておく.
4. IDLE (Python GUI).exeを実行して,参考ページからコードをコピーしてくる.
そしてF5キーでこれを実行してみる.
おそらく,
「TypeError: slice indices must be integers or None or have an __index__ method」とエラーになるはず.

5. y1 = signal.savgol_filter(y, n/4+1, 5)などとなっている部分の「/」が問題で「//」に置き換える必要がある.
(y1 = signal.savgol_filter(y, n//4+1, 5)とするということ)

最近のPythonでは「/」が実数を返してしまうということであった.

≫ 続きを読む

2019/01/10 コンピュータ   TakeMe

C#からExcel操作

C#からExcelを操作して計算結果をエクスポートする例を求められた.
 

レイトバインディングならdllへの依存性が減らせるが,プログラムの作成中に面倒なことにオートコンプリートが効かない.Microsoft.Office.Interop.Excelを参照に追加するアーリーバインディングが普通です.

using Excel = Microsoft.Office.Interop.Excel;

これが参考に載っていたやり方.ただし,using Microsoft.Office.Interopでもよいかも.

namespace app
{
    public partial class Form1 : Form
    {
        private Excel.Application objApp;
        private Excel.Workbook objWorkbook;
        private Excel.Worksheet objWorksheet;


        public Form1()
        {
            InitializeComponent();
        }

        private void button1_Click(object sender, EventArgs e)
        {
            objApp = new Excel.Application();
            objWorkbook = objApp.Workbooks.Add();
            objWorksheet = (Excel.Worksheet)objWorkbook.Sheets[1];
            Excel.Range objRange = objWorksheet.Range["A1..B2"];//objApp.get_Range("A1");
            object[,] tbl = new object[2, 2];
            tbl[0, 0] = "TEST";
            tbl[0, 1] = "TEST2";
            tbl[1, 0] = 0;
            tbl[1, 1] = 0.1;
            objRange.Value = tbl;
            objApp.Visible = true;
        }
    }
}

いったんここまで出来上がってからレイトバインディングに直すのが宜しい.
昔は,VBとC#では大きな違いがあったが,最近はC#でもdynamicなる属性が追加されており,InvokeMember関数を最小限にして置き換えができるようになった.

dynamic属性に付け替えExcel.Application()の部分を汎用のType.GetTypeFromProgID()に変更するだけでよい

    public partial class Form1 : Form
    {
        private dynamic objApp;
        private dynamic objWorkbook;
        private dynamic objWorksheet;

        public Form1()
        {
            InitializeComponent();
        }

        private void button1_Click(object sender, EventArgs e)
        {
            Type t = Type.GetTypeFromProgID("Excel.Application");
            objApp = Activator.CreateInstance(t);

            objWorkbook = objApp.Workbooks.Add();
            objWorksheet = objWorkbook.Sheets[1];
            dynamic objRange = objWorksheet.Range["A1..B2"];
            object[,] tbl = new object[2, 2];
            tbl[0, 0] = "TEST";
            tbl[0, 1] = "TEST2";
            tbl[1, 0] = 0;
            tbl[1, 1] = 0.1;
            objRange.Value = tbl;
            objApp.Visible = true;
        }
    }

参考

【Excelの起動から】 C#でExcelを操作する【保存まで】 ( Milk's Memo Note: https://www.milkmemo.com/entry/csharp_excel)

≫ 続きを読む

2019/01/05 コンピュータ   TakeMe

C#で配列のソート

C#で配列のソートをする方法としてArray.Sortがある.

ただし,引数に与えた配列は順番が破壊されてしまう.

通常は,以下のようにするが,これだと元のtestは書き換えられる.

            Int32[] test = new Int32[5];
            
            for (int i = 0; i < test.Length; i++)
            {
                test[i] = test.Length - i;
            }
            Array.Sort(test);

もとの配列の中でのインデックスを記録しておきたい場合,test1にインデックスを記録しておき,Array.Sortに二つ配列を渡してやると,test1はソート前の配列での順番を表すことになる.

            Int32[] test1 = new Int32[5];
            Int32[] test2 = new Int32[5];
            
            Random rand = new Random(100);
            for (int i = 0; i < test1.Length; i++)
            {
                test1[i] = i;
                test2[i] = rand.Next();
            }
            Array.Sort(test2, test1);
 

≫ 続きを読む

2018/12/28 コンピュータ   TakeMe

Windowsタスクバーを表示したり非表示にしたり

Windows Embedded Standardの装置でタスクバーを非表示にしている例があった.最初はスタートメニューにそのような設定項目があるのかと思っていたが,どうもそうではないらしい.どのように行っているのか確認していると...

サンプルとして以下のようなコードで実施できる.

#include "stdafx.h"
#include <Windows.h>

int main()
{
    HWND hWnd = FindWindow(_T("Shell_TrayWnd"), NULL);

    ShowWindow(hWnd, IsWindowVisible(hWnd) ? SW_HIDE : SW_SHOW);

    return 0;
}

大したことのないアプリである.以下のようなコードでもよいです(Windows 10ではこれで動いていた).

#include "stdafx.h"
#include "WindowsProject2.h"


int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
                     _In_opt_ HINSTANCE hPrevInstance,
                     _In_ LPWSTR    lpCmdLine,
                     _In_ int       nCmdShow)
{
    HWND hWnd = FindWindow(_T("Shell_TrayWnd"), NULL);
    ShowWindow(hWnd, IsWindowVisible(hWnd) ? SW_HIDE : SW_SHOW);
    return 0;
}

≫ 続きを読む

2018/12/05 コンピュータ   TakeMe

Visual Studio 2017ならWindows アプリケーションがテンプレートですぐに

猫でもわかるWindowsプログラミング を買ってからすでに,十数年たっているが,今でも十分に使える説明が多い。買った頃は,BCCでコンパイルしていた.

Windows Embedded Compact 用のアプリケーションもほぼ同じ説明で作れる

Visual Studio 2017なら最初の空のウィンドウ(メニューバー付き)を作るまでテンプレートで出来てしまう。

Visual Studio 2017ならプロジェクトを新規作成するときのダイアログで,Visual C++の項目Windows デスクトップ アプリケーションを選択すると,一気にウィンドウ作成まで行ってしまう.

説明との,違いといえば最近のバージョンのWindows(Windows Embedded Compactを含む)はWCHARが標準になっているところである.

つまり文字の処理関係が変わっている点だけ注意が必要である.しかし,説明自体は非常によくできている.

≫ 続きを読む

2018/12/02 コンピュータ   TakeMe

Visual C++プロジェクトでアセンブリ出力という機能

Visual C++でプロジェクトを作っていれば,アッセンブリ出力という機能が使えることを知った。

プロジェクトのプロパティで[C/C++]->[すべてのオプション]->アッセンブリの出力で「アッセンブリコード、コンピューター語、ソースコード(/FAcs)」を出力するようにすれば,

#include "stdafx.h"


int main()
{
    signed char a;
    short b;
    int c;

    a = -10;
    b = 20;
    c = 30;
    c = a + b + c;
    return 0;
}

このコードが,*.codファイルに出力される。

; Listing generated by Microsoft (R) Optimizing Compiler Version 19.13.26129.0 

    TITLE    E:\start001\ConsoleApplication1\ConsoleApplication1.cpp
    .686P
    .XMM
    include listing.inc
    .model    flat

INCLUDELIB MSVCRTD
INCLUDELIB OLDNAMES

PUBLIC    _main
EXTRN    __RTC_InitBase:PROC
EXTRN    __RTC_Shutdown:PROC
;    COMDAT rtc$TMZ
rtc$TMZ    SEGMENT
__RTC_Shutdown.rtc$TMZ DD FLAT:__RTC_Shutdown
rtc$TMZ    ENDS
;    COMDAT rtc$IMZ
rtc$IMZ    SEGMENT
__RTC_InitBase.rtc$IMZ DD FLAT:__RTC_InitBase
rtc$IMZ    ENDS
; Function compile flags: /Odtp /RTCsu /ZI
; File c:\start001\consoleapplication1\consoleapplication1.cpp
;    COMDAT _main
_TEXT    SEGMENT
_c$ = -32                        ; size = 4
_b$ = -20                        ; size = 2
_a$ = -5                        ; size = 1
_main    PROC                        ; COMDAT

  00000    55         push     ebp
  00001    8b ec         mov     ebp, esp
  00003    81 ec e4 00 00
    00         sub     esp, 228        ; 000000e4H
  00009    53         push     ebx
  0000a    56         push     esi
  0000b    57         push     edi
  0000c    8d bd 1c ff ff
    ff         lea     edi, DWORD PTR [ebp-228]
  00012    b9 39 00 00 00     mov     ecx, 57            ; 00000039H
  00017    b8 cc cc cc cc     mov     eax, -858993460        ; ccccccccH
  0001c    f3 ab         rep stosd

  0001e    c6 45 fb f6     mov     BYTE PTR _a$[ebp], -10    ; fffffff6H

  00022    b8 14 00 00 00     mov     eax, 20            ; 00000014H
  00027    66 89 45 ec     mov     WORD PTR _b$[ebp], ax

  0002b    c7 45 e0 1e 00
    00 00         mov     DWORD PTR _c$[ebp], 30    ; 0000001eH

  00032    0f be 45 fb     movsx     eax, BYTE PTR _a$[ebp]
  00036    0f bf 4d ec     movsx     ecx, WORD PTR _b$[ebp]
  0003a    03 45 e0     add     eax, DWORD PTR _c$[ebp]
  0003d    03 c8         add     ecx, eax
  0003f    89 4d e0     mov     DWORD PTR _c$[ebp], ecx

  00042    33 c0         xor     eax, eax

  00044    5f         pop     edi
  00045    5e         pop     esi
  00046    5b         pop     ebx
  00047    8b e5         mov     esp, ebp
  00049    5d         pop     ebp
  0004a    c3         ret     0
_main    ENDP
_TEXT    ENDS
END

x64の場合には以下のようになる

; Listing generated by Microsoft (R) Optimizing Compiler Version 19.13.26129.0 

include listing.inc

INCLUDELIB MSVCRTD
INCLUDELIB OLDNAMES

PUBLIC    main
EXTRN    _RTC_InitBase:PROC
EXTRN    _RTC_Shutdown:PROC
;    COMDAT pdata
pdata    SEGMENT
$pdata$main DD    imagerel $LN3
    DD    imagerel $LN3+82
    DD    imagerel $unwind$main
pdata    ENDS
;    COMDAT rtc$TMZ
rtc$TMZ    SEGMENT
_RTC_Shutdown.rtc$TMZ DQ FLAT:_RTC_Shutdown
rtc$TMZ    ENDS
;    COMDAT rtc$IMZ
rtc$IMZ    SEGMENT
_RTC_InitBase.rtc$IMZ DQ FLAT:_RTC_InitBase
rtc$IMZ    ENDS
;    COMDAT xdata
xdata    SEGMENT
$unwind$main DD    05051c01H
    DD    010a030dH
    DD    070030025H
    DD    05002H
xdata    ENDS
; Function compile flags: /Odtp /RTCsu /ZI
; File e:\start001\consoleapplication1\consoleapplication1.cpp
;    COMDAT main
_TEXT    SEGMENT
a$ = 4
b$ = 36
c$ = 68
main    PROC                        ; COMDAT

$LN3:
  00000    40 55         push     rbp
  00002    57         push     rdi
  00003    48 81 ec 28 01
    00 00         sub     rsp, 296        ; 00000128H
  0000a    48 8b ec     mov     rbp, rsp
  0000d    48 8b fc     mov     rdi, rsp
  00010    b9 4a 00 00 00     mov     ecx, 74            ; 0000004aH
  00015    b8 cc cc cc cc     mov     eax, -858993460        ; ccccccccH
  0001a    f3 ab         rep stosd

  0001c    c6 45 04 f6     mov     BYTE PTR a$[rbp], -10

  00020    b8 14 00 00 00     mov     eax, 20
  00025    66 89 45 24     mov     WORD PTR b$[rbp], ax

  00029    c7 45 44 1e 00
    00 00         mov     DWORD PTR c$[rbp], 30

  00030    0f be 45 04     movsx     eax, BYTE PTR a$[rbp]
  00034    0f bf 4d 24     movsx     ecx, WORD PTR b$[rbp]
  00038    8b 55 44     mov     edx, DWORD PTR c$[rbp]
  0003b    03 d0         add     edx, eax
  0003d    8b c2         mov     eax, edx
  0003f    03 c8         add     ecx, eax
  00041    8b c1         mov     eax, ecx
  00043    89 45 44     mov     DWORD PTR c$[rbp], eax

  00046    33 c0         xor     eax, eax

  00048    48 8d a5 28 01
    00 00         lea     rsp, QWORD PTR [rbp+296]
  0004f    5f         pop     rdi
  00050    5d         pop     rbp
  00051    c3         ret     0
main    ENDP
_TEXT    ENDS
END

≫ 続きを読む

2018/11/24 コンピュータ   TakeMe