texts: Creating odd harmonics without the fundamental

texts:
I'm interested in creating a series of odd only harmonics, without the fundamental, for audio use.

As an example the input could be a sine wave at specific frequency $f$ and output is $3f$, $5f$, $7f$, etc. - but without the original $1f$.

Since the $1f$ (in a real audio system) is unknown and the input is multi-tone, just filtering the $1f$ out is not an option.

Any idea for replacing an input of $1f$ with an output of $3f$, $5f$, $7f$, etc. ?

Edit: Additional information from the deleted post

Example: Let's assume the fundamental to be anything possible between 100 and 4000 Hz. Whats more, what if the input is multitone (as most music is...). 50 and 425 Hz input simulatenously ??

Im looking for a mathematical function which can recreate any fundamentals odd harmonics, in a multitone input...


texts: Conclusions of sampling around Nyquist Rate

texts:
I'm trying to understand some results of playing around with sampling around a signal's Nyquist sampling rate. For my example, I'm sampling a $B=5\mathrm{Hz}$ wave over a 1 second period.

In the first image, I use a very high sampling rate of 1000Hz to emulate the exact signal (at least for graphing), $\sin(2 \pi 5 t)$. This sampling is subsequently used in the next graphs as a reference.
In the next, I use a rate just below Nyquist's, 9Hz. The resulting wave looks like a sine wave, but doesn't represent the highs and lows of the wave. I'm still tracking with why hitting the Nyquist rate is important from this image.
Next, I get confused with Nyquist: I sample at exactly $2B=10\mathrm{Hz}$, the lower bound. This produces nearly zeros in the data:

[ 0.00000000e+00  1.22464680e-16 -2.44929360e-16 -1.40896280e-15
-4.89858720e-16  6.12323400e-16  2.81792560e-15 -2.69546092e-15
-9.79717439e-16  1.10218212e-15]

(I'm not even sure if these are floating point errors or not). Yes, there is some oscillation, but it is very minor. It seems like we need to further increase the sample rate... Is the rate of $2B$ exclusive?
In the fourth image, I sample above at 11Hz, which above the Nyquist rate. However, this sampling looks very much like image #2 flipped about the x axis.  Again, there is some oscillation, but I don't see how we could reconstruct a 5 Hz wave from this, especially when it has nearly the same qualities as image #2.
Finally, I sample at $4B=10\mathrm{Hz}$ and at last, I can see the 5Hz wave I expect. Why did it take 4x to get here when the theory only states 2x?

I feel like the conclusions I've made above are incorrect (in probably a variety of ways). Would someone please help me explain where I'm wrong?

Jupyter Notebook available here: https://gist.github.com/t-mart/8a46e87938904c7f6a5a102ad6a2ef0e


texts: What does a complex matched filter mean? . A complex bandpass filter can be formed by multiplying a low pass filter by exp(jwt) getting an analytic filter (one sided in the frequency domain)  which is complex, the filter will have complex coefficients but it is implemented as two filters with real coefficients (coefficients multiplied by cosines for the real and the same coefficients multiplied bu sines for the imaginary) but how about the matched filter which is supposed to be a low pass filter centered around DC? how does that one look like?


texts:
For an FIR filter, with symmetrical tap values $h[N-1-n]=h[n]$,

why is the group delay $\frac{N-1}{2} T$ (where $N$ is the number of taps of the FIR filter and $T$ is the sampling time)?
Why is linear phase so important for the filter response?


texts: I have two audio signal that I want to compare using Matlab, my problem is that I can hear the difference between both loud and clear, but when use the function pwelch to compare their PSD ,I don't see much difference between them,any idea how can I compare them with a different methode ! thanks in advance !


texts: The result if order of two filter are reversed

texts:
Apply the Composite Laplacian Filter first, then apply the gaussian filter.
Apply the gaussian filter first, then apply the composite laplacian filter.

My work as below:
Here we assume the original image is function $f(x,y)$
$$ (1): (f-\nabla^2f)\ast G = f \ast G - \nabla^2f \ast G$$
$$ (2): (f \ast G) - \nabla^2(f \ast G) = f \ast G - f \ast \nabla^2G $$
I am wondering the statement $\nabla^2f \ast G = f \ast \nabla^2G$ is true or not.


texts: How to implement Walsh demodulation in matlab?

texts:
1     1     1     1    ==> 0
 1    -1     1    -1    ==> 1
 1     1    -1    -1
 1    -1    -1     1

I have this sequence of walsh spreading sequence, generated with matlab :
N = 4;  % Length of Walsh (Hadamard) functions
hadamardMatrix = hadamard(N)

SNR = 10;
h = [0.4 1 0.4];
s = randi([0 1],N,1);
sig = 2*s -1; % BPSK

r = filter(h, 1, sig);
receiveSig= awgn(x, SNR, 'measured');   %- received signal to equalize

I have one question:
How to make walsh demodulation after signal equalization ?
is it soft combining of LLR ?
thank you


texts: Real-time wavelet decomposition and reconstruction for ECG feature extraction

texts:
I need to locate R-peaks in an ECG signal. I'm using wavelets to extract QRS complexes:
First, I decompose the signal using a maximal overlap discrete wavelet transform with the Symlet 4 wavelet. This wavelet resembles the QRS complexes I want to detect. I then copy scales 4 and 5, because this band maximizes QRS energy. The wavelet coefficients on scales 1-3, and the level 5 approximation coefficients are set to zero. Then I perform an inverse wavelet transform on these modified coefficients. The result is a signal in the time domain that contains mostly QRS complexes, and has very high amplitude at the R-peaks.
I just use a peak detector on this signal to locate the R-peaks.

Everything is well-explained here: https://mathworks.com/help/wavelet/ug/r-wave-detection-in-the-ecg.html

This is the code I use:

wt = modwt(rawECGsignal_buffer,5);
wtrec = zeros(size(wt));
wtrec(4:5,:) = wt(4:5,:);
filtered = imodwt(wtrec,'sym4');

From what I understand, the Symlet 4 is a biorthogonal wavelet, which means that the decomposition and reconstruction can be implemented as a bank of high- and low-pass FIR filters.

Right now, I read the live ECG data into a 5 second buffer, and then perform the wavelet filter on the buffer (as described above). This seems rather inefficient, because I have to filter the entire buffer again when new data comes in.
Is it possible to implement this in a way that allows for real-time filtering of an ECG signal without buffering?
(I.e. single value in → wavelet 'filter' → single filtered/reconstructed value out)


texts: Discrete-time system: getting the true frequency from normalized frequency

texts:
Lеt's say that i have a plot where i can see the spectrum of some discrete signal and let's say that it's x axis is normalized (by normalized i mean $\omega/\pi$ on the $x$ axis). Now, i know that in order to get to the real frequency at some particular point on the plot, i need to know the sampling frequency $f_s$ and then multiply that particular point with $f_s/2$ in order to get the actual frequency.

Now, i know this is maybe something trivial, but i am having struggle to understand why is this so, i mean, why multiplication by $f_s/2$ gives us the actual frequency in hertz?

Any help appreciated!


texts: How to find the envelope of a wide-band signal?

texts:
I would like some feedback on possible techniques that one may use to determine the envelope of a broad-band time domain signal. I have heard anecdotally, that it is not as straight-forward as it seems, but I am not sure. (Some previous discussions about similar subjects exist here and here but I would like an explicit discussion on broad band signals here).

For narrow-band signals, one may compute the envelope via the absolute value of the analytical signal as such:

Let $s(t) = A cos(2\pi ft)$ be the real signal. Then, the analytical signal is:

$$
a(t) = Ae^{j2\pi ft } = Acos(2\pi ft) + jH(s(t)) = Acos(2\pi ft) + jAsin(2\pi ft)
$$

where $H(.)$ denotes the hilbert transform of the real signal. Then the envelope is given by $E = |a(t)|$.

What we are really doing is converting a real signal into phasor form, and taking that phasors absolute magnitude. In other words, since any real signal can be described as two phasors of equal magnitude rotating at both negative and positive angular frequencies, we simply discard of those phasors and work with the remaining one, and take its absolute magnitude.

Now to my understanding, this method works well for narrow-band signals, and not for broad-band signals. Thus my questions are as follows:

EDIT: I have edited and narrowed down my questions based on the feedback I have received:

Question: I am wondering, what ways one may use to compute the energy-envelope of the following broad-band signals, assuming the above analytical approach does not work well. (If it doesnt, why doesnt it?)

i) A chirp signal, at some carrier. (linear, exponential, etc)

ii) A very short duration windowed sinusoid, of length $T_p$. (So bandwidth is $2/T_p$)

iii) An OFDM signal at some carrier. OFDM is composed of multiple orthogonal (I do not care about it as much, but put it for completion).

iv) @Andrey's 'mixture of tones'. (I am not clear on this, but I am guessing it means simple summation of multiple sinusoids around a carrier?)

I cannot think of any other types of broad-band signals, but answering those cases might go far in grounding my understanding here.

Thanks.


texts: Estimate Filter Coefficients from the Result of Linear Convolution with a Known Signal

texts: If I have samples of input say x(1:500) and it passes through FIR filter with 9 taps and some unknown coefficients. The output y(1:508) is also known. The aim is to estimate the filter response in frequency domain. I know using Convolutional theorem $\hat{H} = (X^H*X)^{-1}(X^H*Y)$, but this requires fFT of full 500 samples, However I can only do fft of 128 samples at a time. How to estimate the filter response in frequency domain ?


texts: Which is the best method for comparing the histogram obtained after BOW?

texts:
Euclidian distance
Cosine similarity
Histogram Intersection
Or any other?

Consider the importance of the method in case of images. Eg, Cosine Similarity matches angle of two vectors...
What impact comparing just the angle or distance has on accuracy of BOW?
Detailed explanation is welcome!!!


texts:
Following a past question, I'd like to extrapolate a signal like the one below: (red is signal, blue is the the extrapolated )

t1 = 0:0.1:30;
t2 = 0:0.1:40;
s = @(t) sin(0.2*t)+sin(0.5*t+0.05*t.^2);
figure;
plot(t2,s(t2)); hold on
plot(t1,s(t1));

Linear predictive coding, or AR Burg's method are not successful here.
How can you extrapolate that?


texts: What is the difference between a lag compensation and PI control?

texts:
A lag filter or compensator has the Laplace transform transfer function
$$G_c(s) = \frac{s+z}{s+p},$$
with $-z < -p < 0$. In practice, the effect of lag compensation in feedback control, is to increase the DC gain of the open-loop system by the factor $G_c(0)=z/p$, which improves the ability of the closed-loop system to track low-frequency reference signals. The tradeoff is that the open-loop system incurs a phase shift of $-\pi/2$ at frequencies $z < \omega$.
However, in the limit as $p\rightarrow 0$, we are left with a PI controller, which has the form
$$G_c(s) = K_p + \frac{K_I}{s} = \frac{K_p s + K_I}{s},$$
which has infinite DC gain ($|G_c(0)|$), meaning the closed-loop system now achieves zero steady-state error in response to a step reference. The only tradeoff I can see compared to the lag filter is that there is a phase shift of $-\pi/2$ for $\omega < p$, but I don't see why this is a problem in practice.
Since the steady-state error is improved as the pole moves closer to the origin, why would you ever use a lag filter (with $z > 0$) rather than just place the pole at the origin?


texts: What is the time unit of the converted Channel Impulse Response?

texts:
I don't have a lot of background on wireless communication but as far as I understand, applying IFFT on Channel Frequency Response (CFR) would transform the data into time domain Channel Impulse Response (CIR), in which it shows the delay spread of the different multiple paths.

The CSI measuring tool I'm using measures 30 subcarriers resulting in 30 complex number CFRs per a single CSI trace. So when I apply 30-point IFFT on the CFR, the output is a CIR with 30 time-domain values.

My question is,

What is the time unit of the converted CIR?
How much seconds are they (the 30 impulse responses) apart?
Also, are the impulses like the echoes of the transmitted signal?
The number of multi-path is not determined. So I assume that an impulse response does not represent a single path, since the number of impulse response is fixed to 30. Then how should I understand the CIR value with respect to the multi-paths?


texts: down sampling and a factor of 2 without loss of information

texts:
I am trying to understand the limits of down sampling on a signal, i.e. how far can I go without loss of information or overlapping. It says in the solution manual of the Oppenhiem and Schafer book (on a question where down sampling is 2) that

"There is no loss of information if X(e^(jw/2) ) and X( e^(j(w/2) - pi) ) do not overlap."

Does this mean that my original signal is periodic and when I down sample means that my frequency response X(e^jw) becomes X(e^(jw/2) ) ?

Please help,


texts: Is sampling at double the desired reproduction frequency accurate?

texts:
This is probably a general principles question, though I'm thinking specifically in relation to sound, which is commonly sampled at a rate of 44.1 kHz in part because the maximum frequency the average person can hear tops out at around 20 kHz.

Given that the sampling theorem says the sampling frequency must be greater than twice the maximum frequency one wishes to reproduce, that means you need a sample rate of at least 40 kHz if you're trying to reproduce audio.  However, isn't it the case that the theorem is establishing only the lower-bound at which accurate reproduction becomes possible, rather than the sampling rate at which accurate reproduction is guaranteed?

By way of illustration, here are some examples showing a sine wave with samples taken at twice its frequency (in each example the source data and the interval between samples is the same, and the only thing that varies is the timing of the first sample in the interval):

In the 'In Phase' example, your samples occur exactly in time with when the wave is at a peak/trough, and you get an accurate capture that could reproduce the source wave (or something reasonably like it).  I assume this represents the ideal case, and rarely happens in practice.

In the 'Out of Phase' example, however, a terrible stroke of luck has the samples occurring at exactly the same time that the source data is crossing '0' amplitude.  This results in a capture of [0,0,0,0], which I assume can only reproduce as silence.  If so, that's not accurate at all.  I assume this represents the worst-possible case, and also rarely happens in practice.

In the 'Mid-Phase' example, the samples occur midway between the source data's peaks and troughs.  This captures a sine-wave pattern ([1,-1,1,-1]), but at half the amplitude of the source signal.  That doesn't seem like adequate accuracy for proper reproduction of the source.  I assume some variant of this scenario (where the sample rate is neither perfectly in-phase nor perfectly out-of-phase with the source) is the most common practical case, by far.

Thus it seems like if you sample at double the target reproduction frequency, most of the time your recorded data will have the wrong amplitude.  It may occasionally be accurate, but it's also just as likely to occasionally miss the signal completely.

Is this correct?  And if so, what is the sample rate one needs to use to make accurate reproduction guaranteed rather than merely possible?

If that can only be answered by imposing some constraints upon the definition of "accurate reproduction", let's say there must be zero chance of capturing silence wherever a signal is actually present, and that the recorded amplitude must always be accurate to within 5% of the actual source.


texts: Use samples of Fourier transform as DFT?

texts:
Consider the LTI system given by: $H(z) = 1 - \frac{1}{2}z^{-1}+\frac{3}{4}z^{-2}$

Let $x[n] = (\frac{1}{2})^nu[n]$ be the input to the system. We want to find the output for $n = 0,1,...,N_a$, using the DFT of a relevant portion of the input and samples of $H(z)$ on the unit circle, at $z = e^{j\frac{2\pi}{N_b}k}$, for $k = 0,...,Nb-1$.

Now, because $H(z)$ is FIR, the samples of $H(z)$ on the unit circle are actually the values of the DFT of length $N_b$ of $h[n]$. I can work out the appropriate length of the DFT of $x[n]$ to use and it would work and
using the circular convolution property, find the actual values of the output.

What if instead of a FIR, I'm given a causal IIR such as $H(z) = \frac{1}{1-\frac{1}{2}z^{-1}}$?

Can I use the samples of this new $H(z)$ on the unit circle just as in the first case? My guess is yes, because the system is stable (because it has a pole at $\frac{1}{2}$ and it is causal) and then the Fourier transform exists, so I can obtain the samples of the Fourier transform and use them as DFT of $H(z)$, (I believe) it would be like truncating $h[n]$ and then using its DFT.
Now, what if the system was unstable and causal, for example:

$H(z) = \frac{1}{1-2z^{-1}}$?

The Fourier transform clearly doesn't exist, and I can't truncate $h[n]$


texts: Why are males' pitch lower than females?

texts:
Hi :) I'm not sure this is the right place to ask this question, but I'll try anyway:
Here are two well-known facts, and I don't understand how they work together:

Males usually have deeper voices than females because they have longer vocal folds.

Short and relaxed vocal folds produce a deep voice, long and tense vocal folds produce a high-pitched voice.

So... given that males have longer vocal folds, isn't that supposed to mean that they produced high-pitched voices? What did I get wrong?
Thanks


texts:
Say you start out with $A=\left\{3,−1,0,0\right\}$, run fft() and end up with $B=\left\{2,3+j,4,3-j\right\}$.

What specific steps would need to be done to obtain $A=\left\{3,−1,0,0\right\}$ again using some pre-calculation on $B$, then fft(), then maybe post-calculation steps before obtaining $A$? I have heard that this is possible but have not been able to work out how.


texts: How do you do signal averaging on a realtime data?

texts:
Smoothing algorithms implemented by me

Exponential smoothing - https://www.openprocessing.org/sketch/771368
Moving Mean - https://www.openprocessing.org/sketch/540696

For example - a1+a2+a3+a4+... aN = A/N ..... is an average of all the a's which isn't same as ---> (((a/1+a2)/2 + a3)/3 +a4)/2 .... so on...

can anyone explain me how wrong I am or explain me signal averaging correctly?


texts:
I have files with extension .hpg and I want to perform image processing on those file as they are hp plotter file.
I tried googling but no luck there.
hence hope someone can help me tackle this with some clue.


texts: Distinguish between sinusoidal oscillations and rythmic pulses

texts:
I want to distinguish between a sinusoidal oscillations and rythmic pulses.
Which method can i use to do that?

I want to implement a detector that can distinguish between this two types of signals.

i have calculated the fast fourier transform of this signals. however to determine if it is a rhythmic pulses or sinusoidal oscillations, i have to fixe a threshold manually. It seems not good.


texts:
I am working on a university project in which I need to find a periodic net element in an image.

The net is a set of diagonal lines in a specific angle ( which I don't know in advance ) on a noisy image. They are hard to detect so I would like to use all the information I can get my hands on.

My current algorithm is based on using edge detection + hough transform to find lines at specific angle range.

I was wondering if there is any way to detect periodic signals in image? Something based on ff2 or something like that...

Thanks,


texts: Harmonics below oscillation frequency in the brain

texts:
Apologies if this question is trivial, I'm a neuroscientist and would hugely appreciate your help to correctly interpret my analysis.

I am recording EEG signals and am manipulating the brain using light flickering at 10Hz. I then compute the PSD of this signal using FFT of 4s long epochs. I am interested in what the manipulation does to frequencies around 5Hz.

My problem is that my 10Hz stimulation creates an artifact at 10Hz (light induces current flow). I see this clearly in the PSD and I also see higher harmonics at 20 and 30 Hz.  My question is whether I should expect this artifact to affect the 5 Hz frequency as well? I don't see a clear peak but I don't fully understand the theory so I'm not sure.


texts:
When we have FFT of Square wave of a pure real signal with high frequency(close to f(max), less than half of sampling frequency) we see only 2 symmetric high peaks but when we decrease the frequency we see the two peaks and peaks of harmonics and they are getting smaller and smaller, why? Also the highest peaks are getting away from each other, movig <- and -> , vice versa -> <- when we increase the frequency. why?

FFT of Sin wave has only 2 peaks no matter what and they also get away from each other when the frequency is low and get more and more closer when the frequency gets higher. Why? I would be glad if you can enlighten me!

Lets say our sample count is 64 and sampling frequency is 8kHz and signal frequency 2kHz. We should see 2 peaks one at 16 and one at 48.

The picture above is only the left half of my graph, there is also a symmetrical part on the right. My signal is pure real. These small harmonics dont appear at my above settings but when I decrease the freuency to 600 Hz for example I see harmonics like this for square wave and none for Sine Wave. So on high frequencies they disappear and low frequencies they appear, the lower the more of these harmonics. Also the peaks are  getting away from each other when I lower the frequency and getting closer when I increase the frequency. At 4 kHz I expect to see only 1 Peak, two of them will unite. Why? Also I couldnt test what happens if I run with 5kHz :( and when I run the code with lots of samples like 2048, the graph displayed makes no sense at all. Why?

I know I'm asking a lot questions but I would be glad if you have some time and help me!


texts:
I'm a musician and I'm making a script that takes a wave file and snaps each of its frequencies from the fourier transforms to the nearest musical harmonic. Thanks to help from another question I posted here, that part works, but what I need to do now is to experiment with some interesting ways to interpolate between the before and after version. I don't want to just do a crossfade because I want it to sound like one sound continuously morphing rather than one fading out and in. One idea I have is to snap one frequency at a time, maybe either in order based on the volume of that frequency by the end, from the lowest frequencies to the highest, highest to lowest, or random. Another idea was to interpolate them like image interpolation where if I have a frequency at 100hz with a volume of 255 and another at 120hz with a volume of 215, halfway between the interpolation I should have a frequency at 110hz with a volume of 235. If anyone knows of a function library that lets me do these then please let me know, otherwise I'll be making them and will need some help on getting started. How would I accomplish either of these?

Another thought I had was to reduce the number of harmonics over time so that in the beginning there are enough harmonics for the sound to be indistinguishable from the original but by the end there are only a couple. To get what I mean by the number of harmonics see the part of the script below where I calculate valid_ratios and valid_frequencies.

This is the sample I'm testing with (rename it to missile.wav):
https://my.mixtape.moe/iltlos.wav

Here is the script so far:

import numpy as np
import scipy.io.wavfile

samples, data = scipy.io.wavfile.read("missile.wav")

# fast fourier transform makes an array of the frequencies of sine waves that comprise the sound
data_fft = np.fft.rfft(data)

# the higher MAX_HARMONIC is, the more it sounds like the original,
# the lower it is, the more it sounds like an instrument
MAX_HARMONIC = 2

# generate list of ratios that can be used for tuning (not octave reduced)
valid_ratios = []
for i in range(1, MAX_HARMONIC + 1):
    for j in range(1, MAX_HARMONIC + 1):
        if i % 2 != 0 and j % 2 != 0:
            valid_ratios.append(i/float(j))
            valid_ratios.append(j/float(i))

# remove dupes
valid_ratios = list(set(valid_ratios))

# find all the frequencies with the valid ratios
valid_frequencies = []
multiple = 2
while(multiple < samples):
    multiple *= 2

    for ratio in valid_ratios:
        frequency = ratio * multiple

        if frequency < samples:
            valid_frequencies.append(frequency)

# remove dupes and sort and turn into a numpy array
valid_frequencies = np.sort(np.array(list(set(valid_frequencies))))

# bin the data_fft into the nearest valid frequency
valid_frequencies = valid_frequencies.astype(np.int64)
boundaries = np.concatenate([[0], np.round(np.sqrt(0.25 + valid_frequencies[:-1] * valid_frequencies[1:])).astype(np.int64)])
select = np.abs(data_fft) > 1
filtered_data_fft = np.zeros_like(data_fft)
filtered_data_fft[valid_frequencies] = np.add.reduceat(np.where(select, data_fft, 0), boundaries)

# do the inverse fourier transform to get a sound wave back
recovered_signal = np.fft.irfft(filtered_data_fft)

# amplify to ceiling
recovered_signal /= np.amax(recovered_signal)

scipy.io.wavfile.write("missile_output.wav", samples, recovered_signal)


texts:
In the implementation of block LMS, i need one clarification.

In the 3rd step as shown in the figure attached, the summation over a product of input, $\mathbf{u}$ and error, $e$, associated with each sample value in a block, gives rise to a single value (some integer).

Every first statement (formulae) in a line belongs to Block LMS algorithm and second statement (formulae) belongs to simple LMS Algorithm. What I want is only to consider first formulas in every line.

If all the initial weights $\hat{w}(k)$ are the same (all are zeros), then this constant value will be added to that at the end of each block which results in updates weights having the same value. i.e;

  [0 0 0 0 0] + x = [x x x x x].

If this process continues till the end of input, all the filter weights will be of the same value.

Is this correct?


texts: Energy and SpectralFlatnessMeasure taken for PCM audio VoiceActitivyDetection

texts:
I'm trying to implement the VAD algorithm in this paper:

http://ms12.voip.edu.tw/~paul/Papper/Steganography/iLBC/(VAD)Real-Time_VAD_Algorithm.pdf

I'm polling a microphone on an Android phone, and getting short[]'s which I want to perform operations on.  The audio data is 16bit PCM data.

I've gotten the basics of the algorithm done, but there are two issues I don't really understand yet.

1)  Step 3-1- of the above algorithm says to

  "compute frame energy".

I'm not entirely sure what this means.  Should I take the PCM data and convert it to Db's? Should I generate this off of the magnitude values derived from an FFT?

2) In section 2 (1) of the paper, There is an equation for Spectral Flatness given:

SPM = 10log10(Gm/Am)  or Spectral Flatness Measure = 10*Log10(GeometricMean/ArithmeticMean)

There is a brief description underneath which reads:

  "Where Am and Gm are the geometric means of speech spectrum
  respectively".

I understand what the general terms of this equation are, however I don't understand what I should be summing to generate the Gm and Am terms.  At first I thought it was the energy values computed in part 1) of my question, however rereading it seems to instead relate to the energy present only in the speech spectrum portion of a given frame.  If thats the case what range should I use and how do I generate the value to sum/multiply?


texts: From the IIR filter flow graph below i don't understand how the transfer function is calculated in every node:


texts: Is blending of colour lights a linear system?

texts:
Can we consider recombination process of different light colors a linear system/phenomenon? Like, when you shed two different light sources for example, a blue light and a yellow light on a surface or recombine them through a prism , the resultant light would be green which has a different frequency that the inputs. Producing a new frequency component other than what were in input (Is this assumption correct?) This is contrary to definition of linear systems I assume. So is this process linear, if not, how is it non-linear, I mean the whats the physics?
I've wondered this question from time to time and I know I probably have some wrong assumption some where. I would appreciate if someone clarify things for me or point out my wrong supposition(s).


texts: Applying 2D Image Convolution in Frequency Domain with Replicate Border Conditions in MATLAB

texts:
I have created a function that filters an image (250x250) with a gaussian blur kernel (5x5) using FFT and IFFT. I am trying to get my filtered image to equal exactly the filtered image created by the 'imfilter' function using the 'replicate' option. There are small differences in the two filtered images. They look identical visually but upon inspecting the pixel values, they are slightly different.
This is the general logic my code follows:
Zero pad the borders of the input image to size 254x254
Zero pad the kernel to size 254x254
'clamp to edge' the edge pixels of the input image. The outermost input image matrix rows and columns are extended to the borders, replacing the zeros. I am doing this to emulate the 'replicate' option in 'imfilter'.
Use 'fft2' to convert both the input image and kernel to frequency domain
Element wise multiply the two together
Convert the result to the time domain using 'ifft2'.
Crop away the border pixels, converting the filtered image from 254x254 ---> 250x250
What am I doing wrong here? Thanks.


texts: Amplitude Matching for Exponential Swept Sine

texts:
I am working in the area of aerospace vibration testing and we use swept sine tests on structures and measure the response using the accelerometers. I tried implementing the technique "Simultaneous measurements of impulse response and distortion with a swept sine technique" for analysis of the response acquired for a swept sine test. In the analysis we basically try to plot the accelerometer amplitude response versus frequency. I could plot the straight spectrum. However, a change in the amplitude in the accelerometer response causes a change in the amplitude of the plot.
Is there any way for amplitude scaling in the mentioned technique.  i.e if I increase the amplitude of the swept sine to some X value, then the final FFT spectrum of impulse response also should have the same X amplitude.
Please let me know.


texts:
What are the definitions of symbols and samples related to PCM? For example in binary encoding $\{0,1\}$ with 8 levels of quantization, is one symbol either of $ \ 0 \  or \  1\ $? or is it the collection of  $(x,y,z)$ where $x,y,z \ \epsilon \ \{0,1\}$.

Similary for ternary encoding $\{0,1,-1\}$,suppose with 9 levels of quantization, is a symbol either of $-1,1,or \ 0 \ $ or is it complete $(x,y)$ where $x,y \ \epsilon
 \ \{-1,0,1\} \ $ ?


texts: Construct controller with given crossover frequency

texts:
I have this problem:
I can solve question a) by looking for the phase gain of the bode plot at $\omega_{\pi}$, find the phase gain = 0.2, and divide $\frac{1}{|G(w)|} = \frac{1}{0.2} = 5$, so K can be at max 5.
How do I solve question b) though? I don't understand what is asked. What does $\omega_{c,d}$ stands for?


texts:
I am looking at implementing different window functions in an application I am making.  I have found the formulas for all major symmetrical windows, however I would like to implement half windows.

This may be a stupid question, but does the window literally jump straight from 0 to 1? or is there a small curvature to the start. ie, do I use the symmetrical formula and just use the second half? Or are there specific formulas for calculating half windows?


texts: Pre-Processing Wi-Fi Channel State Information (CSI) Data

texts:
I was successfully able to collect some CSI data using the existing tool(s) on GitHub (https://github.com/StevenMHernandez/ESP32-CSI-Tool). The CSI data is a pair of imaginary and real number which can use to extract amplitude and phase information. However,  I have no prior knowledge of digital signal processing and I'm not sure how can to preprocess the CSI data to remove environmental noise like the Butterworth filter and/or PCA. I'm using Python and would like to know if this is an appropriate tool.
I'm trying to recreate steps from the following papers related to human activity detection using Wi-Fi CSI:

SEARE: A System for Exercise Activity Recognition and Quality Evaluation Based on Green Sensing URL: https://ieeexplore.ieee.org/document/8247258
WiFi Sensing with Channel State Information: A Survey URL: https://dl.acm.org/doi/10.1145/3310194


texts: Why is the frequency measured in radians in filter design?

texts:
This is my first time posting a question, so sorry if my articulation is not proper.

I'v just started out with FIR filter design and I see the domain being -$\pi$ to $\pi$. Now I understand that discrete frequencies repeat after an interval of 2$\pi$ but consider this. I need to design a low pass FIR filter with a cutoff frequency of 100$Hz$. How can I specify this in terms of radians? Is it using $\omega$=$2$$\pi$$f$? But that would be outside the domain. As you can see, I am very confused about this concept so any help would be appreciated.


texts: What is the difference between natural response and zero input response?

texts:
I am new to DSP and was going through different responses of a system subjected to an input. My understanding of zero input response is: it is the response/output of the system when the input signal is set to zero. In other words if a system is described by a linear constant coefficient difference equation the zero input response is the homogeneous solution.

However if the $\mathcal Z$-transform of the input is a rational function $X(z)=N(z)/Q(z)$ and that of the LTI system function is $H(z)=B(z)/A(z)$ and the system is initially relaxed, then $Y(z)= H(z)X(z) = N(z)B(z)/A(z)Q(z)$. Assuming distinct zeros(real  only) and poles(real only) of $X(z)$ and $H(z)$ then

$$Y(z) = \sum_{k=1}^N \frac{A_k}{1-p_kz^{-1}} + \sum_{k=1}^L \frac{Q_k}{1-q_kz^{-1}}$$

which gives

$$y(n) = \sum_{k=1}^N A_k(p_k)^{n}u(n) + \sum_{k=1}^L Q_k(q_k)^{n}u(n)$$

where $p_k$ and $q_k$ are the poles of system $H(z)$ and input signal $X(z)$ respectively and $u(n)$ is the unit step function. Now the first term is referred to as the natural response of the system $H(z)$. It's very confusing to grasp the difference between zero input and natural response.

Edit: The reference of the question is to book DSP : Principles , Algorithms and Applications by John Proakis and D Manolakis
pdf of the book is here
Page no 203 and 204. The two paragraphs after formula 3.6.4 explains the difference between zero input response and natural response

Thank you Peter and Matt for your answers and comments.


texts: fft function in R vs spec function from 'seewave' package? They don't give a similar frequency spectra

texts:
I understand that spec will give me the frequency and the corresponding amplitude of that component, whereas fft will compute the DFT of the signal and throw the complex numbers for each component which I have to plot as per the physical frequencies.

y is the time-series signal
Fs is the sampling freq.

N = length(y);
f1=fft(y)
fa = seq(0,N-1,by=1)*(Fs/N)
plot(fa,abs(f1),type='l',xlim = c(0,3),ylim = c(0,50))
s1=spec(y,f=1000,plot = FALSE)
s1$x<-s1$x*1000   (KHz to Hz)
plot(s1$x,s1$y,type='l',xlim = c(0,3))

Red one is from spec function and blue is from fft function.


texts:
I'm trying to visualise and understand the output of constant q transform. Based on my experience with FFT and the articles I've read so far, I would expect that each output bin would represent one frequency band (logarithmically spaced). Which means a naive way o finding the fundamental frequency would be to find the bin with highest energy in the spectrum.

The constant q algorithm I'm using has complex output, and I assume that the real part is the magnitude and the imaginary part is the phase. So, I compute the spectrum, but odd enough, the result looks exactly the same for every note or chord I play on the piano:

Which pretty much debunks my initial idea of finding the highest energy bin (it's always the last one).

What am I doing wrong? Is this really how the output should look?

EDIT:

Here is the code I use:

chromaSize = 24
frameSize = 8192 * 4 * (chromaSize / 12)

fft = es.FFT(size = frameSize)
cq = es.ConstantQ(binsPerOctave = chromaSize)
ip = es.InstantPower()

audio = pyaudio.PyAudio()
stream = audio.open(format=pyaudio.paFloat32, channels=1,
                rate=44100, input=True,
                frames_per_buffer=frameSize)

plt.ion()
plt.show(block = False)

while True:
    frame = np.fromstring(stream.read(frameSize, False), dtype=np.float32)

    if len(frame) >= 0 and ip(frame) > 0.0001:
        f = fft(frame)[1:]
        r = np.abs(cq(f))

        plt.clf()
        ind = np.arange(len(r))
        plt.bar(ind, r, color='blue', alpha=0.5)
        plt.draw()
        plt.pause(0.0001)

I'm using essentia's implementation. Ultimately, my goal would be to extract the fundamental and it's harmonics, which will then get compared using markov chains with a given data set.


texts: When can we write Heisenberg uncertainty Principle as a equality?

texts:
We know that Heisenberg uncertainty Principle states that
$$\Delta f \Delta t \geq \frac{1}{4 \pi}.$$

But (in many case for Morlet wavelet) I have seen that they changed the inequality to an equality. Now my question is when are we allowed to change the inequality to an equality:
$$\Delta f \Delta t = \frac{1}{4 \pi} $$
why =


texts: buidling a bandpass filter without using inbuilt functions

texts:
I got an assignment which require building a bandpass filter without using any inbuilt functions. I wanted to create a third order Butterworth bandpass filter

%for ecg
f1=5; %cuttoff low frequency to get rid of left over baseline wander
f2=15; %cuttoff frequency to discard high frequency noise
Wn=[f1 f2]*2/fs; % cut off based on fs
N = 3; % order of 3 less processing
[a,b] = butter(N,Wn); %bandpass filtering
f = filtfilt(a,b,ecg);
s1= f/ max( abs(f));

I wanted to get the coefficient without using butter function in MATLAB.


texts:
I am trying to compute the Fourier series representation of a full wave rectifier output.

The equation of the signal is:

$ x_8(t) = | \cos (2 \pi f_o t) $ |

I have tried to find the Fourier series representation of the above signal but I ran with the following difficulties:

1) The solution given in the book is the following which is different from my answer:

2) Surprisingly, my answer matches with the solution of the next problem given in the book which is:

Question: $ x_8(t) = | \cos (2 \pi f_o t)  | + \cos (2 \pi f_o t) $ (Half Wave Rectifier Output)

Solution:

I am not able to understand whether what I have computed is correct or not. If it is incorrect what is the mistake that I have made? I am enclosing a copy of my answer sheet.

My Solution:

Thanks and regards.


texts: How to design a filter with custom group delay -- why is it an IIR?

texts:
Say I have an FIR filter with desired magnitude characteristics. I could imagine using an allpass IIR in series with the FIR to alter the group delay characteristic to provide both a 'magnitude equalization' and 'delay equalization' curve custom for the application.

I imagine this would be simpler in the frequency domain where the imaginary component of each slice of the FFT would adjust the group delay of a bandwidth corresponding to the nyquist frequency of the system divided by the FFT length, and then IFFT'ing to find the impulse response, windowing that impulse response to design the filter.

I see Matlab has an IIR feature for custom group delay in fdatool. Why does it have to be an IIR? Can an FIR have a custom group delay?

Thanks for any pointers!


texts: Degradation of BPSK system using matched filter with a timing error

texts:
Given a BPSK system with:
\begin{align*}
s_1(t) &= Au_T(t) \\
s_2(t) &= -Au_T(t) \\
\end{align*}
where $u_T(t)=1$ for $0<t<T$, and zero otherwise. Given $A > 0$, and Gaussian white noise of power spectral density $\frac{N_0}{2}$.

Design a matched filter for this system.
Draw to scale the output of the filter due to signal only.
Assuming the system has a timing error of $\varepsilon$ with respect to the optimal sampling instant, find an expression for the probability of error.

My attempt:
I have solve the first part and find out $h(t) = \frac{1}{\sqrt{T}}$. However I have no idea how to do the rest.


texts: Shape of structuring elements for morphological gradients

texts:
I am looking to understand recommended shapes of structuring elements used in calculating morphological gradients. According to Pierre Soille: Morphological Image Analysis:

  Only symmetric structuring elements containing their origin are considered. By doing so, we make sure that the arithmetic difference is always nonnegative.

The arithmetic difference mentioned in the quote is referring to three combinations currently used to calculate the discrete gradient:

arithmetic difference between the dilation and the erosion;
arithmetic difference between the dilation and the original image;
arithmetic difference between the original image and its erosion.

But, I think using a SE containing its origin is enough (it ensures anti-extensivity of dilation and extensivity of erosion). In this case, the following holds and ensures nonnegativity in all three cases:

$\varepsilon_B \leq id \leq \delta_B$ (where $id$ is the identity transform)

I am looking for a reason to enforce the symmetry condition. Intuitively, I understand that using a symmetrical SE is better than using a non-symmetrical one (e.g. examining a symmetrical pixel neighborhood). It was also suggested to me that there might be historical reason for this constraint.

However, I would like specific examples, arguments or references that point to desirable properties of symmetrical SEs (or undesirable properties of non-symmetrical ones).


texts: Having problems with understanding the meaning of PSD for colored Gaussian noise

texts:
After reading two articles on signal processing stack exchange:
On coloured Gaussian noise
How the white and colored noise differ in time domain
I do understand that variance do not change over frequency (or time) for colored Gaussian noise.
Let's say I have a random process $y_n$, which is defined as:
$$ y(n) = x(n-1) - x(n), \quad\text{where }x_n\text{~} N(0, \sigma_n^2)
$$
For the clarification, $x_n$ is white Gaussian, and hence $E[x_n\cdot x_m]=0 \text{ for } n\neq m$.
Here, we can see that PSD of $y_n$ has high-pass form since
$R_{yy}(\tau)=2R_{xx}(\tau)-R_{xx}(\tau-1)-R_{xx}(\tau+1).$
we can also easily see that its variance is $2\sigma_n^2$.
Next, we also introduce another random process which has a 2x higher variance than $x_n$:   $w_n\text{~}N(0,2\sigma_n^2)$
Now, I would like to compare the random process $y_n$ to $w_n$.
First, we know that their variance are identical to $2\sigma_n^2$.
Second, although $y_n$ is colored Gaussian and $w_n$ is AWGN, their PDF must look the same, because they are still Gaussian process.
What I do not understand clearly is the meaning of PSD of $y_n$. As shown below, plot with a red color is a PSD of $y_n$. It shows high-pass shape, start having higher power magnitude after passing $0.5 f_s$. But what does it mean for a Gaussian noise having higher power magnitude at high frequency and having lower power magnitude at low frequency?
It cannot be that $y_n$ has higher "variance" at high frequency and lower "variance" at low frequency because I mentioned earlier that variance should not vary over frequency or time.
The next thing I thought about is the difference between certain sequences of $y_n$ over time. For instance, if $y_n$ has a high frequency pattern (i.e., transition density of 1) for $n=1\cdots 100$, then it may temporarily have a higher variance for this particular high frequency pattern. But soon I realized that this argument does not make sense neither since Fourier transform does not care about behavior of $y_n$ for some local pattern.
If anyone can give me any good insights about the meaning of PSD for the example above, I would be greatly appreciated!
Please let me know if anything is not clearly in the above. English is not my first language so there is a very high chance I did not explain things very clearly.
$Y(f)$ vs $W(f)$" />


texts: Tracking object that is missing in some frames

texts:
I am trying to track small object in videos using LK
but the problem is, some time the object does not appear in the frames !

How can I track an object that does not appear in few frames ??
 how can I approximate it is place automatically ,,,

I am trying to find the object using Optical flow LK ... but this technique lost when my   object does not appear

Any suggestion ideas ???


texts: Difference between "Discrete Wavelet Transform" and "Discrete Wavelet Decomposition"

texts: I have a rough overview on Discrete Wavelet Transform (DWT). However, I am confused about Discrete Wavelet Decomposition and did not find a good reference yet which explains this well. What is it actually about? Is it somehow part of DWT or an inverse operation to it?


texts:
In ECE-700 Filterbank Notes the partitioning is $ \Delta = \frac{2\pi}{M}$ with center-frequencies at $ \omega_k = \frac{2\pi}{M}k$ which makes me wonder, considering that frequencies in the range from $\pi$ to $2\pi$ are just the negative versions of the range from $0$ to $\pi$, which divide into  low- and high-pass frequencies.

In the follow-up picture the partitioning then stretches from $0$ to $\pi$ - is there someone so nice as to explain to me what I'm not understanding here?


texts: How do you interpolate between points in an image (2D), e.g. using splines?

texts: I can understand just fine how to use 1-dimensional interpolation on data points where one coordinate is a function of the other: y = f(x). However, when we have an image, we generally have points in the image that are not given uniquely by such a function, e.g. two points (x1,y1) and (x1,y2), so how do we describe this in order to interpolate between points in the image?


texts: What are the effects on filtering upsampled data without interpolation?

texts:
I have an application that displays several signals from different sensors at various sample rates.  In order to display the data, I "stretch" the signal by repeating samples to match the highest sample rate.  For example: if I display a 500 Hz signal and a 100 Hz signal at the same time, every sample in the 100 Hz signal will be repeated 5 times so it appears to have a 500 Hz sample rate.

When I apply filters to the signals (all second-order IIRs: LPF, HPF, notch) I use filters designed for the highest sample rate.  In my example here, I'm using 500 Hz filters on an upsampled 100 Hz signal.

How "bad" is this?  Is my only option to design an interpolation filter during upsampling, and if so, can anyone lead me to good references on how to design such a filter?


texts: I and Q components and the difference between QPSK and 4QAM

texts:
Both 4QAM and QPSK apparently produce the same waveform, but are they the same mathematically?

In a QPSK constellation, are the mapping points at 45, 135, 225 and 315 degrees while the 4QAM is at 0, 90, 180 and 270?

I also struggle to understand the I/Q components of such a constellation diagram. What does "inphase" and "quadrature-phase" actually mean? Are they just another way to specify the real and imaginary part for this type of usage?


texts: Kalman Filtering and space parametrization

texts:
I am familiar with Kalman filtering given a linear (time-invariant) state space model.
However, the state space parametrization is not unique. Given a controllable and observable state space model (A,B,C,D), I can find multiple equivalent forms, e.g. controllable canonical form, observable canonical form, Jordan normal form, etc.

Is there a state space parametrization which is beneficial to transform the model to before running the Kalman filter, especially in terms of numerical stability.

Are there papers you could point to?
I am running a Kalman Filter on an 8th order SSM on a DSP and am running into numerical issues.

Thanks a lot!


texts: How to enhance Pedestrian Dead Reckoning using Particle filter?

texts:
I am working on Pedestrian dead Reckoning using Particle filtering. I am using mobile accelerometer sensor to detect number of steps and stride length and calculating the total distance traveled.

I want to use particle filter to enhance accuracy in my results. I have read the theory of particle filter but stuck at the point that how to apply the model to my case.

In particular my problem, I have the following questions:

How to determine the initial position pdf distribution? I have a Bluetooth beacon to detect the starting position. Particle filter requires a cluster of point at the start, right? I don't think in my case the distribution is Gaussian. Then what should it be?
Should I do the particle filtering on the smartphone or pass the data back and do it on the computer using Matlab? Will the computation very intense.


texts: Can I get some help plotting an IIR filter frequency response in Python?

texts:
The problem I'm trying to solve is described here.

In my script, I used Euler's identity to calculate z^-n, such that z^-n = cos(omegan) - jsin(omegan). I find the magnitude of these values from a range of 0 to pi, and then plot them. However, the resulting graph aperiodically oscillates between 0 and ~ -8 dB, and I'm honestly not sure why.

I'd appreciate any help! The code is as follows:

import math
import matplotlib.pyplot as plt

# define constants
b0 = .125 a1 = .875 x = [] total = 0

while total <= math.pi: x.append(total) total += math.pi/256

return the magnitude of z-n

def find_z(n, omega):

# z = re^jw # e^jw = cos(w) + jsin(w) # z^-n = cos(wn) - jsin(wn)
ret = abs(math.cos(omega * n) - (1j * math.sin(omega * n))) return ret

h = []
for i in x: h.append(b0 / (1 - (a1 * find_z(-1, i))))

for i in range(len(h)): h[i] = 20 * math.log10(h[i])  plt.plot(x, h, color='blue', marker='o', linestyle='-',  linewidth=2, markersize=12)
plt.xscale("log") fig = plt.figure(figsize=(12, 6)) plt.show()

print(h)


texts: Resonant peak frequencies for phaser effect

texts:
I have a conventional phaser setup with allpass filters and feedback.

I can calculate the peaks frequencies of the phaser when the allpass delays are set to 1 (then they work as simple one sample delays). Adding two allpass delays will give an additional peak, and the peak frequencies for peaks i = 1 to N is found by:

(i / (N + 0.5) ) * SampleRate / 2

where N is the number of peaks (and the number of allpass-delays divided by 2).
So if you have a 6-stage phaser at 44100, the 3 peaks will be at:

(1 / (3 + 0.5) ) * 44100 / 2 = 6300
(2 / (3 + 0.5) ) * 44100 / 2 = 12600
(3 / (3 + 0.5) ) * 44100 / 2 = 18900

However, these frequencies start warping when the allpass-delays are not set to 1:

How can I calculate the peak frequencies when the allpass-delays are not set to 1?


texts: Equalize Speaker using Frequency Response Curve

texts:
I have a speaker for which I know the frequency response curve (-5dB at 600Hz, -3dB at 700Hz, ... +5dB at 10kHz). This response curve is not as flat as I would like, and so I want to write a script which takes an audio file as input, and then transforms it such that when this new audio is played out through the speaker, it sounds like what the original file would sound like if it were played through a speaker with an actually flat response curve. (Basically, I want to flatten the frequency response of this speaker in software.)

What would be the best way of going about this? My first instinct is to Fourier transform the audio signal, reweight each Fourier coefficient according to the inverse of the amplitude at the corresponding frequency in the frequency response curve, and then inverse Fourier transform the result. However I imagine that such a direct approach would result in problems such as ringing artifacts, and I know that there are ways of performing similar tasks through the use of specially designed filters.


texts:
In case of using equalizers with ZigBee systems, what type of fading does the channel show, flat or frequency selective, fast or slow?
Is the inter-symbol interference (ISI) problem encountered in ZigBee systems?

I am asking this question because I konw that ZigBee Systems transmit data at low rate (250 kb/sec), and ISI problem and hence frequency selective fading appear with high data rate applications. I found some papers employing equalizers with ZigBee under Rayleigh fading conditions, but the type of fading was not indicated.


texts: How do multirate systems introduce 'efficiency' in computation?

texts:
I've started reading more into multirate systems recently and am somewhat confused. The concepts introduced seem to make sense, that is, how the decimation and interpolation can change certain characteristics of signals (e.g. aliasing and imaging in the frequency domain, etc).

However, I keep hearing the term 'computational efficiency' when subject of multirate systems is introduced. What confuses me is how exactly do we gain this efficiency? What are some applications where implementing multirate introduces efficiency? From a quick glance, the need to downsample and upsample requires more computations in the form of anti-aliasing and anti-imaging. Does the ability to operate at lower sampling rates always counteract this added cost?

Say for example I have a signal that was sampled at 44.1kHz. Let's also say that signal's spectral content is primarily based from 20-2000Hz and there is a disturbance noise at 3000Hz that we would like to remove. Now a single rate system would essentially just apply a digital filter that is capable of removing that 3000Hz disturbance. Let's assume we are using an FIR filter. This filter would have to operate at a 44.1kHz update rate which is pretty fast. In order to create a filter that have an appropriate cutoff frequency of say 2500Hz, the filter length would also have to be pretty long due to the fast update rate.

An alternative method would be to downsample the signal to get to rate where the frequency spectrum ranges from 0 - 3000 Hz (instead of 0 - 22.05kHz). This frequency ranges seems more appropriate for the operations we're trying to do on our signal. If we were to slow down to that rate, the FIR filter length would be smaller(could I get a confirmation on this statement?) and the filter would operate at slower rate. We would essentially downsample, modify the signal, and then upsample.

However, in order to not pollute our lower sample rate signal, it's necessary to introduce anti-aliasing and anti-imaging for the downsample and upsample operations respectively. I realize there are efficient methods for doing integrating the filters with the downsampler/upsampler (i.e. polyphase realization), but how do I know whether these added computations are worth the cost?

TLDR: How does a multirate system introduce computation efficiencies when it's necessary to add more computations when upsampling and downsampling?


texts: How to design a filter with a certain magnitude response

texts:
I am trying to design a filter whose magnitude is the same as that of a given signal. The given signal is wind turbine noise, so it has significant low-frequency content. After designing the filter, I want to filter white Gaussian noise so as to create a model of wind turbine noise. The two signals, that is the original and the filtered noise should sound similar.

I am using arbitrary magnitude filter design in Matlab for that (FIR, Order: 900, Single rate, 1-band, response specified by amplitudes, Sample rate 44100 Hz, Design method: firls). The problem is that, although I design the filter using the values from the original signal's magnitude, the filter magnitude fails to follow the magnitude at higher frequencies. Could you please help me with that?

The idea is that I am using a polynomial curve, to fit the shape of the spectrum of the original sound. Then, I filter the magnitude of the generated noise, by multiplying the extracted polynomial curve with the generated noise magnitude. Finally, after calculating the new magnitude and phase, I get back to the time domain.

[x,fs] = audioread('cotton_0115-0145.wav');  % original noise sample
x = x(:,1);                                  % extract one channel
x = x.';
N = length (x);

% fft of original signal
Y = fft(x,N)/N;
fy = (0:N-1)'*fs/N;

% half-bandwidth selection for original signal
mag = abs(Y(1:N/2));
fmag = fy(1:N/2);

% polynomial fitting
degreesOfFreedom=40;
tempMag=20*log10(mag)';  % desired magnitude in dB
tempFmag=fmag;

figure(1)
plot(tempFmag,tempMag,'b');
title('Spectrum of original signal-Polynomial fitting')
xlabel('Frequency (Hz)');
ylabel('20log10(abs(fft(x)))');
axis([tempFmag(1) tempFmag(end) min(tempMag) 0]);

hold on
p = [];

for i=1:4
    p=polyfit(tempFmag,tempMag,degreesOfFreedom);
    pmag=polyval(p,tempFmag);

    plot(tempFmag,pmag,'r');
    pause

    above=pmag<tempMag;
    abovemag=tempMag(above);
    abovefmag=tempFmag(above);

    tempMag=abovemag;
    tempFmag=abovefmag;

end

hold on
legend('signal magnitude','polynomial');

%
loc1 = find(fmag == 0);
loc2 = find(fmag == 22050);
Nmag = length(mag);

M=((Nmag-1)*max(tempFmag))/(tempFmag(end)-tempFmag(1));
freqFinal=tempFmag(1):max(tempFmag)/M:max(tempFmag);
freqFinal=tempFmag(1):max(tempFmag)/length(mag):max(tempFmag);
magnitudesFinal=polyval(p,freqFinal);

figure(2)
plot(fmag,20*log10(mag)');
hold on;
plot(freqFinal,magnitudesFinal,'g');
title('Spectrum of original signal-Choice of polynomial curve')
xlabel('Frequency (Hz)');
ylabel('abs(fft(x))');
axis([freqFinal(1) freqFinal(end) min(magnitudesFinal) max(magnitudesFinal)]);

%%
% noise generation
Nn  = N;
noise =  wgn(1,Nn,0);
noise = noise/(max(abs(noise)));
Ynoise = fft(noise,Nn)/Nn;
fn = (0:Nn-1)'*fs/Nn;

% polynomial for whole f range
newmagA = 10.^(magnitudesFinal/20);
newmagB = fliplr(newmagA);
newmagB(end+1) = newmagB(end);
newmag = [newmagA newmagB];

%filtering
Ynoisenew = newmag .* Ynoise;
figure(3)
magnoise = 20*log10(abs(Ynoisenew(1:Nn/2)));
fnoise = fn(1:Nn/2);
plot(fnoise, magnoise);

% magnitude and phase of filtered noise
magn = abs(Ynoisenew);
phn = unwrap(angle(Ynoisenew));

% Back to Time domain
sig_new=real(ifft(magn.*exp(1i*phn)));

figure(4)
sig_new = sig_new / max(abs(sig_new));
plot(t, sig_new);

Ysignew = fft(sig_new,Nn)/Nn;
fn = (0:Nn-1)'*fs/Nn;
figure(5); plot(fn(1:Nn/2),20*log10(abs(Ysignew(1:Nn/2))));


texts: Removing a darkened caption with text on it in an image

texts:
Image processing  is new for me, so forgive me if my question sounds vague.

I'll try to supplement it with what I have learned until now and also a couple of images.

See the caption in the image below?

What I'm essentially trying to do, it to remove it and restore the original image in (Python using OpenCV).

Now I have a couple of approaches in mind. First one I read about is a technique called Inpainting. Now I saw a tutorial on inpaiting here but this required me to create a separate mask where the non-zero pixels denote the stuff I want gone.

Now what I noticed is that the caption is not fully opaque. So i was wondering if there is any possible way to restore the original image by first removing the darkened part of the strip. (essentially something very the original image with only the whitened text on it) and then create a mask of the text and then use inpaiting.

Now I have a couple of questions.

What technique do I use to remove the darkened part (let the text be now, we can remove it in the second step using inpaiting)
Does this algorithm even make sense. Is there a better approach I should be looking at?

NOTE: In no way am I looking for any sort of code or specific implementation. I'm just looking for what techniques and procedures I can study up on so as to get the job done. The rest is on me


texts: When lowpass filtering with a Gaussian kernel, how can $\sigma$ be related to cycles/pixel?

texts:
I wish to low pass filter an image with a Gaussian kernel. However, I also wish to select a filter that has the property of producing images that have have no greater than N cycles/pixel. Let's say I wish N to be close to, or exactly .025.

By what method is $\sigma$ selected in order to obtain such an image?


texts: Beginner's book in signal processing with practical examples on fault detection in electrical motors

texts:
I'm an electrical engineering graduate and I'm working on condition monitoring for electrical motors and generators. My goal is to improve reliability of electric rotary machines in such way, that any ongoing faults are discovered before they progress in major failures. Main role in this domain has signal processing (and also statistical analysis along with machine learning), which I would use to develop condition monitoring algorithms. Unfortunately signal processing wasn't included in my education program, so I have to learn this on my own.

So what I think that I need, are some sources which outline  basic concepts, typical tools (FFT, STFT, wavelets, filters...) and problem-solving examples (noise removal, hidden pattern extraction, ripple reduction,...)  with only so much math involved that the main idea can be grasped, and that implementation is possible, which would be demonstrated in lots of MATLAB examples in the context of electrical rotary machines. I truly don't want a source for signal processing, which is aimed for mathematicians or information engineers, but rather sth applicable for electrical engineers.


texts: How to get some law parameters from its FFT figure?

texts:
Given an hidden law $$I(t)=\frac{1}{2}[A+B.cos(2\pi\nu_0t+\phi_0)]$$ I observe its signal with 256 samples and compute its FFT over this interval, and  represent the FFT result from 0 to 1 (normalised frequency).

I get a figure whose pic is at 0.83.

How can I retreive A,B and $\phi_0$ from this graphics ?

A second question is : is it normal that this figure is not symmetric ?

A final question is : doing an analytical computation on classical Fourier transform gives a result. How much does what I observe from FFT and the analytical result differ ?


texts: Extracting audio from picture of time-domain signal

texts:
I saw this nifty art thing where they dress up the time-domain audio signal in some pretty way.
https://www.etsy.com/ArtbloxShop/listing/621561709/sound-wave-art-custom-soundwave-print

My question is whether I can take a photograph of that object and somehow reconstruct the actual audio?  I'm assuming that the printing and my photograph won't have the resolution to directly capture the wave data points.  But I was wondering if using some basic assumptions about frequency if we could reconstruct something crude.  I suspect that too much information is lost, but I was hoping there might be some clever tricks.

There is clearly some resolution of printing that would preserve enough information, but I'm wondering if we're talking about a poster-sized version of a short clip or something.  That wouldn't be very practical.  It would be fun to have time-domain audio hanging as art if I could take my smart phone and take a picture of it, and after some processing, play back the audio.

As an alternative, could this be done using a spectrogram?  Like if that block of art displayed that instead?  Or maybe some combination of that plus time-domain envelope?  If I did want to represent audio as an image as part of art (maybe presume something around 8.5x11 size), is there some way I could have an app photograph it and turn it back to audio?


texts: Real-time DAQ (Data AcQuisition) device and filtering of voltage signals in audio frequencies

texts:
I'm working on a frequency mixer project in which two analog voltage signals at two different frequencies (each <5 kHz) are considered. I'd like to add the two fundamental frequencies and return an analog output in real time. Can I do this digitally and how difficult it is compared to an analog approach?

My idea so far includes an ADC with a fast enough sampling rate, Python scipy (for filtering), and DAC for output.

I'm trying to go low cost (<100 bucks). As far as I'm concerned, I don't see a potential problem in this idea as an ADC with Msps can be obtained, the math process is fast enough, and the DAC can also have Msps sampling rate.

But as a chemistry student, I'm no professional at DSP, so please point out some of the potential flaws/drawbacks of this method and tell me if you'd recommend doing this in analog?

Thanks in advance.


texts: Design $ {L}_{2} $ Norm Optimal Infinite Impulse Response (IIR) Filters

texts:
It is widely known that matching a FIR filter of fixed length to a band model is an unconstrained QP-problem. The MATLAB function firls() implements a solution to this problem. Basically, one symbolically evaluates integrals of the (Z-domain) transfer function over the unit circle, getting astonishingly simple expressions for the match quality. Additionally, that method is stable and will not fail, other than the Parks McClellan Filter Design Algorithm.

Is it possible to fit IIR filters to such a model in the $ {L}_{2} $ Norm Sense?
Several problems arise:

While the cost function (fit quality?) can be written down, is it not a simple quadratic anymore, but a fraction of two quadractics.
You appear to have to optimize nominator and denominator at the same time.
It might be hard to constrain the denominator to "minimum phase", so that the resulting filter is stable.

I am not looking for gradient descent or other local optimization methods, as the cost function might have several local optima. An idea that crossed my mind would be that it appears possible to pose the FIR problem so that it yields the optimal coefficients for an IIR all-pass.


texts: Weird behaviour converting velocity to displacement

texts:
I'm trying to convert accelerometer data from accelerations to displacement,

I assume that the initial velocity is zero and the initial position is also zero.

Here's my code,

Fs = 500;
T = 1/Fs;
t = (0:length(y) -1)*T;

vel = cumtrapz(t, acc);
disp = cumtrapz(t, vel);

here are the respective plots from acceleration to displacement,

what am I doing wrong?

I've implemented a kalman filter but the result is the same. I'm trying to work around with wavelet decomposition, any help would be very appreciated to resolve this problem.


texts: How can I correctly plot an impulse_response() of a discrete transfer function?

texts:
I posted this question on stack overflow and I think it was the wrong place so I moved it here.

Hello, I'm new to python, I just want to use the control library and I'm slowly learning, I did my research and read the docs but as a newbie programmer I don't really understand much. I know it's about the syntax but I don't know what exactly is wrong. I'm moving from matlab to python and I can't get this to work:

from control import *
from sympy import *
from numpy import *
from matplotlib.pyplot import *

Gp = tf(100,[1,0,100])
T = 0.05
Gz =sample_system(Gp, T)

x = linspace(0.0, 0.7, num=15)

T, yout = impulse_response(Gz, T=x)
plot(T, yout)

if I run this I get this output:

ValueError: Time steps ``T`` must match sampling time

and if I change x = linspace(0.0, 0.7, num=15) for x = linspace(0.05, 0.7, num=14) I get this:

ValueError: A value in x_new is below the interpolation range.

No idea what I'm doing wrong, I read the docs on impulse_response() and linspace(), I can't find any examples of similar problems or people plotting impulse responses using python. Thanks in advance

Edit1: Okay I've gotten here so far:

imp = signal.unit_impulse(100, 'mid')
T, yout = signal.butter(4, 0.2)
z = array([0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5])
T, yout = impulse_response(Gz, T=z, input=imp)
T, yout = signal.butter(4, 0.2)
plot(T, yout)
grid(True)
show()

now I get a plot that looks like this:

But I need something like this:

Edit2: Okay seems like the whole problem was that yout had these dimensions: (1,15) and all I had to do was .flatten(), on the other hand I was using wrong linspace so I changed it for arange, and the end result is this:

Gp = tf(100,[1,0,100])
ts = 0.05
Gz =sample_system(Gp, ts)
pprint(Gz)

t = arange(0.0, 15*ts, ts)

T, yout = impulse_response(Gz, T=t)
yout=yout.flatten()
plot(T, yout)
margins(0.1, 0.1)
xlabel('Time [sec]')
ylabel('Amplitude')
grid(True)
show()

and the output is:

Still seems to me that something is missing, but I don't know what, and that this is not a clean solution, anyone got a cleaner way of doing this? Having to flatten() seems wrong to me... in matlab you just had to impulse(Gz, t) and got the plot right away.


texts: What does orthogonality mean in the context of adjacent channels in a filter bank?

texts:
First off, what does orthogonality mean in the context of adjacent channels in a filter bank?

Secondly, I am reviewing some code in which the dot product of the transfer functions of adjacent channels are being used to "test orthogonality". Is this a formal method with a name? How does it work?


texts:
I know the size of DFT is the number of samples used to do the calculation.
But the audio signal (and other signals) normally have many samples.
So, what does really happen, if I choose $2$ seconds of audio, which will have $2 \times 44100$ samples, and $1024$ as the size of DFT ?

First, the result of the DFT is the amplitude and phase of the frequencies which are stepped by $44110 / 2 / (1024 / 2) \approx 43\textrm{Hz}$.
So we get $0\textrm{Hz, } 43\textrm{Hz, } 86\textrm{Hz, } .... \textrm{to } 22050\textrm{Hz}$.
Is this correct?

Second, as I analysis the $2$-second signal, the calculation is done by using every $1024$ samples to calculate $512$ DFT bins, then what ? We can get $44100 * 2 / 1024 = 86$ results of DFT, we calculate the mean of every result to build the final result ?


texts: What is Wavelet packet transform (WPT) time complexity

texts:
I have read in many papers that time complexity of WPT is $O(N \log(N))$, but I would like to know what is the effect of Levels in the time complexity.

Must it be one of the following ?

$O(2^{L} ~N~ \log(N))$
$O(2^{L+1} ~N~ \log(N))$
$O(L ~N ~\log(N))$

$ L$ is the number of decomposition levels

Can anyone helps which is true?


texts:
I'm trying to retrieve a melody from a signal (simple sine with 3 notes).

I've windowed the signal and performed autocorrelation to get the fundamental frequency of each window and then converted frequency to note.

So for a signal with notes "G A E" I am getting something like "G G G A A B E E E E F". It makes sense that notes are repeated and there are some mistakes... But how should I deal with this? If I use a bigger window then I might miss some notes with shorter duration. I can think of a really naive algorithm to get the melody for this example but what are some real techniques that are usually applied?


texts: Mathematical expression of Multicarrier frequency equalizer

texts:
I have a signal $X$ with length of $N$, multiplied with any unitary matrix, i.e  the transpose of DCT matrix as:
$x = D' X$
where $D'$ is the transpose of the of DCT matrix. Then let's add $CP$ guard interval into the signal $x$. The resultant signal after adding the cp into the head of the signal is $x_2$.
We transmit $x_2$ through a channel $h$ leading to have the signal $y$ as follows:
$y = h*x_2$   ......  where $*$ is the convolution operation.
Hint: the convolution is not circular in that case as we used DCT instead of DFT.
At the receiver, after removing the delay of the channel and the CP guard interval, I used the following steps:
$x_3 = ifft(fft(y)./H)$ .....  I used ZF equalizer, $H$ is the frequency-domain channel found by $H = fft(h,N)$
Then, the signal is recovered  by multiplying by the DCT matrix
$x_4 = D * x_3$
My question,  how does it work the step when I get $x_3$?? is there mathematical expression for that,  , I think it's called frequency domain equalizer, but I didn't get its mathematical expression. Most importantly, why does it work however we don't have circular convolution when getting $y$? !!
thank you


texts: Combining two images with different resolutions?

texts:
I am attempting to combine two radar images, each of which is a 3-dimensional cube with a range dimension, an azimuth angle dimension, and an elevation angle dimension. One image has high azimuth resolution but low elevation resolution, while the other has high elevation resolution but low azimuth resolution.

Is there an effective way to combine these images to create a single image, and what sort of resolution can I achieve with the new image?

Thank you in advance for any advice or resources.


texts: frequencies in frequency spectrum with no correlation together

texts:
I have a lack of understanding of the following questions:
If I have a signal from a motor that is recorded with an accelerometer. And the rotating speed of the motor is 150Hz(rpm 9000 ), I can see in the frequency spectrum harmonics of the fundamental frequency(150hz). The harmonics are integer multiples of the fundamental frequency.
But what are the other frequencies for example peaks at other frequencies that are not integer multiples of the fundamental? Shouldn't all frequencies that occur from a motor actually be something that is related to the rotating speed?(and i need to understand the building blocks of a motor to understand the frequencies to find a correlation to the non-integer harmonics?)
Can there also frequencies with no correlation to the fundamentals? That would mean i have multiple fundamental frequencies? and  there is something next to the roational speed of the motor, that also has his own frequency or what?


texts: Expectation of product of Filtered White Gaussian Noise

texts:
Assume I have a random process $X(t)$ that is white gaussian noise with psd $S(f)$.  Now, let's filter that noise through a LPF of bandwidth $B$.  How would I evaluate the following expression:

$$E[S(f_1)S(f_2)]$$


texts: Subtracting signals to re-construct the middle

texts:
Below image represents a sample of some of my captured data. There are three waves if you will. The middle wave has a very large dip ( cluttered by noise ) .

My goal here is to try and resurrect the middle wave to be noise free ( as much as possible ) . I will know I succeeded when the middle wave looks very 'similar' to the other two .

Here is what I am thinking, but I am not sure if the idea I have in mind is mathematically correct.

Given the red portion ( samples between 1920 to 2080 ) and the yellow portion ( samples between 2020 to 2170 ), I can take two FFT's. I can then determine the frequency that is dominant in both sample sets. If I take the average between these two ( FFT from red portion and FFT of yellow portion ), would that give me a mathematically correct estimation of what the middle wave should be?


texts: Are there any realtime voice activity detection (VAD) implementations available?

texts:
I'm looking for a realtime VAD implementation like Discord/Mumble/Teamspeak have.

My goal is to have a system that let's me use the same VAD for all my VoIP applications, as well as something I can use to facilitate VAD for push to talk only applications.

Mumble is open source though I wasn't able to find what I needed when I had a look.


texts: How do I find the phase response of a filter from the filter equation?

texts:
I am new to this forum. I have been teaching myself DSP for the last two years, so my experience and knowledge has a long way to go.

I am trying to find the phase response of a specific filter that I am using on Csound for any given fundamental frequency of the input (and I guess I'll need to be able to do that for any cuttoff frequency and ideally any sampling rate as well).

The filter is a 1st order lowpass IIR, and the filter's equation is as follows:

$$y[n] = c_1  x[n] + c_2 y[n-1]$$

$h_p$ - half power point
$b = 2 - \cos(2 \pi h_p/s_r)$
$c_2 = b - \sqrt{b^2 - 1}$
$c_1 = 1 - c_2$

Three questions:

Can I find the phase response for any frequency based on this equation alone?
Can I derive a phase response equation from this equation alone?
Is there a standardised way of finding the phase response for any lowpass filter, or does it depend on what kind of filter it is?

I have found some info about finding phase response both in these forums and elsewhere online, but without knowing the answer to question 3, it's hard to know if they are relevant or not.

I have also read most of the scientists and engineers guide to DSP (which I notice is also the boards DSP guide, awesome!), and I have been going back over it the last few days to try and get some clues. The best solution I can find is to fft the impulse response, and I know how to find the impulse response for a cuttoff frequency, but it'll take a long time to calculate the impulse response for all possible cuttoff frequencies.

EDIT:

I have been reexamining ths online resouce from analog.com:

http://www.analog.com/library/analogDialogue/archives/43-09/active_filters.html

The website tells me this about the phase response of a 1 pole lowpass filter:

Is this true of any 1st order low pass filter, or are they refering to a specific filter type? If this is true of my filter, I think I might be able to do this (though I won't say for sure yet!)


texts: Response function for digital camera pixel

texts:
To help me with doing some signal processing on an image, I'd like to model how a digital camera works.

The camera has a bunch of pixels.  What is the relationship between the value of that pixel in the resulting image, compared to the intensity of the light entering that pixel sensor?

Is it a linear function?
Is there some nonlinear response function?
If it is nonlinear, how  can I recover or model that response function?


texts: Is the phase of a digital filter always an odd function from $[-\pi, \pi]$

texts: Trying to prove something about the integral of the phase of an FIR from, does this hold true at all? Or maybe for some?


texts: How to find SNR of received DSB:AM signal with white noise after demodulation & LPF'ed?

texts:
This is a homework problem so please do not give a full solution. I intend to keep this at discussion-level.

Problem Statement

I am given the Power Spectral Density (PSD) of a DSB-Amplitude modulated signal and am told that it is corrupted with additive noise with PSD = N/2 "within the passband of the signal". That signal (with noise) is demodulated and passed through a Low-Pass Filter (LPF) whose frequency spectrum is unity between -W and W. I am tasked to find the SNR at the output of the LPF.

Attempt

SNR is defined to be the ratio of (1) the power of the message component of the demodulated received signal to (2) the power of the noise component of the same demodulated received signal.

We just need to take the ratio: (1)/(2).

Assuming the answer to question 2 (below) is yes. The power of the white noise component between -W and W is $2W*N/2 = WN$.

Questions

I'm given the (plot of the) PSD of the modulated signal. Is this (1)?
How does LPF affect the PSD of a signal? (Suppose we're given the plot of both the LPF (unity gain) passband and the PSD of the signal. Do we just point-wise multiply the PSD by the LPF function?).


texts: Computing the pseudo inverse filter from the SeDDaRA blind deconvolution process

texts:
I'm trying to apply the non-iterative method in the research paper Efficient blind deconvolution of audio-frequency signal to my ultrasound data.  I'm at the point where I believe I've calculated the magnitude coefficients
I'm not sure how to take the values in D(v), which are all reals and convert them to a pseudo inverse filter to apply the de-convolution.


texts: Special case: Band-limited in time domain and band-limited in frequency domain

texts:
I know except for some special cases, aliasing is unavoidable. Assume we time-limit a function, $f(t)$, so that it is zero outside an interval say $[0,T]$ to form $y(t)$. Then, in the frequency domain, infinite frequency components are introduced.

If the system that makes this time limitation has impulse response $h(t)$, we have its Fourier representation $H(u)$. We can then calculate $Y(u)$ using convolution to obtain:
$$Y(u)=F(u)*H(u)$$
Obviously, it has frequency components that are replicated to infinity. This phenomenon is because of sliding it across the original function to calculate Convolution.

Could anyone possibly let me know in what condition we can have band-limited input and output.


texts: periodic noise detection in image's frequency domain

texts:
Since I'm learning frequency domain and I'm really curios about this topic, I'm trying to solve a problem (that I have already solved it using morphology [detecting lines by their thickness]) in frequency space (I know it's possible). here it is the picture:

I'm trying to build a notch filter to remove lines frequencies.

I'm not sure but I think peaks on horizontal and vertical line are the noises (paper lines) and I've tried so many ways to remove them and none of them works. I'm not looking for a code and just wanna learn how it works and some tips to solve the problem. I really appreciate your kind contribution in advance.
Results after removing peaks:
Now I'm pretty sure I'm doing something wrong in detecting peaks (I've no idea what kind of peaks I'm looking for)


texts: Use GNU Radio to extract information from a VOR signal

texts:
I want to design a GNU Radio project to extract the information of a VOR signal, then find the angle to the north and the Morse code of the station.

Is there any example that shows how to do it?


texts:
I have axial slices of chest CT scans. Now, I want to use all these slices to construct a planar structure like an X-ray with a coronal view (Posterior to Anterior or Anterior to Posterior view).
I have no idea how to proceed with this problem. One way I think is a weighted average of the slices along the coronal plane with more weights to the frontal slices for AP view.
Please share your ideas on how to proceed with this problem. Thanks in advance.


texts:
I was confused about the value of group delay, I am not sure whether my opinion is correct. For example, if I want to add a group delay of $0.25$ time units to a zero-phase signal, does it mean that I just need to multiply the signal by $0.25\omega$ in frequency domain?

This is a part of my lab of DSP. In this lab, I was asked to design a LPF by least squares method. In the first question, $h_d(\omega)=1$, and it is asked to have a $0.25$ time units group delay in the next question. I tried it but I found nothing different happened, so confused.


texts: invfreqz() and frequency domain filter specification

texts:
I am trying to go from a time-domain description of a filter, via frequency domain, generating a filter based on that frequency domain description, then seeing how far from the original I end up.

L = 10-1;
h = zeros(L+1,1);
h(2) = 1;
L2 = 1024-1;
H = fft(h, L2+1);
H1 = H(1:(L2/2+2));
W = 2*pi*(0:(L2+1)/2)/(L2+1);
n = 6;
m = 2;
[b,a] = invfreqz(H1,W,n,m);

ir_in = impz(h);
ir_out = impz(b,a);
[H2,W2] = freqz(b, a);

For a simple input like this, I get a warning about singular or badly scaled matrix, and the solution (as viewed in the filter time domain) is quite different from the original:

b =
    0.0000    1.0000   -1.4219    0.2578   -0.0000    0.0000    0.0000

a =
    1.0000   -1.4219    0.2578

The frequency magnitude and phase response of the new filter looks great, but the numerical properties not so much.

I realize that setting m = 2 is kind of silly in this example, but having an automated and generic process would have been nice.


texts: Besides computational costs, is filtering, decimating, then filtering different from just filtering?

texts:
Background: I have a variety of BPSK which uses a raised cosine transmitted pulse shape:

$$ h(t) = {1 \over 2}\, (1+\cos(\pi t))\, \Pi(t/2) $$

After matched filtering in the receiver, the end result has ISI. Decimated to 1 sample per symbol, the pulse shape works out to [1/6, 1, 1/6].

So, I might filter and decimate this signal to 1 sample per symbol, then apply an equalizing filter to mitigate that ISI.

My question: modulo the potential consequences to computational complexity, is this any different from combining the matched filtering and equalization into a single filter?

I know filtering by two successive FIR filters is equivalent to one filter with the impulse response of the two convolved together. But with the decimation step between, I'm less sure.


texts:
I'm implementing an Speech Recognition System over a Cortex-M4 processor (STM32F4-Discovery) and I'm facing a problem when a compute the RIFFT.

I believe that I'm passing to the function the vector in the right form, but when I compare the output to the one obtain by MATLAB, it's quite different.

I've also tried to perform an FFT and right after that computing the IFFT, but it gives me a totally different signal as the one I gave it.

I don't know what I'm doing wrong, can anyone help me on this?


texts: Interpolated impulse response for fraction delay?

texts:
I need to find a way to create some fractional delay in a signal processing application that I am working on. Separately to this, I have been messing around with basic filter design (despite reading alot about them, including in the Scientists and engineers guide to DSP, I have still only basic understanding of filters).

During the course of this experimentation, I noticed that one of the random filters I was messing around with caused a delay of 3 samples, and that the impulse response was centred at 3 samples. I have gone back over relevant pages in the DSP guide and now understand what it says better, and that this should be the case.

So, by my current understanding, if I was to create a filter with an impulse response of:

1, 0, 0, 0, 0, 0   (starting at sample 0)

it would pass the signal unchanged.

if I had an impulse response of:

0, 0, 1, 0, 0, 0   (starting at sample 0)

it would pass the signal unchanged aside from a delay of 2 samples.

So, could I create a filter that gives fractional delay by interpolating the impulse response? As an example, if I wanted a delay of 0.75 samples, could I create a filter with this impulse response:

0.25, 0.75, 0, 0, 0, 0   (starting at sample 0)

and if I wanted a time delay of 2.4 samples, an impulse response of:

0, 0, 0.6, 0.4, 0, 0   (starting at sample 0)

I am aware that these filters would be very much sample rate dependant in terms of the real time shifts they created, but that could be easily accounted for.

Cheers in advance for any help,

Pete

EDIT:

I'm sorry if the impulse response notation is a bit crummy, I'm not really sure how the best way to write it would be


texts:
I did the following things:

Generated a 2 Hz sine wave (s1) with sampling frequency 1000 hz
Generated a 8 Hz sine wave (s2) with sampling frequency 1000 hz
combined both of them to get a new waveform (s3)
created a filter with cutoff 4 Hz and gave s3 as input.

filter used is :

[b,a] = fir1(7,4/fs);

filter_out = filter(b,a,s3);

where fs is sampling frequency .

I expected that the filter will remove the 8 Hz component from the composite waveform (s3) and the filter output will be a 2 Hz sinewave.  ie both filter output and s1 will look alike.
But this is what i got.

What is wrong with my understanding/expectation?

Thanks in advance


texts: Why are the bass components of music more audible than other components?

texts:
Consider a situation where I am listening to loud music in a room. If I move into another room and close the door (so music is now played in another room behind the closed doors), I will hear the bass component of the music much better than other components.

I would like to know why the bass voice is more audible (can be better heard through the closed door from another room) than other components in the music ?


texts:
I a recent topic, it has been told that zero-padding in the STFT can improve it, and avoid some circular convolution related things.

Here is my Python code for STFT :

fftsize = 8192;  overlap = 4;  hop = fftsize / overlap  # ie hop = 2048 here
w = scipy.hamming(fftsize)       # hamming window

# x is the input signal, I assume it is already of length 2^k, ie no zero padding required for x

# STFT
x_stft = scipy.array([scipy.fft(w*x[i:i+fftsize]) for i in range(0, len(x)-fftsize, hop)])

# ISTFT
y = scipy.zeros(len(x))
for n,i in enumerate(range(0, len(x)-fftsize, hop)):
  y[i:i+fftsize] += scipy.real(scipy.ifft(x_stft[n])) * w    # overlap-add

How can I do the zero padding with this code ? Which array should I enlarge ?

PS : in order to simplify the code, I haven't normalized the window, so for perfect reconstruction, one has to add a multiplying constant, etc. (it's not important here)


texts: STFT amplitude normalization, librosa library

texts:
So librosa.core.stft returns a complex single sided spectrogram.

My question is:

What normalization of the amplitude values should I perform afterwards?

I believe I have to multiply the amplitude outputs by 2 in order to preserve the energy that was assignated to the negative frequencies. I have not seen in the code that librosa.core.stft does this, although I might have missed it.

And then the other normalization step would be to divide the outputs by the length of the window size right?

However, after doing these two amplitude normalization steps, I still lack a factor of 2 somewhere.
This can be seen in the picture below, where in the coloured picture the spectrogram has been calculated and an amplitude value of 0.5 is observed. And in the second plot, the FT of the whole signal
has been computed and, applying the same normalization, the expected value of 1 is observed.

Picture: The frequency axis are not in Hz, the only important thing here are the amplitude values.

Here is the code:

%Import packages
import scipy.fftpack as fft
import librosa
import numpy as np

%Create sine wave of frequency 1000 Hz and 20000 samples
frequency = 1000
fs = 4 * frequency
time_vector = 1/fs * np.arange(20000)

sine_wave = np.sin(2 * np.pi * frequency * time_vector)

window_size = 1024
hop_length = 512

out  = librosa.core.spectrum.stft(sine_wave, n_fft = window_size, hop_length = hop_length)

%Take magnitude and normalize
out = 2 * np.abs(out) / window_size

fig = plt.figure()
axes = fig.add_subplot(111)
im = axes.pcolormesh(out)
fig.colorbar(im)
plt.show()

%And plot also the fourier transform of the whole signal
transform = 2 * np.abs(fft.fft(sine_wave)) / len(sine_wave)

%Plot only positive frequencies
plt.plot(transform[10000:])


texts: What is the physical meaning of the s-space?

texts:
In my audio processing lecture, the term s-space came up.

My instructor basically told us that it's a space where we perform calculations that couldn't be performed otherwise in the time space. Although he also mentioned that the s-space doesn't have a physical meaning or at least none that he is aware of, I'm curious if it has one. How exactly one can describe the s-space physically?

Also do you have a broader explanation for the s-space with some easy to grasp examples?


texts: Some questions about the intuition of the DWT

texts:
Assuming a DWT of a signal of length 8 with Haar filter taps. At the lowest level, I end up with a3 and d3 both of length 1, d2 of length 2 and d1 of length 4 which is the same number of coefficients of the original signal and which I can plot on a dyadic grid.
In contrast to the WPT in the DWT there is only basis representation for a given decomposition level ([a3,d3,d2,d1] for the 3 level decomposition in the example above). From this sequence one can choose to recover the original signal or zero out some coefficients and then reconstruct for denoising purposes and the likes.
The dyadic grid representation is analogous to the CWT and so displays signal energy over a non-uniformly discretized 2-dimensional grid of time and scale/frequency.
I often see in books a representation in which the coefficients are just concatenated and plotted on a 1-dimensional frequency axis just next to the original signal axis.

I find this representation quite confusing since it suggests to directly compare both graphs although they are displayed over time (figure a) and frequency (figure b) respectively.
Now my questions:

What is the intuition behind the stacked frequency plot of the coefficients, i.e. what is the best way to read it? Going from right to left there is a decreasing time-resolution but increasing freq resolution. In reference to the signal above, I would have the 4 d1 coefficients on the most right followed by the 2 d2 coefficients and so on. Is the correct way of interpreting this to say, that in the spectrum of [pi/2T, pi/T] I have frequency information of the original signal at just 4 different points in time (due to downsampling with a factor of 2)  for [pi/4T, pi/2T] one has freq information at only 2 points in time but being narrowed down to half the frequency interval than before and so on?

Does this mean, that the signal representation I described above has now been transformed from a time-series (the original signal) to a pseudo tf-representation and only becomes a pure time-signal again by sending it again through a synthesis filter bank? Can both representations still be considered equivalent? If one would not have the axis description, how would it be possible to recognize, if one is represented the original signal or a decomposition?

As for a WPT there are now a large number of possible signal representations and one has to keep track of the chosen one, if a later reconstruction is desired. Limiting ourselves to a level basis of level 3 in the above example, what is now the meaning of the approx and detail coefficients in every freq subband? For the DWT one can argue, that d3 are the signal details at the coarsest scale that are necessary to go from a3 to a2 and so forth. So the details capture the missed-out fine-grained signal details that were neglected along the decomposition.
Is there a similar meaning to a WPT decomposition? I have a hard time transferring this intuition to the more general signal deconstruction of a WPT

In all these transforms, usually boundary problems arise whenever using a filter that has a larger number of taps than 2. So decomposing a length 8 signal with dB2 filters already yields more than 4 values both for a1 and d1 at the first decomposition level. If I were to plot the signal at this level, would I have to cut off the additionally introduced values somehow or would the signal at this level just naturally comprise then more values than the original time signal?

Thanks a lot for any help to deepen my understanding of this complex matter


texts: sequence to perform filtering? which filter first?

texts:
Hi all Is there a generally accepted rule as to which order to perform filtering in?

for example if I have three iir filters, a low pass filter a notch and a high pass filter, which sequence would it be best to perform these in? or does it not make any difference at all?

cheers, Graham.


texts: Compressive Sensing: What Class of Signals Are Exactly Model Sparse?

texts:
In traditional compressed sensing, one would require $ m = O ( K.logN)$ measurements to recover a signal of size $N$, and sparsity $K$.

In this paper  Richard Baraniuk and his co-authors show that it is actually possible to get rid of the $log N$ factor in the number of sampled measurements, provided the $K$-sparse signals follow some additional characteristics. In traditional CS, in general, a $K$-sparse signal ( or a $K$-compressible one) lives on (or near) one of the canonical $K$-dimensional subspace of the basis of measurement. In model sparse signals, the additional requirement is that only a subset of these subspaces are allowed. These subspaces are such that they obey some well-formed characteristic. For instance, a Tree-sparse signal is one for which the top wavelet coefficients are to form a connected and rooted subtree in a wavelet tree decomposition.

I have to admit that I am still trying to understand their paper; so it may be a bit naive to ask it, but here is the question.

What sort of signals can one represent with these tree-sparse forms?  The above paper gives a few sample signals to that effect. But is there a general characterization of these signals?
For example, consider a signal comprising of, say $O(N/d)$ gaussians, with peaks randomly located among the N points ( $d < N$ is a constant in this case).  Can this be represented/approximated by one of these Tree-sparse models?

Any general pointer to literature that studied these objects?


texts: how to eliminate frame boundary discontinuity?

texts:
I'm applying at each 2048 samples the pitch synchronous algorithm to pitch scale some audio files, everything seems to occur well, but I've noticed some glitches, artifacts, clicks between every boundary frame!

I'm testing off-line in large áudio files(I'm processing frame-by-frame), what I tried to do was keep the old, current and the next frame to pitch synchronous computation (then I pitch scale 2048*3 samples), but does not seem to be working as expected.

What are the known ways to eliminate frame boundary discontinuity?