फूरियर स्पेस फ़िल्टरिंग

Question

मेरे पास एक वास्तविक वेक्टर टाइम सीरीज़ एक्स लम्बाई टी है और लम्बाई टी < < का फ़िल्टर एच है। टी एच चौकोर स्पेस, असली और सममित में एक फ़िल्टर है। यह लगभग 1/एफ है।

मैं y प्राप्त करने के लिए x के साथ x फ़िल्टर करना चाहता हूं।

मान लीजिए कि टी == टी और एफएफटी की लंबाई टी स्मृति में फिट हो सकती है (इनमें से कोई भी सत्य नहीं है)। अजगर में मेरी फ़िल्टर किए गए एक्स पाने के लिए, मुझे क्या करना होगा:

import numpy as np 
from scipy.signal import fft, ifft 

y = np.real(np.ifft(np.fft(x) * h))) 

की स्थिति के बाद से पकड़ नहीं है, मैं करने की कोशिश की निम्नलिखित हैक:

1. एक गद्दी आकार पी < टी/2 का चयन करें, का चयन एक ब्लॉक आकार बी ऐसी है कि बी + 2P एक अच्छा FFT आकार है
पट्टी प्रक्षेप के माध्यम से
2. स्केल ज आकार बी के होने के लिए + 2 पी> टी (h_scaled)
3. y = []; लूप: लंबाई बी के
   • लें ब्लॉक + x से 2P (बुलाया x_b)
   • प्रदर्शन करना y_b = IFFT (fft (x_b) * h_scaled)
   • ड्रॉप गद्दी पी y_b के दोनों तरफ और y के साथ जोड़ से
   • का चयन अगले x_b पी द्वारा पिछले के साथ अतिव्यापी

यह एक अच्छी रणनीति है? मैं पैडिंग पी को एक अच्छे तरीके से कैसे चुनूं? ऐसा करने का सही तरीका क्या है? मुझे बहुत सिग्नल प्रोसेसिंग नहीं पता है। यह सीखने का एक अच्छा मौका है।

मैं चीजों को गति देने के लिए क्यूएफएफटी का उपयोग कर रहा हूं, इसलिए यह बहुत अच्छा होगा अगर ऑपरेशन का बड़ा हिस्सा एफएफटी है। वास्तविक समस्या 3 डी है। इसके अलावा, मैं एक आकस्मिक फ़िल्टर से कलाकृतियों के बारे में चिंतित नहीं हूं।

धन्यवाद, पॉल।

Answer 1

आप सही रास्ते पर हैं। तकनीक को overlap-save processing कहा जाता है। t इतना छोटा है कि उस लंबाई के एफएफटी स्मृति में फिट हैं? यदि ऐसा है, तो आप अपना ब्लॉक आकार B चुन सकते हैं जैसे B > 2*min(length(x),length(h)) और तेज़ रूपांतरण के लिए बनाता है। फिर जब आप प्रक्रिया करते हैं, तो आप दोनों सिरों से छोड़ने के बजाय y_b के पहले भाग को छोड़ देते हैं।

यह देखने के लिए कि आप पहली छमाही क्यों छोड़ते हैं, याद रखें कि स्पेक्ट्रल गुणा समय डोमेन में परिपत्र रूपांतरण के समान है। शून्य-गद्दीदार h के साथ कनवॉलिंग परिणाम के पहले भाग में अजीब गड़बड़ाने वाले ट्रांजिस्टर बनाता है, लेकिन दूसरे छमाही तक सभी ट्रांजिस्टर चले गए हैं क्योंकि x में परिपत्र रैप बिंदु h के शून्य भाग के साथ गठबंधन है। "Theory and Application of Digital Signal Processing" by Lawrence Rabiner and Bernard Gold में चित्रों के साथ, इसकी एक अच्छी व्याख्या है।

यह महत्वपूर्ण है कि आपका टाइम डोमेन फ़िल्टर कम से कम एक छोर पर 0 पर टेंडर करें ताकि आपको रिंगिंग कलाकृतियां न मिलें। आप उल्लेख करते हैं कि h आवृत्ति डोमेन में वास्तविक है, जिसका अर्थ है कि इसमें सभी 0 चरण हैं। आम तौर पर, ऐसा संकेत केवल एक गोलाकार फैशन में ही होगा, और जब फ़िल्टर के रूप में उपयोग किया जाता है तो आवृत्ति बैंड के माध्यम से विकृति पैदा होगी। एक उचित फ़िल्टर बनाने का एक आसान तरीका है इसे आवृत्ति डोमेन में 0 चरण, उलटा परिवर्तन, और घुमाएं।उदाहरण के लिए:

def OneOverF(N): 
    import numpy as np 
    N2 = N/2; #N has to be even! 
    x = np.hstack((1, np.arange(1, N2+1), np.arange(N2-1, 0, -1))) 
    hf = 1/(2*np.pi*x/N2) 
    ht = np.real(np.fft.ifft(hf)) # discard tiny imag part from numerical error 
    htrot = np.roll(ht, N2) 
    htwin = htrot * np.hamming(N) 
    return ht, htrot, htwin 

(मैं पाइथन के लिए काफी नया हूं, कृपया मुझे बताएं कि कोड को बेहतर तरीका है या नहीं)।

ht, शीर्ष पर, लहर के बहुत सारे है:

आप ht, htrot, और htwin की आवृत्ति प्रतिक्रियाओं की तुलना करते हैं, तो आप देखते निम्नलिखित (x- अक्ष pi अप करने के लिए सामान्यीकृत आवृत्ति है) । यह किनारे पर असंतोष के कारण है। htrot, बीच में, बेहतर है, लेकिन अभी भी लहर है। htwin थोड़ा अधिक आवृत्ति पर फ़्लैटनिंग की कीमत पर अच्छा और चिकना है। ध्यान दें कि आप एन

के लिए एक बड़ा मूल्य का उपयोग कर सीधी रेखा अनुभाग की लंबाई बढ़ा सकते हैं, मैंने असंतोष के मुद्दे के बारे में लिखा था, और यदि आप अधिक विवरण देखना चाहते हैं तो another SO question में एक मैटलैब/ऑक्टेव उदाहरण भी लिखा है।