How do Anti-Aliasing Filters filter the signal? I have programmed a anti-aliasing filter that uses the fourier-transformation and makes a perfect cutoff on my desired cutoff frequency for downsampling without aliasing.
I have read about other types of filters for anti-aliasing. Here i have a nice picture of different filters from this site:

The dotted line seems for me the perfect or ideal filter, but in practical it is not possible to filter a signal perfectly right? And therefore we can use different filters with different orders.
Does this graph mean, the different filters do not filter completely on the perfect cutoff frequency and always contains a little aliasing?

The programmed FFT-Filter is used as followed in python:

    def lowpass(x, binmax):
        N = len(x)
        return np.fft.irfft(np.fft.rfft(x, axis=0)[:binmax], N, axis=0)

x48hr is a Signal sampled with 48kHz and do a cutoff at 4kHz
lr_signal = lowpass(x48hr, binmax=len(x48hr)*4//48) #hr

I use this filter for generating low-resolution samples to train an convolutional neuronal network to learn the mapping between lr and hr to make a hr signal reconstruction.

In general I am actually trying to understand how much aliasing my measurements have and how I can use my measurements to study different types of filters and the aliasing in them. I want to increase the resolution of my signal, therefore I need aliasing in my low-resolution signals. I have asked another question regarding my problem here if you want to know more about it.

Yes that is correct. The width of the transition band is inversely proportional to the length (or memory) of the filter which means to say that you would need an infinite amount of time to achieve your perfect filter. Therefore for practical reasons you decide how much aliasing would be tolerable (similar in many ways to deciding how many decimal places you need to include for a number describing the length of something) and then design your filter accordingly. The tighter the rejection the longer will be the filter delay and higher will be the overall complexity.


What is the use of perception calculation in Speaker Identification? I'm currently trying to find audio features to train a neural network that distinguishes people, and I came across MFCC which is used a lot for speech verification. As far as I understand the frequency cepstral coeffients are measured on the mel scale which is logarithmically and not linear because of the human hearing being logarithmical.

So human perception says that the difference between $100\textrm{ Hz}$ and $200\textrm{ Hz}$ is higher than between $10100\textrm{ Hz}$ and $10200\textrm{ Hz}$.

How does this this "less perceived difference" help me to distinguish people?

Well, I don't know whether it'll actually help you – you just said it would!

Now, in any case, using an algorithm to extract features from a signal that mimics or resembles human perception should inherently give those feature vectors a higher "mathematical" resemblance when a human would find the original signals similar, too.

This is obviously what you want in speech recognition – after all, the computer is supposed to do something that is very human.

Now, the question is whether that helps in identifying speakers.

In my first, intuitive, reaction to that, I'd say, no, applying what is used for speech recognition to speaker identification might be a bad idea: A good speech recognizer should be as invariant as possible to the speaker – and thus, "lose" the information of who spoke somewhere on the way.

In reality, something less than that extreme is true: You still want to analyze speech (as opposed, to let's say, interpret engine sounds) – but under a different aspect. So, using the same  basic steps to extract feature vectors sounds pretty clever. You might want to tweak parameters; for example, though not overly important for speech recognition, higher-frequency sound might be helpful in telling high-pitched female teens from smokey-voiced old men.


Error in power spectral density computation I've written a program to compute the power spectral density of a signal. The problem I'm having is that the computed spectrum is quite a bit more noisy than it should be (comparing to the spectrum computed by Matlab), especially at low frequencies.

http://i.imgur.com/fGOGT.png

The red in this plot shows the spectrum computed by Matlab's cpsd command (0 overlap, NFFT=window=8192, 200Hz signal), and the green shows my own. There's a constant offset (in the log domain) between them, but that's okay for now. What I'm more concerned about is that my implementation is noisier, especially in the 0-3Hz range. The algorithm I've implemented is

Zero out a vector of length NFFT/2+1, for averaging.
Until there are no more samples:

Read NFFT samples (0-fill the remainder if there are less than NFFT)
Compute the 1-dimensional real NFFT
Compute the magnitude of each complex result, and divide by NFFT
Add this to the corresponding position in the averaging vector

Divide each sum by the twice the sampling frequency times the number of windows read

Sizing a window different than NFFT, and overlap, aren't included in the algorithm, but when compared to the Matlab output for window=NFFT and no overlap it shouldn't matter.

Can anyone speak as to why this results in low frequency instability?

Matlab splits your signal into 8 overlapping segments (by default). The difference I see between Matlab and your algorithm is that Matlab also multiplies each segment by a Hamming window. You should get pretty close results if you do that.


Why do I see a frequency shift of the peak between the S11 and impedance response? I have a resonator I am trying to characterize, and the S11-Plot gives me a slightly different center (resonant) frequency than that of the Z-Plot.
These measurements were taken simultaneously and using the same VNA. Yet the Z (magnitude and phase) are shifted to the left by something like 5% (Way more than VNA resolution).
Is there a physical reason, or a calibration error reason, or any other for this?

[This is more of an electronics related question and perhaps should be migrated to Electronics.SE?]
This can occur with a series parasitic inductance or capacitance in the resonator; the case with a series inductance is demonstrated below.
The Smith Chart visual provides immediate insight of a shunt L-C resonator (tank circuit) to ground. In the first plot where we have no parasitic inductance in series with an exaggerated parasitic resistance in the resonator (reducing the input impedance at resonance when the reactive inductance cancels the reactive capacitance). At resonance $X_{Lp}$ cancels $X_{Cp}$ and we would be at the large dot shown (with just the shunt resistance due to finite Q remaining, in parallel with the 50 ohm load and thus reducing the input impedance along the resistance axis). As we move off resonance lower in frequency the tank circuit is inductive and ultimately a short to ground at DC, following the shunt inductance trajectory as shown by the path marked $L_p$ indicating a parallel inductor to ground, as we sweep the test frequency toward DC terminating on the left origin in the Smith Chart consistent with a short circuit. As we go higher in frequency the tank circuit will be capacitive and ultimately a short to ground at $f=\infty$, following the shunt capacitor trajectory as shown by the path marked $C_p$ as we sweep the frequency higher.   In this case for all other frequencies other than at resonance the net impedance of the tank would smaller with a maximum tank impedance at resonance, and the resonance peak would match the S11 null.

Observe if the tank circuit had a series inductance effectively prior to the shunt L and C to ground. This inductor in series would follow the trajectory marked by $L_s$ (and would ultimately end at the right origin where $Z=\infty$, terminating at the large dot shown at resonance with $X_{Lp} = -X_{Cp}$. As we move off resonance upward in frequency, the resonator adds a shunt capacitance that in this case brings us closer to the center origin (where Z=50 ohms, perfect match), providing a lower S11 as it maximizes the net impedance of the resonator's connection to ground; essentially the resonator is tuning out the series inductance as it becomes capacitive. Not detailed here is that series inductance would also increase as frequency increases but we can see how if the capacitance relative to the inductance is significantly larger that the trajectory would move down faster than up and to the right resulting in the minimum as the OP is seeing. Further the inductance itself would have a point where it resonates out and becomes capacitive, so we could also easily envision the possible case of a net circular path of the series inductance and shunt capacitance of the resonator passing right through the center of the Smith chart (or closer to it).

Below is a basic simulation of this compliments of https://www.will-kelsey.com/smith_chart/
At resonance the Smith chart would be at DP5 (where the red line intersects the black line to the upper left of center. As we sweep frequency, the trajectory follows the red line where we see a minimum impedance is reached which is not on the resonance of the shunt L and C.


Sparse CWT matrix-shape data lossless compression Is there a specialized lossless compressor for sparse matrices, i.e. large fraction (>40%) of values == 0? LZ4 worked well on 1D signals, but now whether I feed a sparse 2D or an all-zeros array makes no difference, still ~2x compression.
Slightly lossy compression (<5% MAE) also suffices. If relevant, data is float32 not uint8 (.mat/.npy), and it's not actually an image but a 2D array resulting from absolute value of synchrosqueezed CWT of a signal. float16 might be an option, I've not decided, but ideally float32 is preserved in compression. Data ultimately fed to a neural net.

Start simple: that calls for any Run-Length coding (RLE).
Also, floating point makes very little sense for image data, so the first step that will probably save you 50% of space is conversion to unsigned fixed point (multiply every value with 2¹⁶, then divide through the largest value of the original image, round to nearest integer, keep the result as uint16).
5% mean absolute error means: you don't need floating point; rounding small values to zero is absolutely¹ acceptable to you. So, drop the float idea and convert to integer, as described above.
If for some (somewhat untypical for image data) reason you want to keep it floating point, convert it to a format where you store the mantissa and exponent bits separately and the latter at all only for non-zero mantissa values. Since pixel data typically only has a positive sign, you don't even need to store the sign bit.
After reducing the number of exponents you've stored, maybe any further compression isn't necessary to achieve your goal (don't overengineer!).
Then, look at your data: It's also quite likely that your exponents don't need the full 2⁻¹²⁶ to 2¹²⁷ range, so dropping bits on your exponents makes a lot of sense, too. Make sure you're not totally killing it with the number of mantissa bits, too. Wild guess: half as many as in 32 bit IEEE754 will do.
If after doing these thing it still is necessary to compress further, RLE or Lempel-Ziv would be appropriate solutions.
(LZ4 is kind of exotic; the more established LZMA algorithm as e.g. in XZ should work sufficiently well, and it really doesn't sound like you're CPU bound here...)
In all seriousness, though, sounds like a job PNG, high-quality JPEG, or if you have access to that, JPEG-XL, could do out of the box, so I'd simply look into converting my data to a format that makes sense to readily available lossless image codecs instead of trying to roll my own.

¹ hurr durr, a pun!


Extracting Peak Frequencies Using FFT vs. Time Domain Peak Finding I am trying to compare two processing algorithms that I believe should be producing very similar results. I am processing acceleration signals over the course of the day and we expect to see a decrease in dominant frequency midday and an increase in the frequency at night. The frequency should be ~0.3 Hz and increase and decrease (by ~0.05) depending on the time of day.

My sampling frequency is 64 Hz and I am trying to run this analysis on the data we already have, so unfortunately we cannot run this experiment again with a lower sampling frequency.

First method using FFT:
1) I apply a bandpass filter on about a ~4 hour chunks of the signal allowing 0.1 Hz to 0.5 Hz. This leaves me with a signal of about 900,000 samples.
2) I then run a cross-correlation on this signal
3) I'm using MATLAB, and run the pwelch function on the cross-correlated signal. I run this with a gaussian window of length 3,500. 50% overlap.
4) I then attempt to smooth out the FFT spectrum with a Savitzky-Golay filter.
5) I fit the peak with the peakfit function
6) I extract the peak frequency within the 0.2 - 0.4 Hz range and call that the dominant frequency for that 4 hour chunk of time.

Second method using time-domain processing:
1) I run the same bandpass filter as above on the signal, same time chunk. I do not run a cross correlation.
2) I find all the peaks in the filtered time domain signal using the findpeaks function
3) I then find the difference between the successive peaks and create a histogram, which looks like a normal distribution centered on a certain period
4) I calculate the mean period of this distribution
5) I take the inverse of the mean period and call that the dominant frequency for that 4-hour chunk

I then run these algorithms on 3 weeks of a continuous acceleration signal and obtain very similar frequencies using both methods. However, about 25% of the peak extractions on the FFT seem wildly off and noisy. For example from point-to-point, meaning from one four-hour chunk to the next four-hour chunk, the frequency change is 50-100% different. On the other hand our second method (the time-domain method) never shows this kind of noise or variance and the transition to different time periods is much smoother and what we expect.

I feel like I can keep playing around with the window lengths/different windowing methods with the FFT,  but I had a few questions first.

1) I haven't found any thing in the literature with a similar description to what I am doing with my time-domain method. Am I doing something completely unnecessary here? Is our method valid or just arbitrary? Should I just focus on doing FFT processing? This method makes sense to me, but I cannot back up this method with any relevant literature.

2) Because the time series of extracted frequencies using my FFT method is somewhat noisy, I am leaning towards moving forward with my time-domain method, making my first question more pressing to justify because we would like to publish these results.

3) I'd be happy to hear thoughts on either of these algorithms

Thanks!

You basically have a slowly changing Sine signal where the parameter which changes is its frequency.

The right approach in my opinion would be to use Frequency Modulation processing of the signal.
Since you are dealing with really low frequencies yet your observation period is long you should be OK and have good results.


Can I segment characters from this image with CSV? I am trying to segment my characters on this image with using column sum vector.

This method actually works on the paper I've read about license plate recognition.

it looks useful for my images too but since my threshold values are not very good I am not really sure if I can use CSV for character segmentation

so my question is;

1) Is CSV proper method for my application ?
2) if so what kind of techniques I can try to segment characters ?
3) if not what can be more efficient and robust way to segment characters from some kind of images ?

thank you

Your results (above) looks well for character segmentation. But if you have some problems with calculation of thresholds for total image, you can try slightly different approach. "Column sum vector" method works well for grey images too, but not only for black/white images. So there is no need to calculate threshold for total image, but you must choose some value for segmentation in vector. May be it will slightly simpler to do segmentation in this case.

You can use "line sum vector" for vertical segmentation   too.


Confusion about hough transformation for ellipse I saw the following sentence from Feature Extraction & Image Processing, 2nd Edition

  $$
x = a_0 + a_x\cos(\theta)+b_x\sin(\theta)
$$ $$
y = b_0 + a_y\cos(\theta)+b_y\sin(\theta) $$

  This equation corresponds to the polar representation of an ellipse.

I don't know this meaning. If you need extra information, I'll add it.

It's a little mis-leading and confusing.

The usual way to describe an ellipse is using cartesian ($x,y$) coordinates.  Another way to represent an ellipse is using polar coordinates ($r,\theta$) (where $r$ is the distance (radius) from the origin, and $\theta$ is the angle.

What the equation is saying is that the cartesian coordinates of the ellipse are parameterized by $\theta$ in the two equations. That is, as you change $\theta$ from 0 to $2\pi$, the values for $x$ and $y$ will lie on the ellipse, for constants ($a_0,b_0,a_x,b_x,a_y,b_y$).


Will noise cancellation with playstation eye improve my spectral analysis? Im doing some beginners spectral analysis for a homemade reactive RGB LED project and was thinking of ways to improve the accuracy of the audio analysis. Im using a Playstation eye with a Raspberry Pi Zero W (limited performance) and ive read several articles to try to better my understanding. According to this article, the playstation eye has 4 channels, where the 2nd and 3rd are reversed for noise cancellation.
I assume that the 2nd channel is a reverse of the 1st and the 3rd of the 4th. Am i correct to assume that the way to go about implementing these extra channels would be to average the two sums of 1+2 and 3+4 ie. voltage = (c1+c2+c3+c4)/2?

Below is a script i use for testing, only the basics. Note that i use a sampling rate of 20khz instead of 44.1khz due to the limited cpu resources of the pi zero. I also use int16 as wiki states that the ps eye uses 16bit ints. Am i wrong not to use the default float?

The reason i want to improve it is that i havent recorded any frequencies over 7000, which seems a bit odd to me given the vast variety of music ive tested with. Im wondering if the playstation eye is not suitable for such precise analysis, or if my methods are wrong. Is it possible to improve this with noise cancellation, or will it worsen my results? Should i rather just use a single channel or sum channel 1+4?

#!/usr/bin/env python3

import sounddevice as sd
import numpy as np

RATE = 20000
lowest = 200
highest = 0
h_m = 0

def print_sound(indata, frames, time, status):
    global lowest, highest, h_m
    if status:
        print(status)
    if frames == 0:
        return
    data = np.sum(indata, axis=1)
   # data = data * np.hanning(len(data)) # smooth the FFT by windowing data
    fft = abs(np.fft.rfft(data)) # calculate magnitudes from positive fft coefficients
    freq = np.fft.rfftfreq(frames, 1.0/RATE) # calculate positive frequency bins
    max = np.max(fft) # highest magnitude
    freqPeak = freq[np.where(fft==max)[0][0]] # find frequency with highest magnitude
    if freqPeak < lowest:
        lowest = freqPeak
        print("lowest: %d"%freqPeak)
    if freqPeak > highest:
        highest = freqPeak
        print("highest: %d"%freqPeak)
    if max > h_m:
        h_m = max
        print(max)

with sd.InputStream(dtype='int16', channels=4, samplerate=RATE, callback=print_sound):
    while True:
        response = input()
        if response in ('', 'q', 'Q'):
            break

As I read your prompt, it seems that the real question you’re asking is why you can’t record frequency peaks greater than 7kHz. The second question seems to be whether you can improve this performance through noise cancellation techniques.

If that’s the case, and if you’re right about the 2nd and 3rd channels being reversed, then I think your summation may introduce frequency filtering effects in your output and not necessarily true noise cancellation. This can either work for or against you. But this would depend on the orientation between the microphone and speakers, since in your link it appears that can affect the delays between the channels. Either way, I’d go with the single channel to simplify things.

And if you’ve only been recording music and not individual tones at greater than 7kHz, then it also makes sense that you’re not grabbing high frequency peaks since audio energy tends to be concentrated at lower frequencies. There’s a reason why telephones and AM radio cut off at around 5kHz. You might also be attenuated because of sinc roll off during the sampling process, although I suspect that’s not the big issue here.

Have you tried recording the entire spectrum, instead of just the peaks? And have you tried comparing that against the FFTs of the original audio files?


Applying a time delay with Z-Transform I am trying to take a signal $F(t)$ that has been sampled at some time DelT, I then wish to pass this signal through a channel $H(s)$. To do this I am sampling my signal $H(s)$ at DelT time intervals then finding the Z-transform $H(z)$. How would I go about applying this in matlab?

Currently I am using c2d to convert s-domain to z-transform for my channel but am unsure how I can go about the convolution with my sampled signal which is just an array of numbers. I have tried to take the fft of sampled $F(t)$ then use freqz of $H(z)$. To begin with I am just trying to apply a simple time delay $e^{-st}$ to the signal.

There's a lot of confusion in your question. You refer to "a channel $H(s)$", then to "my signal $H(s)$", so for the sake of clarity, let's properly define:

Input signal - continuous time $f(t)$, impulse sampled $f[n]=f(n\Delta T)$, Fourier transform of continuous time signal $F(s)$, Z-transform of impulse sampled signal $F(z)=F(e^{s\Delta T})$.
Continuous time channel - impulse response $h(t)$, transfer function (Laplace transform of impulse response) $H(s)$.
Discrete time channel - discrete impulse response $h[n]$, transfer function (Z-transform of discrete impulse response) $F(z)$.

Note the definitions above follow a convention:

time-domain signals/functions are lower case ($f, h$); frequency-domain signals/functions are upper case ($F, H$).
arguments of continuous signals/functions are enclosed in parentheses $(t)$; arguments of discrete signals/functions are enclosed in brackets $[n]$.

Then, you write:

  (...) z-transform for my channel but am unsure how I can go about the convolution with my sampled signal which is just an array of numbers

It seems you do not understand the basics of LTI systems. Note that $H(s)$ and $H(z)$ are frequency-domain and $f(t)$ and $f[n]$ are time-domain.

Convolution makes sense when the two convolved functions are in the same exact domain (continuous/discrete time/frequency), although for the usual purposes the two convolved functions should both be in the continuous/discrete time domain.

If you're dealing with a discrete-time signal ("signal which is just an array of numbers"), then you might want to convolve $f[n]$ with $h[n]$. If you already have $H(z)$, then an inverse Z-transform on $H(z)$ would produce $h[n]$.


Python vs Matlab? Which one better for image processing? I want to learn practical image processing from scratch/basics in easy way. Which software will be helpful in this regard? Matlab or python?

[EDIT-20200516: see at the bottom fro comparison with...cars] For learning from scratch, I would not suggest a programming language alone, but instead the couple "teaching materials" (book, lecture notes) + "exercices with a specific  programming language". So if you find a book that you like on "Python for computer vision with exercices" or "Image processing theory and practice with Matlab, that could be interesting starting points.
Also, your programming tastes and skills may evolve, and learning a first programming language helps you learning a second one in general. But laziness sometimes drives you to sticking to your first language, reusing old  librairies.
Last, the purpose is important. In my case, I mostly engineer algorithms as prototypes and proofs of concepts, that can stay as them, or are turned into "solid programs" by people that are better at, and like better, programming with the rules-of-art in lower levels languages, depending on the target.

To that respect, Matlab is great at designing and fine tuning algorithms, possesses a lot of documentation and help that you can follow step-by-step, and enjoys a long list of contributed toolboxes, esp. at MatlabCentral FileExchange.  When the workflow is set, if you care of speed, efficiency, etc., it is time to pass the algorithmic prototyping over to real programmers (C++, or lower level, which I can't do).

Python now has a large community, and has developed toolsets like Scikit-Image, and there is a tutorial for instance at Scikit-image: image processing.
It is also interfaced with OpenCV.

Globally, as long as you grow solid image processing skills, I would think what mostly differ between Matlab and Python are the cost and the trendiness.  On my side, I would switch to Python for machine learning and data science, but I will stick to Matlab for most of my signal processing and image analysis works for a while.

You can also use (free) open source software, with contributed toolboxes and plugins, that often benefit from external publications. I would suggest:

ImageJ: Image Processing and Analysis in Java
Icy: Open Source Image Processing Software
OpenCV: Open Source Computer Vision Library
SciLab: Image Processing & Computer Vision
Julia: Image Processing

For comparisons:

From MathWorks (biased): MATLAB vs. Python: Top Reasons to Choose MATLAB
For RealPython (biased): MATLAB vs Python: Why and How to Make the Switch
Toby Driscoll (more balanced): Matlab vs. Julia vs. Python (compared to BMW sedan, Ford pickup and Testa)

  MATLAB [...] It’s probably still the easiest to learn for basic
  numerical tasks. Meticulous documentation and decades of contributed
  learning tools definitely matter.

  MATLAB is the BMW sedan of the scientific computing world. It’s
  expensive [...]  rock-solid [...] attracts a disproportionate amount
  of hate.

  Python is a Ford pickup. It’s ubiquitous and beloved by many (in the
  USA). It can do everything you want, and it’s built to do some things
  that other vehicles can’t. Chances are you’re going to want to borrow
  one now and then. But it doesn’t offer a great pure driving
  experience.

  Julia is a Tesla. It’s built with an audacious goal of changing the
  future, and it might. It may also become just a footnote. But in the
  meantime you’ll get where you are going in style, and with power to
  spare.


Discretized Beltrami energy In the paper "An Approach Based on Hybrid Genetic Algorithm Applied to
Image Denoising Problem", I encounter this fitness function that by minimizing it, a denoised image of a noisy image is obtained.

It is called Discretized Beltrami energy. According to the paper, $I$ is the image being evaluated, $I_0$ is noisy image, $\nabla I$ is total variation (TV) regularizing term, $\beta$ and $\lambda$ are balancing parameters. $\Omega$ is the set of all points of the image.

What is $\nabla I$? Determinant of $\nabla I$ is a constant, why is there a sum over $\Omega$?

In this particular case, $\nabla I$ is the gradient because the image $I$ is a scalar field. The summation over $\Omega$ would make more sense if $I$ was a vector field. Since the image pixels are samples over fixed locations in space, the $\Omega$ can be dropped.

To understand why this is used as a fitness function, you need to understand Total Variation Minimisation. The assumption there is that a signal has some broad structure that is disturbed by noise. Noise is seen here as large deviations around a central value. But this view would include edges which do not belong to the "noise" part of an image but to what gives it "structure". So, what you are trying to do with total variation regularisation is to estimate the total length of the signal "curve" which largely depends on its structure and then try to minimise the differences between the signal values in areas that it does not "cost" much to do so. If the denoising algorithm tried to suppress an image edge, that would mean a huge change in the length of the curve. But if the denoising algorithm tried to suppress some spurious spike (as low pass filtering would do), then that would not present a large error.

What the authors are trying to do is to use the outputs of three denoising algorithms to descend to the lowest error of this fitness function as quickly as possible by selecting the "best" areas from each denoised image and blend them together with low pass filtering.

Hope this helps.


Remove stripes from image using FFT I try to remove from an image all stripes with slope between 110° and 120° (Fig. 1a) and my first idea was to filter my image in the frequency domain.

As far as I remembered, the spectrum of those stripes should occupy a pie segment between 20-30° (i.e. the original slope - 90°), but to be sure I prepared an artificial black-and-white image with similar pattern (Fig. 2a). It seemed my speculations were correct (Fig. 2b).

Now, I tried to create a mask and simply "cut off" part of the spectrum to supress high frequencies associated with 10-120° stripes.

Only after several trials I came across a mask (Fig. 3a) which satisfies my needs to some degree (Fig. 4). However, I can't figure the relation between this empirically-obtained mask and the spectrum.

Could someone give me some advice where are flaws in my code or in my reasoning?

Here's the MATLAB code I've used for filtering:

f = imread('img/crystal-grain_sample.jpg');
figure, imshow(f); title('f');

m = max(size(f));
P = 2^nextpow2(2*m);
PQ = [P, P];

Fp = fft2(f, PQ(1), PQ(2));

Fc = fftshift(Fp);
Fc_abs = abs(Fc);
Fs_log = log(1+Fc_abs);

%figure; imshow(Fc_abs, []); title('Fc (abs)');
%figure; imshow(Fs_log, []); title('Fs (log)');

Hp = im2double(imread('spec-mask.png'));
Gp = Hp .* Fp;
gp = real(ifft2(Gp));
gpc1 = gp(1:size(f,1), 1:size(f,2));
figure; imshow(gpc1, []); title('g');

And here's the original image:

The problem was I've masked the shifted spectrum. After I've used ifftshift on my mask, I obtained the desired effect.

PS. geometrikal - thanks for putting me on the right track!


Is spectral leakage normal for FFT? I record a sample of 2000 data points, with a sampling rate of 40khz, so freq. resolution is 20hz. I multiplied the chunk of data with a Hanning window prior to doing the FFT.

I recorded a 1khz test tone, and the result of FFT showed a peak at 1khz, as it should. But the magnitudes at 980hz and 1020hz is also significantly higher than the rest, although not as high as the magnitude at 1khz (about 50%).

Is this really bad? Did I apply the Hanning window wrong or is this unavoidable? Will using finer frequency resolution help?

For the Hanning window, the zero crossings occur at a much wider bandwidth (trading for lower side-lobe height) compared to a rectangular window (sharper main-lobe width but higher side-lobe). You need to be careful with the requirements regarding the choice of window given the trade off between side-lobe height and main-lobe width.


why eigenvalues concerned in Harris corner detection? In Harris corner detection, as I see, the goal is to find out those points $(x,y)$ which makes $S(x,y)$ a large value for any direction of shifting, $$ S(x,y) = [\Delta x , \Delta y]M\left[ \begin{array}{c}
\Delta x \\
\Delta y \\
\end{array} \right]$$, while $$M = \sum_{(x,y)}\left[ \begin{array}{cc}
I_x(x,y)^2,& I_x(x,y)I_y(x,y) \\
I_x(x,y)I_y(x,y),& I_y(x,y)^2 \\
\end{array} \right]$$.
Here is my problem, as the wiki tells me that $M$ should have 2 large eigenvalues in order to make $S(x,y)$ large in any direction. WHY large eigenvalue? HOW?

Basic linear algebra. M is a symmetric matrix, so there's a decomposition:

$M=Q^{\mathsf{T}} \Lambda  Q$

where Q is a rotation and $\Lambda$ is a diagonal matrix.

We can rotate $[\text{$\Delta $x},\text{$\Delta $y}]$ by $Q^{\mathsf{T}}$

$Q^{\mathsf{T}} [\text{$\Delta $x},\text{$\Delta $y}]^{\mathsf{T}}=\left[\text{$\Delta $x}',\text{$\Delta $y}'\right]^{\mathsf{T}}$

giving

$S(x,y)=\left[\text{$\Delta $x}',\text{$\Delta $y}'\right]\Lambda \left[\text{$\Delta $x}',\text{$\Delta $y}'\right]^{\mathsf{T}}=\Lambda _{1,1} \left(\text{$\Delta $x}'\right)^2+\Lambda _{2,2} \left(\text{$\Delta $y}'\right)^2$

So $S(x,y)$ is large in both directions $\text{$\Delta $x}'$ and $\text{$\Delta $y}'$, if both eigenvalues are large.


Pattern matching methods: Based on Area and Based on features Good morning, I'm new to image processing and I'm currently working on a project for pattern matching and OCR more specifically.

In the state of art I've found the matching process can be classified into 2: Area based matching and feature based matching. The ABM approaches involve techniques like Cross Correlation which is basically and sliding window, meanwhile the FBM is based on features extraction in order to do the matching of the pattern based on a feature vector...

And here comes my doubt: The extraction of the characteristics can be generally done by algorithms like SIFT and SURF and RANSAC for the selection of the best matching points, but I came across the HOG descriptor which basically would be of good use for pattern matching without having to use any classifier at all, and it would be very simple to implement for the recognition of the characters by just comparing feature vectors.

So my question is, could this HOG descriptor be accurately classified as a method for FBM? OR where in the Computer Vision/Image Processing Area could this be classified?
And what would be the best and correct approach to using HOG in the OCR?

Thank you so much in advance!

Yes, if applied properly, HOG is a good feature extractor even for OCR. I could point you towards an OpenCV sample and a MATLAB digit classification tutorial which do exactly that. However, if you like to get state-of-the-art performance, I would suggest deep learning methods, more specifically convolutional neural networks. Many works demonstrated the power of these approaches on MNIST dataset. You might like to check this one or this one for some beginner's tutorials.


What is wrong with this GNURadio flow graph? I am simply trying to modulate a wav file sampled at 22.05kHz, then demodulate it.  The audio output is choppy and I get errors in the console output ("volk warning, no arch found... though I've run the profiler.  And aUaUaU....).  I guess the issue is with matching sample rates but I've gone through the graph several times and haven't found it.
Also maybe this is related: after I add the bottom row of components (channel model and beyond) the constellation diagram for the transmit side changes to flat.  Meaning, instead of a circle, it's a line (second screenshot below).  What's going on?

I wasn't alone it seems.  It is apparently a known issue and the solution here solved this for me as well: https://github.com/gnuradio/gnuradio/issues/1516


Time-spread multiple echoing of audio signals I have been reading papers related to adding echos into digital audio signals recently. I guess I understood the logic behind, yet I guess I am still having lack of DSP knowledge to figure out everything.

In order to add a single echo delay $\Delta$ into audio signal $S$ using convolution operator, echo kernel $k(n)$ is defined as:

$$ k(n)= \delta [n] + \alpha \cdot \delta [n-\Delta]$$

where $\delta [n]$ is Kronecker Delta Function, and $ 0 < \alpha \ll 1.$

Echoed sound signal $S'$ calculated using convolution operator "$*$" as:

$$S'=S \ * k(n)$$

I have no problems until this part. My question is based on this publication, [Ko, Nishimura, Suzuki - 2005]. It is based on echo data hiding for audio watermarking (steganography), but a very good source to understand sound echos. Authors say that multiple echos provide better quality for HAS than single echos (see Fig. 1.)

A pseudorandom sequence is generated to spread multiple echos in time. I generate it in MATLAB with the following code:

rand('seed', key);         % set seed of random sequence
a = rand(length, 1);
prseq = (a > 0.5)*2-1;     % convert into -1 and +1

In next step, time-spread echo kernel is defined using PN Sequence as: (see Fig. 2.)

$$ k(n)= \delta [n] + \alpha \cdot P[n-\Delta],  0 < \alpha \ll 1 $$

Schemetic explanation is made in Fig. 3. below:

I tried this tecnique in MATLAB and I actually have two questions:

Audio signal is being corrupted when the number of echos increases no matter how small $\alpha$ value is chosen. In the same article length of PN Sequence is chosen as $1023$, and delay = 1 ms which is about 44 samples in 44.1 kHz audio signal. So how is it possible to provide "a real room" echos as it is said in the article? Do I miss anything?

I am adding my summarized MATLAB codes:

delta = 44;
alpha = 0.02;

% Single echo kernel
for n=1:delta+1
    k_single(n) = kdelta(n-1) + alpha*kdelta(n-delta-1);
end

% Time-spread echo kernel without using kronecker delta:
kernel = [1 zeros(1,delta-1) alpha*prseq'];

echoed = conv(audio_data, kernel);

function [ out ] = kdelta( x )
    out = double((x==0));
end

Alternatively:

single_echo = signal + alpha*[zeros(1,delta) signal(1:length(signal)-delta)];

I was adding echos in an audio signal by shifting and adding as I mentioned above "alternatively". I figured out that it would take so much time to use it for multiple echos, so I thought I had to use convolution. But conv(x,y) function in matlab changing data size up to length(x)+length(y)-1 according to definition of convolution. Is it a good idea to just to crop it like echoed(1:length(audio_data)) to keep the length fixed? What is best option to accomplish multiple echoing as fast as possible and using less memory?

Excuse me for such a long question. I wanted to explain my problem in details, so it made my question longer.

Thanks in advance.

I'm not that familiar with the technique but reading up it looks like you basically just want an echoic convolution kernel of some sort which you can apply. I don't think the kernel you are building there is correct - if you want, for example, a single echo then the convolution kernel would need to contain two impulses, one at zero offset and one offset into the kernel by the time delay of the echo. For further echoes you would want further impulses, again offset into the kernel by their delay amount. You can have more sophisticated echo kernels but I would suggest this as a starting point

In response to question 2 I would not crop your signal unless it is particularly important to keep the length the same. Keeping the tail will make the deconvolution more accurate and unless your kernal is huge it won't be a very noticeable extension of the signal.


A filter with low phase shift or no phase shift at all? Which one? I'm seeking a good advice for a filter. Can you help me?

Right now I made my own simple filter. Not sure what I'm doing, but it seems to work.

function testFilter (t, y, p)

    % Forward filtering
    % Do first filtering by jumping k/2 elements in the future
  l = length(t);
  ye = y;
    for i = 1:l
      sum = 0;
      k = 0;
      for j = 0:p-1
        if(i + j <= l)
          sum = sum + ye(i + j);
          k = k + 1;
        end
      end
      if(k != 0)
        ye(i+floor(k/2)) = sum/k;
      end
    end

    % Do the first p elements!
    for i = p:-1:1
      sum = 0;
      for j = 1:i
        sum = sum + ye(j);
      end
      ye(i) = sum/i;
    end

    % Done!
  figure
    plot(t, ye, t, y);
    legend("Filtered", "Noisy")
end

It's some kind of moving average, ich.
The reason why I selected this is beacuse a low pass filter caused phase shift and I don't like that.

Here is a plot that shows a example how my filter algorithm works.
I think it's OK, but can it be done better without phase shift?
Do you think Fast Fourier Transform is a much better filter in this noisy case?

If there is no filter that cause phase shifting, can you recommend a good filter for me then with low phase shift?

If your filter is $h[n] \leftrightarrow H(e^{j\omega})$, then the DTFT of the time reversed filter $h[-n]$ is $H^*(e^{j\omega})$. Consider the result of passing your signal through the combined filter $h[n]*h[-n]$, which has a DTFT of $H(e^{j\omega})H^*(e^{j\omega})=|H(e^{j\omega})|^2$ and has zero-phase. In this setup, I've assumed that $h[n]$ is causal so that makes $h[-n]$ anticausal, and this makes the overall filter $h[n]*h[-n]$ noncausal and not applicable for real time applications but can be used for offline processing.

This is the idea used in the MATLAB command filtfilt. You pass filtfilt the filter coefficients and it runs the input signal through the filter and then the time reversed filter. Of course, you can also just as well do that by hand.

Here is a dumb example of a noisy sine wave and a moving average filter just to illustrate.


Signals in real world?especially natural signals? Do all the signals in real world ,that are of interest in engineering(for measurements or processing)are only continous/analog or there any other alternatives also?
For example birds make voice,animals make voices are examples of continous signals

Few natural discrete signals in my opinion
Number of teeth and number of bones in human
Number of mountains in a range,height of a mountain/hill
Please kindly give me other examples from natural world

Please kindly guide me if i am wrong in my understanding of continuous and discrete signals

Eventhough most engineering signals are continuous, not all are such. For example, statistics of many kinds provide discrete data such as daily stock market indicators, computer network stats, hourly number of customers in a queue, etc...

Of course, audio signals are transmitted as acoustic pressure waves through air (typically) and pressure is a continuous modeled physical variable.


Difference between Leaking Factor and Forgetting Factor I am using a Recursive least square adaptive filter to process electromyography signals and it is working decently so far. I decided to implement an LMS adaptive filter as a noise cancellation, so that I can compare the results, however, going through the matlab documentation for the LMS filter and seeing the LeakingFactor $$0<\alpha<1$$ and the ForgettingFactor of the RLS filter$$0<\lambda<1$$ I am now confused if there is an actual difference between those two parameters or is the difference only in the names?

The play similar role in those algorithms - the ability to forget the past and adapt to current reality.

In the LMS, the classic implementation has $ \alpha = 1 $.
Namely the optimal weights at any point are function of all inputs.

The Leakage factor allows to weigh the past differently in a damped manner which over times means the far past has practically no significance on the current result.
The other side means the ability to integrate data against noise is decreased which is the same balance the user stands when selecting the parameters of the RLS.

Namely, as usual, the balance between high bandwidth filter which adapts fast yet is sensitive to noise or the low bandwidth filter which is slow yet if the data is stationary can handle high energy noise.


Mean delay for a M tap harmonically weighted filter I am trying to derive the mean delay for $M$ filter taps harmonically weighted:

$$y[n] = \frac{M\,x[n] + (M-1)x[n−1]... + 1 \cdot x[n−M+1]}{\tfrac12 M (M + 1)}$$.

For a uniformly weighted filter, I used group delay formula and solved it but I am not able to solve this question in a similar way. If I find frequency response $H(e^{j\omega})$ for the above mentioned equation, I end with a arithmetic-geometric progression

$$H(e^{j\omega})=M+(M-1)e^{-j\omega}+...e^{-j\omega(M-1)}$$.

I don't know whether my approach is right or not.

Also, if this is a non-linear phase filter, how should I solve it?

Thanks in advance.

The coefficients of this FIR filter do not posses any type of symmetry, hence it cannot be a linear phase filter. However, you can still compute its phase response either analytically (or numerically) as:

$$\phi(\omega) = \tan^{-1}( \frac{ \mathcal{Im} \{ H(e^{j\omega}) \} }{\mathcal{Re} \{ H(e^{j\omega}) \} } ) $$

Then you can compute its derivative (again either analytically or numerically) to obtain the group delay function

$$ \tau(\omega) = -\frac{d\phi(\omega)}{d\omega} $$

from which you can compute the mean delay according to its definition.


Example of Entropy and Channel Capacity Computation Can you help me on verifying if this computation of entropy is correct and on understanding its meaning?

I am not sure of the result especially because it is equal to 0: it means that we cannot transmit any bits/s on that channel. Why?

Since all transitions are equally likely, there is no information in the channel output about the source. Hence, the mutual information, and, consequently, the channel capacity, are zero in this case. You might as well look at the output of a random number generator and from that guess the source symbols.

Expressed in formulas you have

$$H(Y)=H(Y|X)\tag{1}$$

because $X$ and $Y$ are independent, and, consequently,

$$I(X,Y)=H(Y)-H(Y|X)=0\tag{2}$$

i.e., there is no information in $Y$ about $X$ (which makes sense because the two random variables are independent).

Since mutual information is symmetrical, i.e., $I(X,Y)=I(Y,X)$ we also have $H(X)-H(X|Y)=0$, i.e.,

$$H(X)=H(X|Y)\tag{3}$$

which means that the uncertainty about the outcome of $X$ is the same whether or not the channel output $Y$ is given.


How Is the Bandwidth of an Image Reduced by an Octave in Gaussian Pyramid Procedure? Reference paper 1. For arguments sake, let's say that the maximum and minimum frequency in an image are 256Hz and 1024Hz. When I perform filtering using the gaussian function, the "bandlimit" reduces by an octave (1/2). Does this mean the new freq. range will be from 128Hz to 512Hz? Or does it mean that it will be 256Hz to 512Hz, but the bandpass signal (While calculating the Laplacian pyramid) will now have a freq. range 512Hz to 1024Hz (which is 1 octave)?

Also, apart from the low pass filtering, since downsampling is also performed, shouldnt the image bandwidth reduce by more than an octave? I'm guessing that the LPF operation is what causes the reduction in bandwidth

The model is simple.
Assume we have data I0 on the frequencies range of $\left[ 0, 1000 \right] $ [Hz].

We have 2 Black Boxes:

BlackBox 1 - LPF which always has a cut off at the middle of the spectrum.
BlackBox 2 - Decimation by 2.

What happens if we apply BlackBox 1 on the data?
We'll get an LPF version of it we'll call it L0 with data on the range $ \left[0, 500 \right] $ [Hz]..
We can also have H0 = I0 - L0 which is vasically all data on the range $ \left[500, 1000 \right] $ [Hz].

Now, be decimating I1 = decimate(L0) we have new data on the range $ \left[ 0, 500 \right] $ [Hz].
Applying the same tricks as above (Using the same black boxes which is what's nice about it) can yields the following series:

L# - Always contains the lower half frequency of the # Level.
H# - Always contains the higher half frequency of the # Level.

If you see the amount of data is reduced by factor of 2 each level (Octave).
Another nice property that you can reproduce the data using the lowest level LPF and the HPF of all other levels.


How to detect decaying oscillations in a signal I have a sampled signal which contains (among other things) some damped oscillations. From the underlying theory I have reasons to believe that it can be characterized by some specific frequencies of oscillation $\omega$ and rates of decay $-\sigma$. I believe the signal contains components of the form $f(t)=e^{\sigma t+i \omega t}$. It might be that these are "there" only in some part of the signal, because at late times other components win over (due to exponential damping of these oscillations). I am interested in finding these dominating timescales. I need a hint how to proceed. There are plenty of tools available, like Fourier analysis, wavelets etc. However, I don't know which is the most appropriate for this task, as I have no experience with any of them, except maybe the Fourier transform. Below I include a graph of the signal that I want to analyze.

http://i.imgur.com/fZwhdkq.png

Your signal model lends itself particularly well for so called parametric frequency analysis. The basic idea is that you find the autocovariance matrix of your signal and identify its eigenstructure. The eigenvalues correspond to the exponents of a sum-of-exponentials model. In your case you will find conjugate pairs of eigenvalues with non-zero imaginary part describing the exponential decay.

This thesis looks like a good introduction to the topic:
http://zet10.ipee.pwr.wroc.pl/record/11/files/Leon_mon.pdf


Image Gradient: Light to Dark, Dark to Light From link:

  Black-to-White transition is taken as Positive slope (it has a
  positive value) while White-to-Black transition is taken as a Negative
  slope (It has negative value). So when you convert data to np.uint8,
  all negative slopes are made zero. In simple words, you miss that
  edge.

So all i want is to write a function so i can keep negative gradient values.
I know how to use Sobel.

Here is the sobel link where I got to know light-to-dark, dark-to-light are gradient directions, which are also referred as $+\text ve$ and $-\text ve$ directions.

IMPROVED Question:
Please Suggest me a logic of the function to retain only $-\text ve$ values and clip $+\text ve$ values to zero?

My code:

Sobel(grayImg, grad, CV_16S, 0, 1, -1, BORDER_DEFAULT);
         gradCpy = grad.clone();
         for(uint64_t i=0; i <(gradCpy.rows); i++)
    {
        for(uint64_t j=0; j<(gradCpy.cols); j++)
        {
            if ((gradCpy.data[((gradCpy.step/gradCpy.elemSize()))*i+j]) > 0) {(gradCpy.data[(gradCpy.step/gradCpy.elemSize())*i+j]) = (int32_t)-0.0;} ///nega gradient img
            //if ((gradCpy.data[((gradCpy.step/gradCpy.elemSize()))*i+j]) < 0.0) {gradCpy.data[((gradCpy.step/gradCpy.elemSize()))*i+j] = (int16_t)+0.0;}   ///pos gradient img
        }
    }
    convertScaleAbs(gradCpy, opImg2);
    imshow("negative grad", opImg2);

But i am not getting expected results. As per Sobel link, I want to get edge at black-to-white that is $+\text ve$ gradient and another one at white-to-black that is $-\text ve$ gradient

The UINT8 type is limited to integer numbers on the range [0, 1, 2, ..., 255].
Hence negative values are clipped into 0.

A solution could be wither use other types (Floating Points) or scale and shift.

For instance, given a Floating Point image with values imageMax and imageMean you can scale it into the [0, 1, 2, ..., 255] by doing:

Create only positive numbers: image - imageMin.
Scale it into the range: round(255 * (image - imageMin) / (imageMax - imageMin)).

Now the conversion won't loose data.

Others have different approach, set 128 to be zero.
Namely, add 128 to the image, clip it into [0, 255] and then convert to UINT8.

Enjoy.


What's the unit of spectral intensity in calculation of RGB tristimulus? According to CIE color space and relation of RGB with XYZ, we can express

$$
R = \int_0^\infty I(\lambda)r(\lambda)d\lambda
$$$$
G = \int_0^\infty I(\lambda)g(\lambda)d\lambda
$$$$
B = \int_0^\infty I(\lambda)b(\lambda)d\lambda
$$

where $r(\lambda),g(\lambda),$ and $b(\lambda)$ are the corresponding color matching functions and $I(\lambda)$ is spectral intensity function.

Now, I am trying to calculate RGB values of hyperspectral image captured by hyperspectral camera. And from the hyperspectral image, I can get spectral intensity values between 0 - 65535. But I am not sure about unit or range of the intensity function $I(\lambda)$ above. Is it values between 0 - 1?

You are probably working not with CIE RGB, but CIE XYZ, so these are properly denoted the XYZ tristimulus values, not RGB. Their absolute values do not matter if you want to determine the chromatacities ("color", as we commonly understand it), because you take the ratio of the components to their sum.

Once you have the XYZ tristimulus values or xy chromaticities, you can convert to something like sRGB to yield usable RGB values. You might need to scale the result to 0-255.


In OFDM, does $N N^H$ equal $I$? I'm still new in OFDM and reading about it, but I have a question.

In OFDM system, when having $N$ data symbol assigned to sub-carriers, $N$ should be orthogonal with each others, which means $NN=0$. That's ok.  What's about that matrices $N$ itself?

Is there any property in OFDM saying that $N N^H$ should equal $I$ (where $I$ is the identity matrix)?

You're conflating two different definitions of orthogonal.

Orthogonal in the OFDM sense means that the subcarriers are mutually orthogonal to one another from the perspective of a correlator-based receiver (the optimum receiver for the AWGN channel).
An orthogonal matrix $\mathbf{A}$ is one whose columns are all unit vectors that are jointly orthogonal to one another, such that $\mathbf{AA^H} = I$.

These two concepts aren't related; there is no constraint on the modulation applied to the various subcarriers in an OFDM symbol. The orthogonality is imposed via careful choice of the symbol period and subcarrier spacing, typically implemented using DFTs since it conveniently provides the desired time/frequency grid.


Periodogram Averaging FOR CONSISTENT ESTIMATION OF PSD

Let, $x_i(n)$ be a sequence of random variables for $i=1,2, ... , K$ uncorrelated realizations and $\widehat{P}_{(i)}^{per}(e^{jω})$ is their corresponding periodograms. Then, the average of these periodograms is -

$$\widehat{P}_x(e^{jω})=\frac{1}{K}\sum_{i=1}^K\widehat{P}_{(i)}^{per}(e^{jω})$$

Let, $x_i(n)$ be a WSS process also.

Question:

Are the following relations correct?

$E\left\{\widehat{P}_x(e^{jω})\right\}=E\left\{\widehat{P}_{(i)}^{per}(e^{jω})\right\}$
 and $Var\left\{\widehat{P}_x(e^{jω})\right\}=\frac{1}{K}Var\left\{\widehat{P}_{(i)}^{per}(e^{jω})\right\}$

I found these two equations in my textbook (Statistical Digital Signal Processing and Modeling, Monson H. Hayes, Page 412).

I am failing to understand these relations (meaning I don't get how you derive them). Can someone please explain them for me?

Since the data records $x_i$ are uncorrelated realizations of the same random process, $\hat P^{per}_{(i)}(e^{j\omega})$ are all uncorrelated random variables with identical means and variances (given by the mean and variance of the Bartlett Method's PSD estimator). So, for any fixed $\omega$, $\mathbf E[\hat P^{per}_{(i)}(e^{j\omega})]$ are equal $\forall i$ and also $\mathrm{Var}[\hat P^{per}_{(i)}(e^{j\omega})]$ are equal $\forall i$. So,
\begin{eqnarray}
\mathbf E[\hat P_{x}(e^{j\omega})] &=& \mathbf E\left[\frac{1}{K} \sum_{i=1}^K\hat P^{per}_{(i)}(e^{j\omega})\right] \\
&=& \frac{1}{K}\sum_{i=1}^K \mathbf E\left[\hat P^{per}_{(i)}(e^{j\omega})\right] \\
&=& \frac{1}{K} K \;\mathbf E\left[\hat P^{per}_{(i)}(e^{j\omega})\right]\\
&=& \mathbf E\left[\hat P^{per}_{(i)}(e^{j\omega})\right]
\end{eqnarray}

and similarly,

\begin{eqnarray}
\mathrm{Var}[\hat P_{x}(e^{j\omega})] &=& \mathrm{Var}\left[\frac{1}{K} \sum_{i=1}^K\hat P^{per}_{(i)}(e^{j\omega})\right] \\
&=& \frac{1}{K^2}\sum_{i=1}^K \mathrm{Var}\left[\hat P^{per}_{(i)}(e^{j\omega})\right] \\
&=& \frac{1}{K^2} K \;\mathrm{Var}\left[\hat P^{per}_{(i)}(e^{j\omega})\right]\\
&=& \frac{1}{K}\mathrm{Var}\left[\hat P^{per}_{(i)}(e^{j\omega})\right].
\end{eqnarray}


What is the advantage of weighted median filter over median filter? I could not find any proof or comparison of weighted median and median filters. Yes the more occurring pixel contributes much but what should be the criteria of choosing weight function. And any explanation of cases where weighted median is strongly recommended to be used.

Median Filter replaces pixel value c with p where p is the median of pixel values in neighborhood of c

In the case of weighted median there are N $\left[ I_1,I_2,..,I_N \right]$ neighbor pixels,for each pixel there is also weight. Weighted median of that neighborhood is $k^{th}$ pixel  where $k$ is the minimum integer with $\sum\limits_{i=1}^k w_{i}\geq\frac{1}{2}\sum\limits_{i=1}^N w_{i}$

Weighting is a common approach to control the importance per pixel. In other words, one might require or desire different importance for different pixels. Then, the 50% percentile is estimated via taking the weights into account. Instead of treating each pixel equally, the sorting function is tweaked to consider the weights.

Imagine, we want to process pixel $p$ in image $I$, in the neighborhood $R(p)$. The radius of this local window $R$ is $r$. Number of pixels in such locality is $(2r+1)^2$. For each neighbor $q \in R(p)$, a weight $w_{pq}$ is associated. The typical choice of weighting is the affinity  of $p$ and $q$ in the feature map $\mathbf{f}$, which can be any feature, but in general selected to be intensity, color and etc. So we write

$w_{pq}=g(\mathbf{f}(p), \mathbf{f}(q))$

A reasonable choice for function $g$ is Gaussian (common preference for affinity measures):

$g=exp(-\frac{\lVert \mathbf{f}(p)-\mathbf{f}(q) \rVert}{\sigma^2})$

The weighted median operates as follows:

$p^*=\min{k}$ s.t. $\sum\limits_{q=1}^k w_{pq} \geq \frac{1}{2}\sum\limits_{q=1}^n w_{pq}$

This means that for all pixels before the median point $p^*$, the sum is roughly half of all weights summed together.


Given an image how should I detect this is a day picture or night picture? I have a 2D image which contains R G and B values. I think this should be purely based on averaging of pixels and selecting the threshold to decided whether this is day picture or night one.

Can you provide any kind of different ideas to solve this problem?

I am shifting this question from the https://stackoverflow.com/questions/28157035/given-an-image-how-should-i-detect-this-is-a-day-picture-or-night-picture?noredirect=1#comment44685520_28157035

This will work only if you're dealing with natural scenes (no artificial light sources).

I think there are two other possible ways:

1) looking for the percentage of dark pixels (define a threshold) that will tell you how many black areas are present in the image. Again, this will work only if we're not considering photos close to artificial light sources. Doing an average of the pixel value will lower the result if you have many lights in the scene, considering only dark areas could lead to better result (if there are enough in the scene).

2) Natural light have some difference in frequencies than artificial light, you could search which one is more present in the photo. I'm thinking at night pictures of cities for example, where lights are often in the yellow range. But i think this seems easier to say than implement, and looks more conditioned to the type of photo (natural panorama, city images, etc.).

P.S: i'm assuming the picture has been taken with the correct exposure time


Does such filter exist? Assume the following operation on input $x$ with filter $w$ of size $2L+1$:

$y[i] = \frac{\sum_{k=-L}^{L} x[i+k] w[k]}{\sum_{k=-L}^{L} x[i+k]}$

In other words, my desired filter does a convolution and divides by the sum.
Does this operation has a name (or implemented in conventional software such as MATLAB or python's Scipy)?

In case you are wondering why I need that, this is my problem: I'm having a 1-D temporal signal in which some textures (say edges) are to be discovered. The noise has attenuated locally (in time) some of the textures by some scaling factor which is unknown to me. Now I want to discover the potential time locations where a certain texture happens by normalizing the output of usual filter to the sum of absolute values of my signal (I've dropped the abs() since my signal is positive always).

It is not really a convolution, more a weighting. Your expression is $0$ homogeneous in $x$, as you divide $x$s values by $x$s.  Specifically:

it is not defined if $\sum_{-K,K} x[i+k] = 0$
if $x[k]$ is locally a constant, the output does not depend on it (because it will be equal to $\sum h[k]$)
if $h[k]$ is an averaging filter ($h[k] = \frac{1}{2K+1}$), the output will be $\frac{1}{2K+1}$

So globally I would not consider it as a  classical filter for $x$. But it could be seen as an adaptive weighting for $h$, based on a unit-normalized window $\chi$ with $\chi[n] = \frac{x[n]}{\sum_{|k|\le K} x[n]}$. As such, it is easily implemented in Matlab or Python. So far, I don't have a use for it. What was your purpose?


Multi-Time Window FFT One can achieve better resolution results by taking FFT of different sizes of the input signal. FFT size decreases as frequency increases, i.e. longer FFT length for lower frequencies and shorter FFT length for higher frequencies. I have tried to find papers on this topic but did not find any so far. Rational Acoustics has few brochures where it mentions MTW - Multi-Time Window FFT, but there is no mathematics behind it. Anyone can help me with the answer about underlying mathematics or some code implementation (C++, C, or Java)? In other words, how I can apply longer FFT, in my software, for lower frequencies and shorter FFT for higher frequencies to get uniform resolution as a result?

If your goal is to plot a magnitude spectrum with the frequency axis on a log scale, but with a roughly even visual resolution along that axis, then a single FFT might provide too low a density of plot points (without interpolation) at low frequencies, and more plot points than can be plotted on a line (without averaging) at high frequencies for a given print or pixel display resolution or smoothness.  If you use different FFT lengths for different octaves or sub-octaves in frequency, and select subsets of the results from each FFT, then you can maintain a lower delta in density of FFT result points when plotted on such a log scale.  How many FFTs you might want to use depends on the maximum variance you want in log frequency resolution of the final joined plot.

However, since the FFTs are of different lengths, then they are for different sets of data.  A sequence of longer FFT windows can be done with larger offsets if more time resolution isn't needed at low frequencies.  You will have the problem of how to join or blend all these different FFT results (usually subset segments of the results) so that the magnitude response might be uniform across FFT boundaries.  In total, you also end up with an overdetermined set of FFT results (but so will highly overlapped windows).

Even better results (in terms of even log frequency plot resolution) might be obtained by using some form of wavelet transform (a Morlet or Gabor wavelet or constant-Q transform, for example) instead of a bunch of semi-redundant FFTs.

But libraries for optimized FFTs might be more available on some platforms.  Each basis vector of a windowed FFT in a "multi-time" set can be considered a non-optimally-sized wavelet.  So, in some ways, using these so-called "multi-time window" FFTs is a "poor man's" wavelet transform.  But some of the many FFT results might be useful for other forms of more traditional DFT analysis of the signal input (in conjunction with the log plot), thus serving a dual purpose.

A partial psycho acoustic justification for increasing FFT size at lower frequencies (or using wavelets) is that, over a certain mid-frequency range, the human ear/brain takes a length of time to determine frequency or pitch of a sound roughly proportional to the period of the frequency.


Energy of a signal with decreasing amplitude I read that a infinite extension continuous time signal with decreasing amplitude will always be an energy signal. Though if you take the integral of, for example, f(t)=1/(t+1) we can see that is diverges to infinity. Doesn't the saying that if a signal has finite area then it is an energy signal? But that doesn't go the other way, meaning that if a signal is power signal then it has a finite area? I suppect something happens when you square the signal to calculate it's energy but I don't want to dive into caclulus and series convergence and such

I read that a infinite extension continuous time signal with decreasing amplitude will always be an energy signal.

That is wrong. It suffices to show a single counter example to show it's wrong: Let
$$s(t) = \frac1{\sqrt{\left\lvert t\right\rvert}},$$
which clearly is infinitely extending and always getting smaller the further you go from $t=0$. Then its energy is
\begin{align}
E_s &= \int\limits_{-\infty}^{\infty} \left(\frac1{\sqrt{\left\lvert t\right\rvert}}\right)^2\,\mathrm dt\\
&=2\int\limits_{0}^{\infty} \frac1{\left(\sqrt{t}\right)^2}\,\mathrm dt
\end{align}
and as "$\sqrt\infty$" doesn't ever exist, doesn't converge, and hence, this is clearly not an energy signal.
Your example works, too, and you get an "$\log\infty$" as a result. So, don't bother too much basing questions on a wrong claim.

Doesn't the saying that if a signal has finite area then it is an energy signal?

Well, "area" doesn't work great for signals that aren't real, or that aren't piece-wise continuous, or…
But yeah, if you can assign a finite area to it, then it's an energy signal, as that requires bounded absolute values. But thinking in "areas between curve and axis" really isn't the best level of abstraction for signal processing.


How to implement a $1$-bit quantizer in Simulink? I am trying to implement a delta-sigma modulator block diagram in Matlab Simulink but can't find a $1$-bit quantizer (which is required in the modulator) anywhere. Could anyone please tell me how to go about it?

You can see this example of Sigma-Delta A/D Conversion.


Closed form of $\mathcal Z$-transform : decomposition signal $x(n)$ The text of my exercise ask :

  Determine the closed form of the $\mathcal{Z}$-transform for this $x(n)$  $$ x(n) = \begin{cases} |n-N|  & \text{if 0<$n$<2N} \\ 0 & \text{elsewhere} \end{cases} $$

I see $x(n)$ like a composition of two signals

$$ x(n)= x'(n)+x''(n) $$
I use for greater semplicity $x_1(n)$
$$x_1(n)= (u(n)-(u-(N+1))$$
whence
$$x'(n)= (n-N) x_1(n) $$
$$x''(n)= x'(n-N)$$

Determination $X(Z)$ it easy, but if
the above equations are correct, for this reason i need post this question here.

$$X_1(Z)= \frac{z^{-(N+1)}-1}{Z^{-1} -1}$$
$$X'(Z)= z^{-1} \frac{{d}}{dz^{-1}} X_1(Z)- N X_1(Z)$$
$$X''(Z)= X_1 (\frac{1}{Z})z^{-N}$$

This is my solution.

$$X(Z)= X'(Z)+X''(Z)$$

Or you can proceed this way: \begin{align}\mathcal Z\left\{x(n)\right\}= X(z) &= \sum_{n=-\infty}^\infty x(n)z^{-n}\\
& = 0 + \sum_{n=0}^{2N} \left|n - N\right|z^{-n} + 0\\
& = \sum_{n=0}^{N} \left(-n + N\right)z^{-n} + \sum_{n=N}^{2N} \left(n - N\right)z^{-n}\\
& = N\sum_{n=0}^{N}z^{-n} - N\sum_{n=N}^{2N}z^{-n} - \sum_{n=0}^{N} nz^{-n} + \sum_{n=N}^{2N} nz^{-n}
\end{align}

And changing the boundaries for the second and fourth term, putting $k = n - N$, you get :

$$X(z) = N\sum_{n=0}^{N}z^{-n} - Nz^{-N}\sum_{k=0}^{N}z^{-k} - \sum_{n=0}^{N} nz^{-n} + z^{-N}\sum_{k=0}^{N} kz^{-k} $$

You can continue from there.


How do I identify a notch filter given the transfer function? I am studying for an exam, and a problem that my professor gave produces a transfer function that looks like the following:

$$H(z)=\frac{z^2-j}{z^2-\frac{1}{4}j} = \frac{\left(z-e^{j\frac{\pi}{4}}\right)\left(z+e^{j\frac{\pi}{4}}\right)}{\left(z-\frac{1}{2}e^{j\frac{\pi}{4}}\right)\left(z+\frac{1}{2}e^{j\frac{\pi}{4}}\right)}
$$

How am I able to determine:

That this is a notch filter,
At which frequencies the notches occur at?

His solution says that the notches occur at $\frac{\pi}{4}$ and $-\frac{3\pi }{4}$.

Remember that a notch filter must have transfer function zeros on the unit circle $z=e^{j\omega}$, i.e., there are certain frequencies that are completely attenuated. But note that zeros on the unit circle are not sufficient to characterize a filter as a notch filter (even though any frequency with perfect attenuation is called a notch frequency). E.g., most frequency-selective filters will have most, if not all, of their zeros on the unit circle (to get a good stop band attenuation). However, they are no notch filters. So for IIR filters you also have to consider the pole locations. And for notch filters, the pole angles are identical to or at least very close to the angles of the zeros. So for each zero on the unit circle you get a pole at the same or a similar angle, just with a smaller radius, because for causal and stable filters the poles need to be inside the unit circle.

That's why the filter in your example qualifies as a notch filter. All its zeros are on the unit circle and its poles are at the same angle as the zeros. The notch frequencies are simply the angles at which the zeros occur:

$$z_{01}=e^{j\omega_{01}}=e^{j\pi /4},\qquad
z_{02}=e^{j\omega_{02}}=-e^{j\pi/4}=e^{-j3\pi /4}$$


Time averaging denoising signal I would like to use time averaging technique for denoising vibration signal, using the below function, how do we choose the appropriate parameters D and N for optimal denoising  !
% sigav.m - signal averaging

%
% y = sigav(D, N, x)
%
% D = length of each period
% N = number of periods
% x = row vector of length at least ND (doesn't check it)
% y = length-D row vector containing the averaged period
% It averages the first N blocks in x
function y = sigav(D, N, x)
  y = 0;
  for i=0:N-1,
      y = y + x((i*D+1) : (i+1)*D);       % accumulate i-th period
  end
  y = y / N;

Your approach will only denoise anything, if your signal is actual periodic. That's typically the case if you use a sine wave or periodic random noise or sweep as an excitation signal for a measurement. It can also work, if you vibration source is periodic (e.g. an engine running an constant rpm).
In this case $D$ should be chosen to represent EXACTLY one period. $N$ is simply determined by the length of your data set. The longer the better.
This works best if the excitation clock and acquisition clock are phase locked, i.e. derived from the same clock source. If that's not the case, it may be required to do some sort of clock recovery or clock tracking first.


Understanding Nyquist rate I'm starting to learn signal processing, and am trying to understand Nyquist rate a bit better.

From what I understand, if I sample at a rate > Nyquist rate I'm supposed to have no data loss.

I'm testing on $\displaystyle f(t)=\sin(2t\pi)$.

To my understanding, its frequency is $1\ \rm{Hz}$, so the Nyquist rate is supposed to be $2\ \rm{Hz}$.

I tried FFT-ing the same function with sampling rates of $3\ \rm{Hz}$ and $10\ \rm{kHz}$ and got different results: (function + FFT displayed on same plot)

Why are the results different? Doesn't the fact that I'm over the Nyquist rate means that the FFTs are supposed to be the same?

This is the code I used in Octave:

x= [0:1/3:4]
x2 = [0:0.0001:4]
f = sin(2*x*pi)
g = sin(2*x2*pi)
subplot(2,1,1)
plot(x,f,x,fft(f))
subplot(2,1,2)
plot(x2,g,x2,fft(g))

1) The sampling rate has to be above twice the bandwidth, or twice the highest frequency.  The shorter the time you sample, the higher the sampling rate has to be ABOVE 2x highest frequency.

2) The signal has to be strictly bandlimited to at (or within) some frequency range below 0.5 sampling rate for no loss in information during reconstruction to occur.  Note that any sharp corners or steep edges in your signal imply an extremely high bandwidth.

Note also that any time-limited signal (e.g. not infinite in duration) has infinite bandwidth, thus will not be strictly bandlimited, so there might always be some error in reconstruction from samples.  This error can be reduced to the noise floor by sampling longer at a higher sample rate, which is why (in the real world) one should always use some reasonable margin above the theoretical Nyquist rate (which might be OK for imaginary, perfect, and infinite signals, as written on some chalkboard).


Confused about anti-aliasing in analog filters I'm a newbeginner when it comes to signal processing. I'm reading up about anti-aliasing, Nyquist and band-pass filters. There's something I don't get: Anti-aliasing in analog filters.

Lemme try to explain. I use a 8th order analog butterworth band-pass filter (centered at 17.5 Hz with a 25 Hz 3dB bandwidth) on the analog signal. So it cuts at 5 and 30 Hertz and ofcourse there's the transition-band before the stop-band. It is supposed to be an anti-aliasing filter. Where is this anti-aliasing happening?

I asked my teacher before he left why the butterworth was called anti-aliasing filter when it didn't really do it. He answered that since I was working on the analog signal, it was just a band-pass with no anti-aliasing because this was before ADC. The aliasing only occur after the ADC when the analog signal is made to digital and digital filters are used. WUT?

I mean, the anti-aliasing still happens when the signal goes through ADC, right? (well, I sample at 512Hz so I oversample the signal but stil...) So why is the filter called anti-aliasing filter?

Could someone explain this to me? I totally don't get it. Why is it called anti-aliasing filter and where does this anti-aliasing happen?

Aliasing does not happen in the analog world. It only happens when you sample a signal.

If you have an ADC with sampling frequency $f_s$, then you'll have aliasing if the input signal has frequency components larger than $f_N=f_s/2$. To avoid aliasing, it's common to low-pass filter the signal before it is input into the ADC. That way, you make sure that the digital signal out of the ADC has no aliasing.

Such an analog filter is called an "anti-alias filter", not because it will experience anything related to aliasing itself, but because its purpose is to avoid aliasing in the ADC.


What is the explanation for the inverse FFT of a cross-correlated signal being shifted by some amount in the time domain? I created two signals in the time domain: signal A which is a square wave, and signal B which is a right-angled triangle wave, as shown in the picture below.

I took the FFT of these signals, multiplied channel 1 FFT with the conjugate of channel 2 FFT. Then I inverse-FFT'd to observe it back in the time domain, and I observed this:

Now, instead if I did the correlation directly in the time domain using signal in Python like this:

txcorr = signal.correlate(wave1, wave2, mode = 'same', method='direct')

I see this output:

So, my assumption is taking the conjugate of a signal in the frequency domain reverses the time axis? Could someone give me a more mathematical background or where I could read to understand what is happening exactly in the frequency domain? Of course, when I simply do the multiplication in the frequency domain without taking the conjugate of a signal, it produces the shape expected

Additionally, this is probably the biggest confusion: Why is the "maximum point" of the cross-correlated signal in the time domain not at 0.3s? This is what I expected from what cross-correlation page in wikipedia showed:

The issue is that the FFT version always assumes that the indices start numbering from $n=0$.  The result of the FFT version will be periodic, with period the same as the FFT length.

To show the same result as the correlation version, you need to use numpy.fft.fftshift.

Here is my attempt to generate your original signals:

Then the two versions of correlation:

which seems to do what you're seeing.

Whnn I apply fftshift to the FFT version, I get:

Code Only Below

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
t = np.linspace(-0.2, 0.5, 500, endpoint=False)
t2 = np.linspace(-0.2, 0.5, 250, endpoint=False)
t3 = np.linspace(-0.4, 1.0, 999)
square = (signal.square(2 * np.pi * t - np.pi/2 , duty=0.2) + 1)/3
triangle = np.zeros(500)
triangle[250:500] = (signal.sawtooth(2 * np.pi * t2 - np.pi, 0 ) + 1)/2

txcorr = signal.correlate(square, triangle, mode = 'full', method='direct')
txcorr2 = np.fft.ifft(np.fft.fft(square, 999)*np.conj(np.fft.fft(triangle, 999)))

plt.figure(0)
plt.plot(t, square)
plt.plot(t, triangle)
plt.figure(1)
plt.plot(t3, txcorr)
plt.plot(t3, txcorr2)
plt.figure(2)
plt.plot(t3, np.fft.fftshift(txcorr2), '.')
plt.plot(t3, txcorr)


How to Remove a Glare / Clipped Brightness from an Image? I have an image

Is there any way of removing the bright white spots ? Please help thanks

Edit:

After operating with gaussian and then displaying using imagesc get the following output which clearly shows the bright red spots How do i get rid of them

Red Channel :

Green channel:

Blue channel:

Edit 2:

Defect detection using Gabor filter

Its Histogram :

How to calculate its appropriate threshold adaptivily.?

Let's assume glare portions are the only saturated areas in the image.
Detection can be performed by thresholding the intensity (code in Mathematica):

saturated = Binarize[ColorConvert[img, "Grayscale"], .9]

Then we need only to replace the portions of the image around the saturation mask (enlarging the mask is done by the morphological function Dilation). Inpainting using texture synthesis (using the function Inpaint) seems to work well in this example, though I can not test it as input to your defect-detection algorithm:

Inpaint[img, Dilation[saturated, DiskMatrix[20]]]


Shazam / SoundHound Music Analysis - Recommendation for Reading I need to explore some area around music frequency analysis.
A good example is the way Shazam / SoundHound app works (mapping frequency domain, comparing, creating some vector map, handling noises and errors).

I am looking for a while for a practical book or articles about this. By practical i mean with a simple language/code examples, etc.

No Sigmas, and integrals for infinity that takes you back to university.

It's really hard to find any music/audio analysis practical material.

The basic idea is called - Music Finger Printing.

Searching for it will yield many results.
I'm attaching few good ones I found:

How Does Shazam Work?
How Does Shazam Work? Music Recognition Algorithms, Fingerprinting and Processing.
Creating Your Own Shazam (Identify Songs) with Python Through Audio Fingerprinting (The YouTube Video).
Audio Fingerprinting with Python and Numpy.
Fingerprint Audio Files & Identify What's Playing.
Playing with Shazam Fingerprints.

Reading all those won't teach the exact rigorous Math, but it will give you intuition and tools to start with.


Causal system and Physical Systems According to the Paley-Wiener criterion, a system is causal if satisfies:

$$\int\limits_{-\infty }^{+\infty }{\frac{\ln (|H(f)|)}{1+{{f}^{2}}}}df<\infty$$

So I want to know

This equation is related to LTI system only?
Only casual system is realizable?for linear and non-linear systems?
Which kinds of systems can we create?(only Paley-Wiener criterion can be enough?)

The Paley-Wiener criterion defines a condition on the magnitude spectrum of a causal time-domain function. So if the Paley-Wiener criterion is satisfied for a given $A(\omega)=|H(\omega)|$, we know that there is a causal function with magnitude spectrum $A(\omega)$. It should be noted that the Paley-Wiener criterion is only applicable to square-integrable functions $A(\omega)$. If $A(\omega)$ is not square-integrable then the criterion is neither necessary nor sufficient.

A frequently occurring misunderstanding is that, given a frequency response $H(\omega)$, we can check the causality of the corresponding system using the Paley-Wiener criterion. This is generally not the case. If the criterion is not satisfied we know for sure that $H(\omega)$ is not the frequency response of a causal system. But if it is satisfied we only know that there must be a causal function with magnitude spectrum $A(\omega)=|H(\omega)|$, but we still don't know if the given $H(\omega)$ corresponds to a causal system. This depends on the phase of $H(\omega)$, and the only thing we know is that it is possible to find a phase response $\phi(\omega)$ such that $H(\omega)=A(\omega)e^{j\phi(\omega)}$ corresponds to a causal system.

To answer your questions:

Yes, it only applies to LTI systems because only LTI systems are fully characterized by a frequency response.
Yes, only causal systems can be realized because we cannot look into the future.
This question is a bit broad. There are realizability conditions for certain restricted sets of systems (e.g., passive LTI systems), but in general, causality is definitely a necessary condition. However, for the above mentioned reason we cannot just take a given $H(\omega)$, check the Paley-Wiener criterion and conclude that we can realize $H(\omega)$ as a causal system, because we didn't take into account the system's phase.


Output gain when upsampling and downsampling Do the processes of upsampling and downsampling affect the magnitude of the transform of a signal? And if not, why am I seeing everywhere that a filter with gain different from 1 is applied after up/downsampling a signal?

EDIT: For example in Discrete-Time Signal Processing by Oppenheim-Schafer (Section 4.6, 2nd edition):

Also, in the summaries of the subject that my teacher did himself, up- and downsampling are always followed by a lowpassfilter with gain different from 1.

All of the answers above are good in their own way. The reason that the Oppenheim-Schafer book and other similar resources use a gain of $L$ for upsampling and a gain of $1$ for downsampling is to maintain a constant amplitude for the signal before and after the operation.

  Do the processes of upsampling and downsampling affect the magnitude of the transform of a signal?

The short answer is that it depends on whether we are upsampling or downsampling. Upsampling does alter the magnitude of the transform (more on this in a minute). Downsampling doesn't alter the magnitude of the transform (assuming a perfect anti-aliasing filter).

For me, this is easiest to understand with a concrete example. So, let's consider a simple sinusoid, $x(t) = A \cos ( 2 \pi f_0 t + \phi )$. Let's sample this at some sampling rate $f_s$ higher than the Nyquist rate (i.e. $f_s > 2 f_0$). This results in a sample sequence

$$ x[n] = A \cos \left( \frac{2 \pi f_0}{f_s} n + \phi \right). $$

The process of upsampling consists of inserting zeros in between every other sample. A mathematical way to express this is

$$ x_{\uparrow 2}[n] = \left\{\begin{array}{cr} A \cos \left( \frac{2 \pi f_0}{2 f_s} n + \phi \right) & \mbox{$n$ is even} \\ 0 & \mbox{$n$ is odd} \end{array}\right. $$

and another equivalent way is

$$ x_{\uparrow 2}[n] = \frac{A}{2} \cos \left( \frac{2 \pi f_0}{2 f_s} n + \phi \right) + (-1)^n \frac{A}{2} \cos \left( \frac{2 \pi f_0}{2 f_s} n + \phi \right). $$.

If we note that $\cos(\pi n) = (-1)^n$, then we obtain

$$ x_{\uparrow 2}[n] = \frac{A}{2} \cos \left( \frac{2 \pi f_0}{2 f_s} n + \phi \right) + \frac{A}{2} \cos\left( \pi n \right) \cos \left( \frac{2 \pi f_0}{2 f_s} n + \phi \right). $$

From here, we apply the identity $\cos(\alpha) \cos(\beta) = \frac{1}{2} \left[ \cos(\alpha + \beta) + \cos(\alpha - \beta) \right]$ and use the fact that $\cos(\pi - \alpha) = \cos(\pi + \alpha)$. This yields the final form that we care about

$$ x_{\uparrow 2}[n] = \frac{A}{2} \cos \left( \frac{2 \pi f_0}{2 f_s} n + \phi \right) + \frac{A}{2} \cos \left( \frac{2 \pi \left(f_s - f_0\right)}{2 f_s} n + \phi \right). $$

This is the signal prior to passing through the interpolation filter. Notice that now we have $2$ sine waves. In general, if we upsample by $L$ then we have $L$ sine waves. Note that each one is weighted by $1/L$, this ensures that the summed amplitude of all sine waves has the same value as the original sine wave at the correct sampling instants and is zero elsewhere. After applying the interpolation filter, consider the transform of the sequence. It is the same obtained by sampling the original function at the higher rate, but it is scaled by $1/L$. Thus, to preserve the signal amplitude while filtering out the aliased images (in this case the additional sine wave component), we must scale the passband of our interpolation filter by a factor of $L$.

For the second part, assume that we are downsampling by $2$. We start with the same $x[n]$, but this time we choose $f_s > 4 f_0$ so that it isn't rejected in the anti-aliasing filter. The decimated signal is now given by

$$ x_{\downarrow 2}[n] = A \cos\left( \frac{4 \pi f_0}{f_s} n + \phi \right).$$

Unlike the upsampling case, this remains a single sinusoid within the Nyquist bandwidth. So, the magnitude spectrum is unscaled compared to the spectrum we would have arrived at by sampling the original signal at a slower rate.


How to make bode plot when output signal changes amplitude? When I do frequency analysis on my feedback controlled system and the controller is really tightly tuned, I get a frequency response that looks like this:

Blue is excitation signal and green is system output signal. Does this mean that my system has nonlinear behavior? Or this just a condition that would be visible when hitting a zero or a pole?

Smaller frequencies work fine and produce a sine wave output that is easy to interpret. It is higher frequencies that end up looking like this.

Since I intend to use this to derive a bode plot, I'm also curious how do I interpret this data and what is the actual magnitude and phase on this plot for this frequency? It seems to me like it is fairly random. How would I go about sweeping frequencies and making a bode plot when many frequencies in this region look like this? FFT maybe on both signals? But how exactly?

You are probably still seeing other frequency components from the transient of the system. I assume you that you simulated your system for a fixed number of periods of the excitation signal. If this is the case then it would also explain why you only see the transient at high frequencies. This would be because a fixed number of high frequency periods take less time than the same number of low frequency periods. Therefore there is not enough time for the transient to die out for the high frequency excitation.


How to filter out everything but a single frequency in the time domain? I'm new to signal processing and would like some insights about the best way to filter simple data without being too computationally intensive.

I have an audio signal and I want to extract a single frequency (21khz) from it.
This frequency represents whether something is on or off.
I figured a simple amplitude threshold would do the trick to guess the on or off state but I can't think of a good way to single it out so I can threshold it.

I thought about doing it in the frequency domain but using the FFT seems a bit overkill.

I'm sorry if this is obvious but I'm not familiar with signal processing.

Depending on the exact characteristics of the signal and implementation requirements (if any) I can think of a number of approaches to extract content at a "single" frequency. The most common:

Apply an FFT to calculate the spectrum of the signal and extract the frequency of interest.
Apply Goertzel's algorithm to effectively extract signal content at a single spectral line.
Apply a bandpass filter around the frequency of interest.
Apply an inverted notch filter around the frequency of interest.

Depending on the amplitude and frequency content of the signal, the additional noise present and other implementation requirements (e.g. should the signal be processed in real time or offline, should the implementation be digital or analog, etc.) a specific approach might be preferred.


Does direct-digital RF signal generation have a "baseband?" Semiconductor technology has advanced to the point that the digital representation of a modulated RF signal can be applied directly to a DAC for transmission. Is there, then, any distinction between a "baseband" signal and the RF transmitted signal?

Yes there is an important distinction which is independent of the semiconductor technology. In radio terminology, "baseband signals" are typically centered on "DC" and as such are required to be complex (such as I and Q) to represent a waveform that does not have a symmetrical spectrum about the carrier (basically nearly every modern waveform where modulation of amplitude and phase is involved) while IF signals and "RF Transmitted" signals can be real to represent these same waveforms.
In general "baseband" is the modulated waveform with DC (f=0) as the carrier while "IF signal" and "RF Transmitted signal" are representations with a non-DC carrier.


Signal magnitude area I'm confused by the SMA. I'm reading a paper which states that "SMA is the square root of the summation of the filtered acceleration signal". So I understand that it's $sma = \sqrt{\sum{samples}}$

I check the referenced paper, which says that "SMA is calculated after dividing area under the acceleration curve with the time interval". So $sma = {1 \over T} \sum samples$, which is some kind of average, which looks like what is on this StackOverflow question.

And the Wikipedia article is so new I'm not sure it's reliable.

In the end, what is exactly SMA as I couldn't find a proper definition online?

I've found a proper definition thanks to the following article:
http://www.rehab.research.va.gov/jour/2013/509/jrrd-2012-12-0233.html

This is what I've edited Wikipedia with:
$$f_{sma}= {1 \over {T}} {\int_{0}^{T} |x(t)-a_x|+|y(t)-a_y|+|z(t)-a_z| dt}$$

(It wasn't because of text-search but images-search. Algorithms are sometimes surprising)


Estimating signal strength from PCM-encoded signal I am processing rtl_fm output (in time domain), extracting Manchester encoded bits. I can detect preamble and, therefore, decode the message which follows the preamble.
This is how my signal looks like (here I have 2 messages, the rest is noise):

And this one of the messages zoomed in:

Now I want to know how strong was the signal that I've just decoded. I know what is the average amplitude for my symbols: "1" is 2900 and "0" is 1500 (for the amplitude range of -32768 .. 32767). Is it possible to know the signal strength in dB in this case?

In general, it is very hard to estimate the power of a signal after it has gone through an analog front-end and sampled. The front-end filters and mixes the signal, which can amplify or attenuate it. It's also common to have an automatic gain control (AGC) stage, which adjusts the signal's amplitude to match the ADC's input range and maximize dynamic range. In many cases, the digital back-end will also normalize the signal amplitude to a pre-defined range, such as -1 to 1.

If you understand every piece in the system very well, you may be able to estimate the signal's power from its samples. Another alternative is to use a scope or a spectrum analyzer to measure the signal as close as possible to the receiver's input.


How to determine initial values in Kalman Filter I have EEG signal and exract PSD feature from it, then must apply it a noise reduction algorithm, I used Kalman filter,

but the output signal in the paper is much smoother than my output and
vertical axis is different in my case

I think that is because of Kalman filter initialization; initial state and predict.

Here is my output and paper output.

Papers output:

Mine :

Here is MATLAB code:

N = length(z);          % number of Klamn filter iterations
Qfactor = 1;            % process noise mult factor
Rfactor = 1;          % measurement noise mult factor
F = [ 1   2            % update matrix
      0   1 ];
H = [ 1   0 ];            % measurement matrix
sigmaQ = 5e-5;
sigmaR = 1;
Q = sigmaQ^2 * [ 8/3  2     % process noise covariance matrix
                  2   2 ];
R = sigmaR^2 * [ 1 ];         % measurement noise covariance
P = zeros(2, 2, N);
x = zeros(2, N);
x(:,1) = [ 0
           0 ];
P(:,:,1) = Q;

for i=2:N
    [xpred, Ppred] = predict(x(:,i-1), P(:,:,i-1), F, Q);
    [nu, S] = innovation(xpred, Ppred, z(i), H, R);
    [x(:,i), P(:,:,i)] = innovation_update(xpred, Ppred, nu, S, H);
end

plot(x(1,2:N),'b');

The functions:

Prediction

function [xpred, Ppred] = predict(x, P, F, Q)
 xpred = F * x;
 Ppred = F * P * F' + Q;
Inovation :

function [nu, S] = innovation(xpred, Ppred, z, H, R)
 nu = z - H * xpred;                   %% innnvoation
 S = R + H * Ppred * H';                 %% innovation covariance
and Inovation update :

function [xnew, Pnew] = innovation_update(xpred, Ppred, nu, S, H)
 K = Ppred * H' * inv(S);                 %% Kalman gain
 xnew = xpred + K * nu;                  %% new state
 Pnew = Ppred - K * S * K';              %% new covariance

but the output signal in the paper is much smoother than my output and

I don't know which paper we're referring to, but it's far from uncommon that you'd smoothen your output for plotting – it would be adequate to specify how you did that, but notice that this sadly doesn't mean everyone adheres to that practice.

  vertical axis is different in my case

that might just be a scaling error, and is not necessarily bad. Again, we don't know which paper is about and whether you're working on exactly the same data, but the fact alone that your line is much thinner and looks like it has more points indicates that they applied some smoothing and/or normalization.

  I think that is because of Kalman filter initialization; initial state and predict.

I don't see indications for that. On the contrary, assuming things converge, I'd say your initial state is fine.


Find the frequency response of a system I'm trying to find the frequency response $$H(\omega) = Y(\omega)/X(\omega)$$
for this system- the signal equations are given:
$$y[n] = v[n - M] - g * v[n]$$
$$v[n] = x[n] + g * v[n - M]$$

I've tried using the z transform, and the fourier transform, but I just don't know how to get rid of the recursive terms to get it in terms of $Y(\omega)$ and $X(\omega)$.

The system is illustrated as follows:

I'll go through it in the $z$-domain. First, we find the transfer function $H_1(z) = \frac{V(z)}{X(z)}$. As you noted, in the time domain, $x[n]$ and $v[n]$ are related as follows:

$$
v[n] = x[n] + g * v[n - M]
$$

Take the $z$-transform of the above and you get:

$$
V(z) = X(z) + z^{-M}G(z) V(z)
$$

taking advantage of the convolution property, which states that $x_1[n] * x_2[n] \Leftrightarrow X_1(z)X_2(z)$, as well as the time-delay property, which states that $x[n-M] \Leftrightarrow z^{-M}X(z)$. Now, we can rearrange terms in the above to get:

$$
H_1(z) = \frac{V(z)}{X(z)} = \frac{1}{1 - z^{-M}G(z)}
$$

Now, find the transfer function $H_2(z) = \frac{Y(z)}{V(z)}$. They are related in the time domain as follows:

$$
y[n] = v[n-M] - g * v[n]
$$

Take the $z$-transform of the above to yield:

$$
Y(z) = z^{-M} V(z) - G(z)V(z)
$$

Rearrange terms again to get:

$$
H_2(z) = \frac{Y(z)}{V(z)} = z^{-M} - G(z)
$$

Now, you can get the overall transfer function $H(z)$ by multiplying the two:

$$
\begin{align}
H(z) &= H_1(z) H_2(z) = \frac{V(z)}{X(z)}\frac{Y(z)}{V(z)} = \frac{Y(z)}{X(z)} \\
&= \frac{1}{1-z^{-M}G(z)} \left(z^{-M} - G(z)\right) \\
&= \frac{-G(z) + z^{-M}}{1-z^{-M}G(z)}
\end{align}
$$

which has some nice symmetry to it. The overall response will be dependent upon the system $g$ in the feedback path. If you want to evaluate the frequency response, just let $z = e^{j\omega}$:

$$
H(\omega) = \frac{-G(\omega) + e^{-j\omega M}}{1-e^{-j\omega M}G(\omega)}
$$


Obtaining real numbers from FFT with same length as original signal In this article the real valued time domain signal is transformed to the frequency domain to extract some features like mean or variance.
But after transform to the frequency domain I calculate average hourly values (original measurements are 1/hour), so I need the signal in the frequency domain with the same length as the signal in the time domain.
Later I need features as real numbers (e.g. mean, variance), so I also need real numbers in the frequency domain. To this end, I think they use PSD. But after using Scipy's PSD (Welch estimation) I get much shorter array.
Question:

How can I get real valued array in the frequency domain from the signal with the same length as the original signal?

How should I use FFT / real FFT / PSD to get this (what arguments' values should I use)?

Would taking absolute values of the FFT output make sense?

Since you are reproducing the paper's results, read it carefully.
At the bottom of page 4 it says

This representation is formed by complex numbers, eliminating the imaginary part of each number in the frequency-domain signal. For this transformation, it is needed to calculate the power spectral density (PSD), as shown in Equation (3),
$$
P = \lim_{T\to\infty} \frac{1}{T} \int_0^T |x(k)|^2 dt \tag{3}
$$
where $P$ represents the energy from the signal, $T$ represents the length of the signal lapse and $x(k)$ represents the frequency-domain signal. The spectrum is normalized by the length of the signal.

The PSD is a real function with respect to frequency $k$, so it has the same length of FFT result $x(k)$ if you keepace with Eq. (3).
There are many methods to estimate PSD, but if you want to recreate the results in the paper, you should use exactly the same equation.