ある日、友達から問題が送られてきました。
一枚の画像(pngファイル)がLINEで送られてきましたが開けません。
きっと画像の中に別のデータのファイルが埋め込まれているのだろうと思い、バイナリエディタで開いてみました。
するとマジックナンバーからwavファイルだとわかり、拡張子をwavに変えましたが何も聞こえませんでした。
解けないため聞いてみると、
このサイトで開くとわかるよとのことでした。
実際に開いてみると
なんと文字が出てきました。
なぜ音から文字が！？
どうやら
このサイトでダウンロードしたようです。
音で文字や画像を作れるなんて面白い！と思い、自分も作れないかやってみることにしました。
ちなみにこのwavファイルをpngに変えたのは、LINEで直接wavファイルを送ると余計な周波数の部分がカットされてしまうからですね。


まずは何をしようか。
とりあえず一定の周波数の音が作れないかやってみました。
上のサイトではpythonを使っていたのでpythonでやってみようと思いました。
最初はgoogleで「python　音を作る」とかいろいろ検索して頑張って作ろうとしていたのですが、ふと気が付きました。
この時代にはchatGPTがあるじゃん！！
それでchatGPTに聞いてみると案の定コードを書いてくれました。
実行するとちゃんと音を作ってくれました！！！
すごい！chatGPT！！
試しに5000 Hzの波形を作ってwavファイルとして保存してspectrum analyzerに入れてみました。
すると
バグりました。
なんじゃこれ。
周波数が全体に広がっている・・・！？
いきなり躓いた。。
友達に聞いてみると、wavファイルはint16だとヒントをくれました。
16bitというのは2の16乗＝65536
符号付きなら半分の32768で、0の分を考慮して最大値が32767となる。
要するに整数で最大値32767にしなければいけませんでした。
最大値で割って絶対値の最大をにし、それに32767をかけた整数にすると、

ちゃんと5000 Hzの線が出ました！
実はこの正規化とかだけでもすごい時間がかかったんですよねー。
しかしまだ最初にすごい幅広い周波数の線が入ってる・・・
実はこれが一番最後の課題となってきます。

一定の周波数の波を作れたので、今度は帯を作ろう。
帯を作るには異なる周波数を足し合わせればできそうだ。
じゃあ5000～10000 Hzの波を足し合わせて帯を作ってみよう！

！？
帯がない！？
うっすら四角が見えるけれど中身は空洞だ。
このときはこの周波数と時間のグラフの名前もどういう原理なのかも知りませんでした。
というか波形を扱うこと自体が高校の物理以来でした。。

友人によると、多くの異なる周波数の波を足すと打ち消しあっちゃうから、打ち消される前の最初の部分だけを取り出せばうまく行くよとのことでした。
なるほどと思い、0~0.01秒の範囲のsin波をつなげてみることにしました。
おおお。
なんかすごい広がっているけれど、一応5000～10000 Hzの波が強く見えてる。

広がっているものはどうすればいいのか。
友人はバンドパスフィルターを使えばよいと教えてくれました。
バンドパスフィルターもchatGPTが書いてくれました。
おおお。
ついに帯ができた。
実はこの帯を作るだけでもすごい時間がかかりました。
波形についての基礎的な知識が足りなかったり忘れていました。
そもそもsin(2πft)ってなに？なんで2π掛かってるの？とかそんなレベルでした。
ずっとtのところにπを含めた値を入れていて、なんでπの倍数入れてるのに周期が変なの？どゆこと？とか考えてました。
時間と位相がぐちゃぐちゃになっていたんですよね。
そのためいろいろ波形を作ってみて感覚をつかみました。


この帯を組み合わせれば「Hello World」の最初の「H」が描ける！
太くて短い波と細くて長い波を組み合わせて

ついに「H」が描けた！！

ここまでは帯の組み合わせで行けましたが、「e」はそうはいきません。
「e」からは小さい帯を同じ時間に複数を作らないとできません。
また「e」を手打ちで作るのは面倒なので、画像を読み込んで作ることにしました。


ここまで自分がどうやって作ってきたかを思い出しながらプログラムを作り直してたのですが、気力がなくなったのであとは言葉だけで説明します。。。


一番最初の”Hello World”の音をフーリエ変換すると17000～20000 Hzの範囲でした。
そのため同じ範囲で作ろうと思いました。
縦幅が3000 Hzで、横幅はどうしようかなと考えて0.5秒にしました。
これを一文字分の四角だとして、どのくらいの解像度で描くか考え、50×50にしました。
すると小さい帯は3000/50=60 Hz ×　0.5/50=0.01秒になります。

最初は60 Hzを1 Hzずつすらした波60個を足そうと思いましたが、60 Hzおきに1個の波でも十分そうでした。

白黒の画像を用意し、そこに文字を書いてRGBが(0,0,0)のとき波を足す感じです。
波は0.01秒で異なる周波数を足したものを追加していく感じで作りました。

なんかここでもいろいろ苦戦した気がしますが忘れました。
とりあえず頑張って
文字を書くことができました！！！
ここまで結構長かったです。
分かんないことだらけでやる気が全然でないときもあったのですが何とかここまで来ました。
これは17000～20000 Hzでバンドパスフィルターを掛けています。

まだ文字の周りに波が広がっているので、位相をずらしたらいい感じにならないかなと思いやってみました。
50種類の波を使っているので0.01秒を50分割してずらしました。


ずらすと波が打ち消されにくくなるのか強まって文字が少し見やすくなりました。
それでも文字の周りに波は残っています。

この邪魔な波をどうにかして消せないものか。
特定の波を消すにはどうすればいいんだろう。
バンドパスフィルターを組み合わせればできそう？
でも普通にバンドパスフィルターを使うと、同じ時間においては一つの帯しか残せないよな。
小さい帯ごとにバンドパスフィルターを掛けたものを足し合わせばきれいになるかも。
やってみたけど変わらなかった。
同じ時間に2つ以上の帯を残せるマルチバンドパスフィルターとかライブラリないかな。
ないっぽい？

他に方法はないものか。
フーリエ変換してどの周波数の波が含まれているかを確認して、余計な周波数の部分を0にして逆フーリエ変換すれば綺麗な文字を取り出せるんじゃね！？
てか最初から周波数を用意して逆フーリエ変換すればできるかも？
でもせっかくここまで作ったし、作った波に対してフーリエ変換しよう。

こんな感じで考えてフーリエ変換によって文字をきれいにしようと考えました。
chatGPTに書いてもらって何となく実装しましたがうまくいきません。
そこで単純な波においてちゃんと特定の周波数を消せるか実験しました。

フーリエ変換：
import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
duration = 1 # 音声の長さ（秒）
sampling_rate = 44100 # サンプリングレート
t = np.linspace(0.0, duration, int(duration * sampling_rate), endpoint=False) # 時間軸

# 周波数成分の定義
frequencies = [17000, 17060, 17120] # 取り除く周波数のリスト

# 複数の周波数成分を持つ波形を生成
waveform = np.sum([np.sin(2 * np.pi * frequency * t) for frequency in frequencies], axis=0)

# FFTを計算して周波数領域で操作
fft_data = np.fft.fft(waveform) # FFTを計算
freq = np.fft.fftfreq(len(fft_data), d=1.0 / sampling_rate) # 周波数軸を生成
amp = np.abs(fft_data / (len(waveform) / 2))



remove_frequencies = [17060]
# 特定の周波数を取り除く（ゼロにする）
for frequency in remove_frequencies:
index = np.abs(freq - frequency).argmin() # 指定周波数に最も近いインデックスを取得
 fft_data[index] = 0.0 # インデックスに対応する周波数成分をゼロにする

filtered_waveform = np.fft.ifft(fft_data)

# FFTを計算して周波数領域で操作（フィルタリング後）
fft_data_filtered = np.fft.fft(filtered_waveform) # FFTを計算
freq_filtered = np.fft.fftfreq(len(fft_data_filtered), d=1.0 / sampling_rate) # 周波数軸を生成
amp_filtered = np.abs(fft_data_filtered / (len(filtered_waveform) / 2))


# オリジナル波形の周波数スペクトルをプロット
plt.subplot(2, 1, 1)
plt.plot(waveform)
plt.title('Original Spectrum')
plt.xlabel('tate')
plt.ylabel('yoko')
plt.xlim(0, 2000) # x軸の範囲を0~2000Hzに限定
# フィルタリング後の波形の周波数スペクトルをプロット
plt.subplot(2, 1, 2)
plt.plot(filtered_waveform)
plt.title('Filtered Spectrum')
plt.xlabel('tate')
plt.ylabel('yoko')
plt.xlim(0, 2000) # x軸の範囲を0~2000Hzに限定



# オリジナル波形の周波数スペクトルをプロット
plt.subplot(2, 1, 1)
plt.plot(freq[1:int(len(waveform) / 2)], amp[1:int(len(waveform) / 2)])
plt.title('Original Spectrum')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.xlim(16800, 17200) # x軸の範囲を0~2000Hzに限定
# フィルタリング後の波形の周波数スペクトルをプロット
plt.subplot(2, 1, 2)
plt.plot(freq_filtered[1:int(len(filtered_waveform) / 2)], amp_filtered[1:int(len(filtered_waveform) / 2)])
plt.title('Filtered Spectrum')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.xlim(16800, 17200) # x軸の範囲を0~2000Hzに限定

# グラフの表示
plt.tight_layout()
plt.show()
実際に適用したいのは0.01秒なので、どんどん時間を短くしてちゃんと波形が消えるか試してみました。
すると0.1秒ではなんか波形が消えてました。
もっと短時間でフーリエ変換したい。
そう考えて調べていると、短時間フーリエ変換(stft)というものを見つけました。
そして今まで音で文字を書いていたグラフのことをスペクトログラムと呼ぶことも知りましたｗ

stftならできそうだと考えてchatGPTに書かせて実験しました。
stftの場合（このコードは未完成！）：
import numpy as np
import matplotlib.pyplot as plt
import librosa

# 元の信号を作成
duration = 0.03 # 音声の長さ（秒）
sample_rate = 22050 # サンプリングレート
before_freq=1000
t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)
signal = np.sin(2 * np.pi * before_freq * t) #+ np.sin(2 * np.pi * 2000 * t)

# STFTを計算
n_fft = len(signal) # フレームサイズを信号の長さと同じに設定
hop_length = int(n_fft / 4) # フレームのオーバーラップを設定（例としてフレームサイズの1/4を使用）
stft_matrix = librosa.stft(signal, n_fft=n_fft, hop_length=hop_length)


# 周波数軸の取得
frequencies = librosa.fft_frequencies(sr=sample_rate)

# 対象周波数のインデックスを特定
target_frequency = 1000
tolerance = 10
target_index = np.argmin(np.abs(frequencies - target_frequency))

# 周波数領域で特定の周波数成分をゼロにする
stft_matrix[target_index - tolerance:target_index + tolerance, :] = 0

# 逆STFTを計算
filtered_signal = librosa.istft(stft_matrix, hop_length=hop_length, length=len(signal))

# 元の信号とフィルタリングされた信号をプロット
plt.figure(figsize=(10, 4))
plt.plot(t, signal, label='Original Signal')
plt.plot(t[:len(filtered_signal)], filtered_signal, label='Filtered Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Signal with Frequency Removal')
plt.legend()
plt.show()
このコードは最初は動いていたんですが、n_fftとかhop_lengthって何？と弄っていたらエラーを吐きました。
しかし、このときのn_fftとかhop_lengthって何だろうと思ったおかげで最後の課題を解決できました。
そのため今から実験用コードを修正する意味ないじゃんと面倒くさくなって放置しました。。

n_fftとかを理解するためにstftを勉強しようと思いました。
ググったり、chatGPTに聞いたり、youtubeの動画を見ました。

もともとstftは短時間フーリエ変換という名前の通り、波を短時間だけ切り取ってフーリエ変換するものだと理解していました。
この波を切り取るときに窓関数を使用しており、ただ切り取るだけではダメなようです。
なんで窓関数を使用しないといけないのか調べると、ただ切り出すだけではその切り出した両端の波が不連続になってしまうからだそうです。
不連続な波をフーリエ変換すると低周波数から高周波数まで幅広く広がってスペクトル漏れが生じてしまう！！！！！

すなわち今まで周波数が全体に広がっていたのは波が不連続だったから。
そしてその不連続な波を連続的にするには窓関数を使えばいい！！
窓関数にはいろいろ種類があるようですが、有名な3つを試してみました。
（もうこのときにはn_fftとhop_lengthを忘れています。。）

hamming window:
ちょっとスペクトル漏れが薄くなった！？
hanning window:
おおおお！！！！
ほとんどスペクトル漏れがない！！！！

blackman window:
おおおおおおおおおおおおおおおおおおおおおおお！！！！！！！！！！！！！！！！！！！！！！！
ついに！！！！！
めっちゃきれいに文字だけを出せたぁあああああああああああああああああああ！
ブラックマンありがとう！！！！！！！！！

この窓関数、たった一行を書くだけでもすごい時間がかかりました。
ほんとにちょっとコードが違うだけで結果が大きく変わってくるのがプログラミングに限らないですけど面白いですよね。

こんな感じでついにきれいに文字を音で表現することに成功しました！！
いやー、すごい時間がかかりましたね。
少しずつ積み重ねてやっと達成できたので達成感がすごかったです！
この方法では文字を画像化して画像を読み込んでいるので、画像をそのまま描画することもできます。
音で表現したい画像を見つけて、ペイントでモノクロで保存して読み込むだけです。
読み込むときに画像によってRGBの戻り値が変化することがあるのでそこだけ注意です。
解像度とかをうまく調整すればもっときれいにできそうですがまあこんな感じです。
あんまりかわいくなくなってしまった・・・


最後にコードを載せておきます。
コメントアウトもそのままだったりとめっちゃ汚いですが、自分が試行錯誤した記録なので残しときます。

調べてみるとCTFにもやっぱりスペクトログラムの問題が出てるみたいですね。

また別の方法として、先に波を作っておいてから画像の白と黒で0と1をかけ合わせて作る方法もあるみたいですね。
確かにこの方法なら波が連続しているのでスペクトル漏れは起きなさそうですね。

結構疲れましたがいろいろ勉強になりました。
友人とchatGPTには本当に感謝です。
夢中になれて楽しかったです！

主なサイトまとめ：
あとリアルタイムでスペクトログラムを見たいときには

import librosa
import numpy as np
import wave
import matplotlib.pyplot as plt
from scipy import signal
from PIL import Image
# サンプリング周波数
sample_rate = 48000
#音の範囲
Freq1=17000
Haba=3000
Freq2=Freq1+Haba

def main():
makemoji("Hello")
image_path = "A.png"
 # 画像を開く
 image = Image.open(image_path)
# 画像のサイズを取得
 width, height = image.size
# 画像のピクセルデータを取得
 pixels = image.load()
w1 = np.empty(0)
f_base=Freq2
dy = 50
 dt = 0.5/dy
#time = np.linspace(0, 1, int(sample_rate*dt), endpoint=False)
 time = np.linspace(0, dt, int(sample_rate*dt), endpoint=False)
#time = np.linspace(0.14, 0.14, int(sample_rate*dt), endpoint=False)
 #time = np.random.uniform(0,np.pi/2, int(sample_rate * dt))
 # 各ピクセルを取り出す
 for x in range(0,width,int(height/dy)):
dw = np.zeros(int(sample_rate * dt))
for y in range(0,height,int(height/dy)):
# ピクセルのRGB値を取得
 #print(pixels[x, y])
 #r, g, b, _ = pixels[x, y]
 r, g, b = pixels[x, y]
# 取得したRGB値を表示（例：x座標、y座標、R値、G値、B値）
 #print(f"({x}, {y}) - R: {r}, G: {g}, B: {b}")
 if r==0 and g==0 and b==0:
num=y/(height/dy)
ff1=f_base-(Haba/dy)*num
sa=dt/dy*num
time = np.linspace(0+sa,dt+sa, int(sample_rate * dt),endpoint=False)
#dw+= makeWave(ff1, time,sa)
 window=np.blackman(len(dw))
dw+= makeWave(ff1, time,0)*window
"""
 for i in range(int(ff1), int(ff1+3000/50), 1):
 sa = dt / (3000) * (num*(3000/50)+i-ff1)
 time = np.linspace(0 + sa, dt + sa, int(sample_rate * dt), endpoint=False)
 dw += makeWave(i, time,0)
 """
 w1 = np.append(w1, dw)
"""
 w2=0
 ##余計な周波数を取り除きたい
 for x in range(0,width,int(height/dy)):
 dw=w1[int(sample_rate * dt*x):int(sample_rate * dt*(x+1))]
 for y in range(0,height,int(height/dy)):
 r, g, b = pixels[x, y]
 if (r==255 and g==255 and b==255):
 num=y/(height/dy)
 ff1=f_base-(3000/dy)*num

 dw+= remove_frequency_stft(dw, sample_rate, ff1, tolerance=0)
 w2 = np.append(w2, dw)
 #"""
 ##
 # 画像を閉じる
 image.close()

f1 = Freq1
f2 = Freq2
fp = np.array([f1, f2]) # 通過域端周波数[Hz]※ベクトル
 fs = np.array([f1 - 1000, f2 + 1000]) # 阻止域端周波数[Hz]※ベクトル
 gpass = 3
 gstop = 40
 #w1 = bandpass(w1, sample_rate, fp, fs, gpass, gstop)

 waveform = np.empty(0)
waveform = np.append(waveform, w1)
#waveform = np.append(waveform, w2)
 #waveform = np.append(waveform, w3)
 """
 a=500
 #print(len(waveform[:a]))
 waveform=remove_frequency_stft(waveform[:a], sample_rate, 18000, tolerance=10)
 #print(len(waveform))
 """
 #正規化
 waveform /= max(abs(waveform))
#16ビットに変換
 waveform = np.array(waveform * 32767, dtype=np.int16)

#plotWave(waveform)

 # STFTを計算
 stft = np.abs(librosa.stft(waveform.astype(float)))
# STFTのグラフを表示
 plt.figure(figsize=(10, 4))
librosa.display.specshow(librosa.amplitude_to_db(stft, ref=np.max), sr=sample_rate, x_axis='time', y_axis='linear')
plt.colorbar(format='%+2.0f dB')
plt.title('STFT')
plt.xlabel('Time')
plt.ylabel('Frequency')
plt.ylim([Freq1-1000, Freq2+1000])
plt.tight_layout()
plt.show()
"""
 N = len(waveform)
 fft_data = np.fft.fft(waveform)
 freq = np.fft.fftfreq(N, d=1 / sample_rate)
 amp = abs(fft_data / (N / 2))
 plt.plot(freq[1:int(N / 2)], amp[1:int(N / 2)])
 plt.show()
 """
 # WAVファイルを作成
 with wave.open('sound.wav', 'w') as file:
file.setnchannels(1) # モノラル
 file.setsampwidth(2) # 16ビット
 file.setframerate(sample_rate)
file.writeframes(waveform.tobytes())

def makedw(f1,f2,dt):
#dt = 0.014
 time = np.linspace(0, dt, int(sample_rate * dt), endpoint=False)
dw = np.linspace(makeWave(f1,time),makeWave(f2,time),int(sample_rate*dt),endpoint=False)
#dw = np.zeros(int(sample_rate * dt))
 #for i in range(int(f1), int(f2), 0.1):
 # dw += makeWave(i, time)
 return dw
def makeWave(frequency, time,b):
return np.sin(2 * np.pi * frequency * time+b)
def bandpass(x, samplerate, fp, fs, gpass, gstop):
fn = samplerate / 2 #ナイキスト周波数
 wp = fp / fn #ナイキスト周波数で通過域端周波数を正規化
 ws = fs / fn #ナイキスト周波数で阻止域端周波数を正規化
 N, Wn = signal.buttord(wp, ws, gpass, gstop) #オーダーとバターワースの正規化周波数を計算
 b, a = signal.butter(N, Wn, "band") #フィルタ伝達関数の分子と分母を計算
 y = signal.filtfilt(b, a, x) #信号に対してフィルタをかける
 return y #フィルタ後の信号を返す
def plotWave(freq):
plt.plot(freq)
plt.show()

def makemoji(text):
from PIL import Image, ImageDraw, ImageFont

# 画像サイズと背景色の設定
 height = 200
 background_color = (255, 255, 255) # 白色

 # テキストの設定
 #text = "Hello, World!"
 font_size = 200
 text_color = (0, 0, 0) # 黒色
 font_path = "C:\Windows\Fonts\yumin.ttf" # フォントファイルのパス

 # フォントの読み込み
 font = ImageFont.truetype(font_path, font_size)
#width, _ = font.getsize(text)
 width, _ = font.getbbox(text)[2:]
# 画像を作成
 image = Image.new("RGB", (width, height), background_color)
draw = ImageDraw.Draw(image)

# テキストを画像中央に配置
 #text_width, text_height = draw.textsize(text, font=font)
 lux,luy,rdx,rdy= draw.textbbox((0, 0), text, font=font)
x = (width - (rdx-lux)) // 2
 y = (height - (rdy-luy)) // 2

 # テキストを描画
 draw.text((x, y), text, font=font, fill=text_color)

# 画像を保存
 image.save("A.png")

def remove_frequency_stft(signal, sample_rate, target_frequency, tolerance):
# STFTを計算
 stft_matrix = librosa.stft(signal)

# 周波数軸の取得
 frequencies = librosa.fft_frequencies(sr=sample_rate)

# 対象周波数のインデックスを特定
 target_index = np.argmin(np.abs(frequencies - target_frequency))

# 範囲の制限
 start_index = max(0, target_index - tolerance)
end_index = min(target_index + tolerance, stft_matrix.shape[0] - 1)

# 周波数領域で特定の周波数成分をゼロにする
 stft_matrix[start_index:end_index + 1, :] = 0

 # 逆STFTを計算
 filtered_signal = librosa.istft(stft_matrix)

return filtered_signal




if __name__ == "__main__":
main()