The easy way is to use Zonal statistics as a table to extract the mean value for all buffers, then join the table to your shapefile. However, there might be a problem if your buffers overlap : in this case, it gets more difficult because you need to iterate on each polygon (e.g. in model builder).
Alternatively, you can use focal statistics with a circular neighborhood of the size of your buffer to create a mean NDVI raster, then extract the value (extract multiple value to point) at the location of each observation point. This should give you the same result on average, with small differences due to the rasterization.
When you perform the analysis, make sure that your environment settings are set to the same pixel size as your image and that you snap your image (Zonal stat will actually rasterize your shapefile internally)
you can use the field calculator with the Python parser
sum([!field1!, !field2!, !field3!, !field4!])/4
EDIT: for accounting for null values
sum([a for a in [!field1!, !field2!, !field3!, !field4!] if a is not None )/4
not that, in this case, it works as if you assume that null values are ZERO as in your comment. Alternatively, you can IGNORE null values in the mean (so you divide with the number of valid field values.
sum([a for a in [0, !field1!, !field2!, !field3!, !field4!] if a is not None )/max([1,len([a for a in [!field1!, !field2!, !field3!, !field4!] if a is not None] )
Here I used a little trick to avoid errors when all fields are Null (will return 0). This is a weak error handling, but in one line.
Best Answer
Do you have access to the Spatial Analyst extension?
ArcToolbox > Spatial Analyst Tools > Neighborhood > Focal Statistics