CNN 特徴マップの可視化

CNNで畳み込み処理後・プーリング後の特徴マップの可視化を行いたいのですが，
調べてみると，VGG16をロードしてから特徴マップを可視化しているサイトばかりで，
いまいち具体的なコードの書き方が分かりません。

import matplotlib.pyplot as plt
import os
import cv2
import random
import numpy as np
import tensorflow as tf
from sklearn.model_selection import train_test_split
import time
import feature_visual
import filter_visual

DATADIR = "aa"
CATEGORIES = ["z2","x2"]
IMG_SIZE = 28
training_data = []
def create_training_data():
    for class_num, category in enumerate(CATEGORIES):
        path = os.path.join(DATADIR, category)
        for image_name in os.listdir(path):
            try:
                img_array = cv2.imread(os.path.join(path, image_name), cv2.IMREAD_GRAYSCALE) 
                img_resize_array = cv2.resize(img_array, (IMG_SIZE, IMG_SIZE))  
                training_data.append([img_resize_array, class_num])  
                
            
            except Exception as e:
                pass
create_training_data()
random.shuffle(training_data)  
X_train = [] 
y_train = [] 

for feature, label in training_data:
    X_train.append(feature)
    y_train.append(label)

X= np.array(X_train)
y= np.array(y_train)

x_train, x_test, y_train, y_test = train_test_split(X, y, test_size=0.20)

x_train, x_test = x_train / 255.0, x_test / 255.0　　　　　　　　　　　　　# データセットの準備



model = tf.keras.models.Sequential()
model.add( tf.keras.layers.Reshape((28, 28, 1), input_shape=(28, 28)) ) 
model.add( tf.keras.layers.Conv2D(16, (3, 3), activation='relu') )  
model.add( tf.keras.layers.Conv2D(32, (3, 3), activation='relu') ) 
model.add( tf.keras.layers.MaxPooling2D(pool_size=(2,2)) ) 
model.add( tf.keras.layers.Dropout(0.5) )           
model.add( tf.keras.layers.Flatten() )                                
model.add( tf.keras.layers.Dense(128, activation='relu') )
model.add( tf.keras.layers.Dropout(0.25) )
model.add( tf.keras.layers.Dense(2, activation='softmax') )
    
    

model.compile(optimizer='adam',loss='sparse_categorical_crossentropy',metrics=['accuracy'])　　　　　　　　# モデルの学習
model.fit(x_train, y_train, batch_size=256, epochs=40)


loss, accuracy1 = model.evaluate(x_test,  y_test, verbose=2)　　　　　　　　　　# モデルの検証

コード