texts: Reconstructing the song using its dominant frequencies

texts:
I have found out 20 dominant frequencies of a song by splitting it into 20 sound-clips and applying Fourier transform on each one of the sound-clips. Now I am trying to reconstruct the song back from these 20 dominant frequencies.

info_ = audioinfo('../wws.wav'); % 11 secs long song

Fs = info_.SampleRate;
chunkDuration = 0.6;
numSamplesPerChunk = chunkDuration*Fs;

freq_vec = [];
for startLoc = 1:numSamplesPerChunk:info_.TotalSamples
    endLoc = min(startLoc + numSamplesPerChunk - 1, info_.TotalSamples);
    [x, Fs] = audioread('../wws.wav', [startLoc endLoc]);
    x = x(:,1);
    dt = 1/Fs;
    t = 0:dt:(length(x)*dt)-dt; % time scale
    L = length(t);
    nfft = 2^nextpow2(L);
    y = fft(x, nfft);
    y = abs(y.^2);
    y = y(1:1+nfft/2);
    [v,k] = max(y);
    f_scale = (0:nfft/2)* Fs/nfft;
    freq_vec = [freq_vec, f_scale(k)];
end

freq_vec is the vector of dominant frequencies of the song. My question is how can I use these frequencies to build back the song.


texts: does aliasing occur always if i sample a vibration in real world applications?

texts:
I was reading about the aliasing effect and nyquist. I understand that aliasing effect occurs if the sampling rate is lower than twice the maximum frequency in the signal I want to sample.

So I was wondering if aliasing is always in signal measurements. Example: If I have a vibrations sensor that has a max sample rate of 8kHz -> It can reconstruct signals till 4kHZ perfectly right?
But what about frequencies which occur also in the measurements with much higher frequencies? If there is a vibration with 16 kHz and because I sample with much lower rate, there could be aliasing in my measurement and false signals I don't know about? is this possible? That means i could have peaks in my measurements which are because of aliasing?


texts: Values (amplitudes) of the CWT and their relationship to waveform amplitudes?

texts:
I'm coding up a scalogram in C++. Right now it calculates a single row of the CWT matrix and immediately maps it to a line of pixels.

Trouble is my colouring scales with the values of the transform and naturally brightest colour = highest value. Yet as it is I don't know how to find the maximum value so that I can generate my spectrum. I know CWT is funny with amplitudes so I'm asking for guidance.

Thanks in advance,
a meme


texts: Modelling an existing DSP filter by analysing the differences between its inputs and outputs

texts:
A DSP filter produces always the same output.

$$\boxed{\textrm{wav file}}{\longrightarrow}\boxed{\textrm{DSP}}{\longrightarrow}\boxed{\textrm{wav file(modified)}}$$

The DSP is a black box as I don't have access to its code, but I have access to original and modified wav files.

Questions:

Can the DSP filter be mathematically modelled (say in MATLAB) by comparing the original and modified wav files?
If so, would the resulting modelled DSP, when applied to any wav file, generate a wav equals to the one that will be produced by the original DSP?


texts:
I'm using a sin wave $s_{sc}(t) = d(t)\cos(\omega_{sc}t + \varphi_{sc})$ to modulate a carrier using phase modulation. I'm working on Matlab so I use the pmmod function to do so. It is not too complicated to understand what this function is doing, probably computing something like $\cos(\omega_c t + h s_{sc}(t) + \varphi_{ini})$ (with h the modulation index).

In order to demodulate it, I'm using the pmdemod function of Matlab. But for this one, I have absolutely no idea of what is happening in background. I tried to analyze the results obtained after a pmmod-pmdemod sequence, and this is what I observed on Matlab:

Regarding the convergence time, it seems that Matlab use a PLL in background. but what kind of PLL is it ? I'm very interested about the digital methods used to demodulate PM modulation and more particularly the maths behind it. But I cannot find any article/document/book about it (Or I just don't know what to search on google, because when I type "phase modulation" I just find article about CPM modulation or digital PSK/FSK modulation methods...).

Any help would be helpful :)


texts: Understanding the mathematical proof for the alias frequencies in a sampled sine wave

texts:
I'm struggling to get my head round the mathematical proof for the alias frequencies in a sampled sine wave.

I understand that sampling a sine wave of frequency $f_0$ every $t_s$ seconds gives you:

$$x[n]=\sin(2\pi f_0nt_s)$$

I also understand that, because the sine wave is periodic every $2\pi$, you can add any multiple of $2\pi$ to the angle and get the same values for the sine, i.e.,

$$\sin(2\pi f_0nt_s)=\sin(2\pi f_0nt_s+2\pi m) \quad\text{(where $m$ is any integer).}$$

The proof I'm looking at then factors out $2\pi$ and $nt_s$ to get:

$$\sin\left(2\pi(f_0+\frac{m}{nt_s})nt_s\right)$$

...but then it says to let $m$ be an integer multiple of $n$ so we can replace the $\frac{m}{n}$ ratio with an integer $k$.

I don't understand how $m$ can go from being "any integer" to "an integer multiple of $n$". If $m$ is any integer and $n$ is an integer then how can the ratio between them be an integer?

I know I'm missing something obvious here and I'm searching for that light-bulb moment but it's not happening. Because this is so fundamental to DSP I don't just want to accept the formula and move on without thoroughly understanding it.


texts:
Let $s(t)$ be a signal that can be approximated by a uniform spline
  function of order $K$ (say $K=2$): $$s(t)\approx\sum_{n\in \mathbb{Z}}c_n\beta_+^{(K)}(t-n) $$

  Suppose that we know the values of $s(t)$ at integer locations
  $s(n)=x[n]$. Develop and implement a filter-based algorithm to recover
  sequence $c_n$ from $x_n$.

  With the coefficient sequence $c_n$ we can easily compute the sampled
  version of the derivative or fractional translation of $s(t)$, i.e.
  $s(n)$ and $s(n+\tau)$ where $\tau$ is a fractional number. Test your
  results with function $s(t)=\sin(2\pi \cdot 0.001\cdot t)$for $t\in [0,100]$.

Check the attached problem please. I am a beginner in spline fitting and have a few questions:

1) How to find the coefficients $c[n]$. Is it by DTFT?

2) I understand how to find the derivative but not sure about how to deal with fractional translation?

3) What's the right function in Matlab to obtain the (discrete) samples of a uniform spline function

Thanks!


texts: Finding the discrete time fourier coefficients to this problem

texts:
I'm trying to find the fourier series to this discrete time signal.
$$x_1[n] =\begin{cases}
+\frac72&\text{if }0\le n \le 4\\
-\frac72&\text{if }5\le n \le 9
\end{cases}$$
My approach:
We see that the discrete period, $N$, is 10, so $\Omega_0 = \frac{2\pi}{10} = \frac{\pi}{5}$
Now using the formula to get the DT fourier coefficients, $X[k]$
\begin{align}
X[k] &= \frac{1}{N}\sum_{n = 0}^{n = 9}x_1[n]e^{-j\frac{\pi}{5}kn} \tag{1}\\
&= \frac{1}{10}\left( \sum_{n = 0}^{n = 4}\frac{7}{2}e^{-j\frac{\pi}{5}kn} + \sum_{n = 5}^{n = 9}\frac{-7}{2}e^{-j\frac{\pi}{5}kn} \right)\tag{2}\\
&= \frac{7j}{10}\sum_{n = 0}^{n = 4}\sin\left(\frac{\pi k n}{5}\right)\tag{3}
\end{align}
However according to the solutions for $k = 2$, we have $X[k = 2] = 0$, but when I substitute $kk = 2$ into the final form of the above expression, I don't get $0$. Thus, my answer must be wrong, but I'm not sure where I went wrong... Can someone please explain where my math is flawed and if there's an easier way to find the coefficients?


texts:
I have attempted question 1 a) below in MATLAB 6.5, but I need my answer to be checked.

How do I do question 1 b) below ?
What does it mean to have a zero mean?

  Generate the following signals with a computer program:

  a)   Random noise in the amplitude interval (0,0  1,0)
  b)   Random noise in the previous question but with a zero mean

function [] = Question1a()
datapoint = 1000;
t=0:1/datapoint:2;
x=sin(2*pi*5*t);
y=sin(2*pi*25*t);
z= x + y;
NoiseAmplitude = 1;
Noise = x + NoiseAmplitude*randn(size(z));
plot(t,Noise);
title('Noise in the Amplitude interval (0,0,1,0)')
ylabel ('Amplitude')
xlabel ('Time(s)')


texts: What image-processing techniques are ideal for this particular shift/scale invariant template matching?

texts:
A problem I had originally discussed here has evolved, and might have gotten a little simpler as I studied it in some more and attained new information.

Bottom line, I would like to be able to detect this pattern shown here, using computer-vision/image-processing techniques. As you can see, the ideal pattern is made up of four 'pings'. The object-recognition should be:

Shift Invariant

Horizontally, the image will be cyclical. (ie, Push to the right, comes out on the left, and vice-versa).
(Fortunately) Vertically, it is not cyclical. (ie,Push to either top or bottom and it will stop).

Scale Invariant (The pings can vary in 'thickness' as you can see.)

I could go on and on about it but I have attached images covering what I mean, please see below:

Of course, they can also be at a different 'scale', as can be seen from this family:

And finally, here are some 'realistic' scenarios of what I might actually receive, where there would be noise, the rows might 'fade' as you go towards the bottom, and of course, the  image will have plenty of spurious lines, artifacts, etc.

And of course, as a grand finale, there is the distinct possibility of this 'extreme' scenario:

So once again, I would like to ask for some guidance on what computer-machine-vision techniques I should utilize here in order to best detect the occurance of my pattern, where I need to be shift and scale invariant as you can see, and also be able to get decent results for the realistic scenarios. (Good news is that I DONT need it to be rotationaly invariant). Only idea I have been able to come up with thus far is a 2-D correlation.

I should add, that in reality, I do NOT have colored images - I will just get a martix of numbers, so I suppose we are talking about 'greyscale'.

Thanks in advance!

P.S. For what its worth I will probably be using open C V.

EDIT # 1:

Based on the comments, I have added the details you requested here:

For characteristics defining the data, we can assume the following:

Horizontal length of each ping can vary, but I know the upper and lower bounds of it. YES for anything within this bound, NO for anything outside. (Example, I know the length of the pings can be anywhere between 1 and 3 seconds, for example).
All pings need to be 'visible' for YES, however, last row can be missing, and still want to say 'YES'. Otherwise NO.
Vertical length, (the 'thickness') of each ping can vary, but again, also know upper and lower bound. (Similar to what you see in those images). YES for anything within that bound. NO for anything outside.
Height between each ping should always be the same for YES. If they are not, then NO. (Example, you can see how all the pings are of the same height relative to each other, ~110 on the vertical axis). So so 110 +/- 5 can be a YES, anything else has to be NO.

I think thats about it - but let me know what else I can add... (Also, everything shown here should register as a YES, btw).


texts: How is it possible to Filter out a person's voice out of 100 of other voices?

texts:
So, I've just learned that the human voice is not a single sine wave its a bunch of unlimited Sine waves each having a different frequencies,

According to Wikipedia,

  The voice consists of sound made by a human being using the vocal
  folds for talking, singing, laughing, crying, screaming, etc. Its
  frequency ranges from about 60 to 7000 Hz.

So if the human voice is a composite signals, it contains various frequencies ranges between 60Hz to 7KHz.

Suppose if there is a group of people singing the same song all together, each person has its own ranges of voice frequencies,

For example,

If a person A has the following frequencies,
100Hz , 250HZ 6KHz, 10Hz, 87Hz, 52Hz, 2KHz.......

and a Person B has the following,
217Hz1, 11Hz, 12Hz, 2323Hz, 839Hz, 4KHz, 100Hz, 10Hz.....

there must be so many frequencies which are similar in both the person A & B, like in above example the frequencies 100Hz and 10Hz are common between two persons.

I was watching a TV Show name "Fringe" where they filter out the particular Man's voice from an audio file while there were other people voice present there too.

So how does exactly they filter someone's voice out of the voice of 100s of people if there are so many frequencies common among all of them does it have to do something with the amplitudes of person's frequencies ?


texts: Inverse of the Gibbs Phenomenon to generate a sine wave

texts: I am a quite new to the area of digital signal processing. I am currently doing a project on sine wave generation. I know that according to Gibbs phenomenon, many sine waves add up to form a square wave. So, square wave would contain many frequencies. If we filter that square wave stage by stage using low pass and band pass filters, will I be able to get a sine wave at expected fundamental frequency? Please help.


texts: How to Detect Gradients and Edges in Images?

texts:
I want to be able to find points in images that are the centre of a radial gradient like the one shown in the left picture below. Any ideas on how I could use a Hough transform or some other computer vision method?

Thanks

example search image:


texts: What is the frequency of $\cos x -\sqrt 2\cos\sqrt 2 x$?

texts:
Question

Considering that $1\over{2\pi}$ is the frequency of $\cos x$ and also of $\cos x - 2\cos 2x$, what is the frequency of $\cos x - \sqrt 2\cos\sqrt 2 x$?

Thoughts

Perhaps "frequency" isn't the interesting quantity. I think the period is infinite. That would make the frequency zero. Maybe what I want to know is bandwidth?

I think the additive combination (superposition) of the two waves produces a component at the higher frequency $1/(2\pi(\sqrt 2-1))$. I would hope that sampling at a frequency of $1/(\sqrt 2\pi -\pi)$ would be sufficient to reveal as many zero-crossings as we may wish. Is that correct? Is it the best answer?

I understood Fourier series and partial differential equations when I took that course and I still get the general idea, but it has been a long time.

Motivation

Mistakenly, I devised and posed this question for my colleagues as a weekly math challenge before I realized that I don't know how to get the answer myself.

I'll provide additionally just some background about my motive for asking this particular question. Maybe it can help you give an answer that is relevant to me. However, this background information in itself is not part of my question. So, I came up with the question about frequency hoping to move us nearer to another question that Mark Kac used in order to elucidate stochastic independence. That is, given any particular value $-2 <\alpha < 2$, in what fraction of the interval $0\leq x < 100$ (just to choose some arbitrary bounds) is $f(x) <\alpha$ where $f(x) = \sin x - \sin\sqrt 2 x$? The interesting destination that I want to reach is that the answer to this further question is (somehow related to) the normal distribution function, yet there is nothing "random" involved.

Of course $f'(x) =\cos x -\sqrt 2\cos\sqrt 2 x$. If I know the "frequency" (or something) of $f'(x)$ then I hope I can find, at least numerically, its zero-crossings. Those include all of the local extrema of $f(x)$. Then, knowing that the sum-to-product formula gives $f(x) = 2\cos(x/2+x/\sqrt 2)\sin(x/2-x/\sqrt 2)$, I think that my favorite computer algebra system and I can find, at least numerically, where $f(x) <\alpha$.


texts: Can you recommend an efficient encoding algorithm for sound with long periods of silence?

texts:
As a side personal project, I'm coding an Android application that I plan to use to record the sounds of our first baby. I'm kind of new to the world of Android development and sound processing, so please bear with me if this question is naive :)

What I'm trying to do is to constantly save the sounds in our baby room, splitting the recordings in 2-hour intervals. I plan to keep this data for future analysis (e.g. plotting crying times, which I'll correlate with other events, and maybe cry identification).

So far the application is working fine. However, I wonder if there is an efficient codec I can use considering that 90% of the time there will be silence. I've tried modifying the PCM raw data before enters the codec setting to zero any values below a specific threshold, but this did not reduce the file size substantially. This was a wild guess, just in case this would "help" the codec (I've tried AAC/HE and MP3, not yet AMR).

Can you please point me in the right direction?...

I don't mind using a solution that requires post processing of the files. For example, I've looked into http://sox.sourceforge.net, but that strips silence periods reducing the recording length (I want to keep the timestamps and length untouched).

Note: using SoX to split the file in parts using a naming scheme that allows me to reconstitute original times it's an option... but I could not come with that script, yet).

Bonus question: will AMR-NB be enough to encode a baby cry, if I want to run some processing later? (for example this or this) or... should I stick to AMR-WB?

Thank you,

Update July 22nd:

Thanks everyone for the responses so far. The idea of using just WAV injecting zeros, and later using bzip2, etc. to compress the repetitive data is really good (thanks @MBaz).

@Fat32: I've erased the files, but I recall sizes around ~20Mb/hr, with no substantial change after replacing low PCM values (approximately 50% of the stream) with zeros.

Another update: AMR-WB is giving me better results than AAC/HE, ~6 Mb/hr. So I'll probably stick to it, however I'm not sure how the low bitrate + sampling rate will affect things if later I decide to run some infant identification analysis later (i.e. if the loss of information will remove data points crucial for the neural net).


texts: How Can Convolution and Deconvolution be Defined for 3D Images?

texts:
I am trying to understand how convolution and deconvolution can be represented for 3D images/ stacks of data. I would prefer it, if you built the these concepts from 1D vectors to 3D matrices in terms of their equations. I am currently writing a report and am trying to represent 3D deconvolution as an equation.

Please do not hesitate to ask me any questions to clarify this question.


texts: Received waveform frequency is way too low than it was sent

texts:
I am using Quartus 12.1 sp1, vhdl and Altera Nios II programmed in C code for DE0-Nano Development. Basically, what I have is, data is sent from fpga vhdl block to NIos II system, Nios II system sends the data to serial port, Matlab access the serial port to real time plot the graph.

In my system, the rate the graph is plotted is much slower (11Hz) than the rate the processor speed (50MHZ). So I do some buffering and send the data for plotting afterwards.

I am buffering data received from fpga vhdl say 1000 sample points and save it in sdram... after saving these 1000 sample points, the C code will send data to uart so that Matlab can access these serial port data and plot the graph real time.

Please refer to this image:

My problem is, a complete cycle of sinewave obtained from fpga is with a frequency of 50Hz, so it should be having a period of 0.02 seconds instead of 0.16 minutes as I got now, I don't know what could be the problem? Forgive me, I am newbie to signal processing.. Thank you.

int max_sample = 1000; //number of sample points wanted

if (a_out_data==1) //check if it is the beginning of a sinusoidal waveform?
{

for(i=0;i<max_sample;i++)//start the process of buffering data
        {
vref_data = IORD_ALTERA_AVALON_PIO_DATA(vref_base); //read a single value of reference voltage from fpga
*(pSDRAM+i) =vref_data; //save this single value in a location of SDRAM specified in SDRAM_U1_BASE
        }
}

for(i2=0;i2<i;i2++) //start the process of sending data to uart for all data saved previously
    {
        data16=*(pSDRAM+i2); //access each value from SDRAM

        //send data to uart before increment i++
    }

EDIT: I posted the Matlab code here:

max_packets = 20000; %maximum number of packets to keep in memory (all are saved to log file)
    update_plot = 1; %how often we should update the plot
    plot_length = 12; %maximum length of plot
    figure('Position', [10, 50, 900, 900]); % new figure
    ax1 = subplot(3,1,1); % top subplot
    ax2 = subplot(3,1,2); % bottom subplot
    ax3 = subplot(3,1,3);
    drawnow;
    s = serial(port);
    cleanupObj = onCleanup(@()cleanupFunc(s)); %clean-up function
    set(s,'BaudRate',115200,'timeout',2); %timeout is in seconds
    fopen(s);
    vref= NaN(12,max_packets);
    j = 0;
    tic; %start timer
    while 1; %do forever...
          if j < max_packets
          j = j+1;
      else
          fprintf('\n Buffer full, reset to 0\n');
          j = 1; %reset counter to stop memory overflow
      end
        %get data line by line
        line = fscanf(s);
        [c,vo] = disect_line(line);
        if c == 1 %start recording from set 1
            vref(1,j) = vo;
            for i = 2:12 %get set 2 to 12
                line = fscanf(s);
                [c,vo] = disect_line(line);
                if c == i
                    vref(i,j) = vo;
                else
                    fprintf('\nError\n');
                end
            end
           dlmwrite(logfile,[toc  vref(:,j).' ],'delimiter',',','-append'); %save data from all 12 sets to csv file
       update plot1 every update_plot packets, and after we have at least 2 data points (otherwise we get errors)
          if mod(j, update_plot)==0 && j>1
              if j < plot_length
                  start = 1;
              else
                  start = j-plot_length+1;
              end
           plot(ax1, start:j, vref(1,start:j).');
                for k = 1:1
                    text(j, vref(k, j), num2str(k), 'Parent', ax1);
                end
                xlim(ax1, [start j]);
                ylabel(ax1, 'Vref');
                title(ax1, ['Time elapsed = ' num2str(toc) ' s']);
    end
    end
    function [c, vo]  = disect_line(line)
    c = str2double(line(1:2)); %set number
    vo = str2double(line(4:7));
    end
    function cleanupFunc(s)
    fprintf('\nClosed serial port.\n')
    fclose(s);
    delete(s);
    clear s;
    end


texts:
I need to compute the sine and cosine of an argument along with n "harmonics"
\begin{matrix} \sin(x) & \cos(x) \\
  \sin(2x) & \cos(2x) \\
\cdots \\
\sin(nx) & \cos(nx)
\end{matrix}
This takes a lot of computing time and since this is a real-time system, I need to optimize this calculation.
Here's my solution so far :
I defined a rotation matrix $$
R = \begin{pmatrix}
    \cos(x) \ -\sin(x) \\
    \sin(x) \ \cos(x) \\
    \end{pmatrix}
$$
I compute $[\cos(nx); \sin(nx)] = R  [\cos((n-1)x); \sin((n-1)x)]$; I simply iterate to compute the missing $n-1$ harmonics.
This method is 3 times fast than computing all harmonics separately.
I'm curious, are there better solutions? I know a look-up table could be even faster but I would prefer not to go this way.


texts: What is sampling frequency in this case?

texts:
I am using Matlab fdatool to implement highpass IIR butterworth filter
But 'Fs' value cell is gray and i am not able to change sampling frequency, what is the actual value in this case of attached photo?


texts: How to set STFT-parameters to visualize woodpecker / Sawtooth Signal?

texts:
I always wanted to work with FFTs and Spectrograms to characterise sounds and their frequencies using python.

Lately I recorded a woodpecker at his work.

Woodpecker Wav File

To know where to look at I compared the sound with a generated sawtooth from Onlinegenerator.com. I think 17 Hz fits well.

17 Hz Sawtooth Wav File

Checking the code with the generated file worked good.

But sadly I cannot find the frequency in the DFT/FFT or spectrogram generated with the code below. No matter which settings I tried, I cannot visualise the frequencies below 20 Hz properly. No major signals cuts on or off. Only a signal at 23 Hz, which seems to be continuous.

from scipy import signal
import scipy.io.wavfile as wav
import matplotlib.pyplot as plt
import numpy as np

# Read WAV File
the_file = 'IMG_0710_short.wav'
samplerate, samples = wav.read(the_file)

#samples_left = samples[:,0]
samples = samples[:,1] # reduce to right channel

timevec = range(len(samples)) # Time vector for plot
timevec = [x / samplerate for x in timevec]
t_max =timevec[-1]

# Lineplot of Signal
dpi = 200
plt.rcParams['figure.dpi']= dpi

plt.plot(timevec,samples)
plt.title(the_file)
plt.ylabel('Amplitude [-]')
plt.xlabel('Time [sec]')

# STFT Settings
nperseg  = 0.5 * samplerate # Window Size 1
noverlap = nperseg*0.95 # Overlap
nfft     = 1 * samplerate # STFT 2
window   = 'hann' # Window Type

timespan = [0, 2] # Calculation Window
fromm = int(len(samples)/t_max*timespan[0])
too   = int(len(samples)/t_max*timespan[1])

f, t, Zxx = signal.stft(samples[fromm:too], samplerate, nperseg=nperseg, window=window, noverlap=noverlap, nfft=nfft)
t = t + timespan[0]
cmap=plt.cm.nipy_spectral
vmin = 10
vmax = 18
fig = plt.figure(figsize=(7, 5))
pcm = plt.pcolormesh(t, f, np.log(np.abs(Zxx)), cmap=cmap, vmin=vmin, vmax=vmax)

plt.ylim(0,50)
fig.colorbar(pcm)
plt.title(the_file)
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')

# FFT
n = len(samples) # length of the signal
k = np.arange(n)
T = n/samplerate
frq = k/T # two sides frequency range
frq = frq[range(int((n/2)))] # one side frequency range

Y = np.fft.fft(samples, norm='ortho')#/n # fft computing and normalization
Y = Y[range(int(n/2))]

fig, ax = plt.subplots(2, 1)
ax[0].plot(timevec,samples) # plotting the signal
ax[0].set_xlabel('Time')
ax[0].set_ylabel('Amplitude')
ax[1].plot(frq,abs(Y),'r') # plotting the spectrum
ax[1].set_xlabel('Frequency (Hz)')
ax[1].set_ylabel('|Y(freq)|')
ax[1].set_xlim(0,50)

How do I need to set the stft parameters in order to visualise the woodpecker spectrum correctly?
Or is there a fundamental problem with stft here?


texts: Why can't I add variable covariance noise to this 2D kalman filter?

texts:
The goal is to include variable measurement noise to my kalman filter.
To begin with I kept the noise covariance matrix as follows:
$$
R=\begin{bmatrix}1&0&0&0\\
0&1&0&0\\
0&0&1&0\\
0&0&0&1\end{bmatrix}$$

This gave good results, the Kalman filter did a good job following the input data.

Next I tried to calculate the measurement noise covariance myself:
https://pastebin.com/7ETMC0Fy

Which gave such results:

Could someone tell me why when I try to calculate my measurement noise covariance matrix myself the result is unsatisfactory? Why doesn't the output of my Kalman filter follow the input data at all?


texts: Analytical expression for salient Blob detection in scale space

texts:
I am working on exact mathematical expression which can be obtained in image processing tasks.(I do not know much about image processing but my work is more mathmatical) In Gaussian scale space, we obtain the salient blob like structures in image by finding the extremum points in scale-normalized Laplacian of the image in the famous SIFT method. Can anyone give me reference to

1) how can we represent this condition mathematically, for e.g. this will be the set of points given as:

$$
S=\{s= (x,y,t): \text {s is the point of extremum of }\sigma^2 \nabla^2 L(x,y,t)  \}
$$
should I write it as:
$$
S=\{s= (x,y,t): \frac{d[ \sigma^2 \nabla^2 L(s)]}{ds}=0 \}
$$

2) How have we reached to the fact that extremum in scale-normalized laplacian will give us the salient blob like structures? I have read some papers e.g. Blostein and Ahuja and T lindeberg "scale space theory in computer vision", where it is mentioned that you can use laplaican operator to obtain such features, but even in that paper I cannot find any mathmatical expression leading to such conclusion.

If you see paper like "scale and differential structure of images" by florack et. al., they have mathematically proven how we obtain Gaussian kernel in scale space. Why is there no such mathematical treatment of image features like salient blob detection in SIFT and edge detection in non-linear scale space by KAZE(where you use similar extrema for determinant of Hessian.). If there is, can someone give reference of the paper or book?

I will be very thankful if someone can answer this.


texts: why LTE doesn't use DC subcarrier in the downlink?

texts:
I know LTE doesn't use the DC subcarrier for OFDM/OFDMA modulation for the downlink, but no author gives a detailed explanation to this stuff.

Can you tell me the reason behind it?

Thank you!


texts:
I have a project in which I have to recognize 2 lines in different "position", the lines are orthogonal but can be projected on different surfaces. I'm using opencv.

I'm not asking for  code, but only for hints about how can I solve this? I tried houghlines function but it works only for straight surfaces.
thanks in advance !


texts: how to find the distance of the root from the origin?

texts:
Good Day, I have 9 polynomial roots. i have plotted the roots. my problems are:

1)how to find the distance of the roots from origin? I want to find the distance of all the roots from the origin.

2) I have 3 models and I want to differentiate between them using the polynomial roots. I want to observes how the roots are moving from each models. Is it possible for me to arrange the roots in ascending order? the one that is closer to the origin and farther from the origin.

Any idea? Looking for your reply. Thanks!


texts:
Does anyone know of a filter to attenuate non-speech? I am writing speech recognition software and would like to filter out everything but human speech. This would include background noise, noise produced by a crappy microphone, or even background music. I have already implemented a first order filter that compensates for the 6 dB roll-off of the power spectrum, but I'm still hearing noise (though the speech sounds a lot clearer). I have thought to use a low-pass filter, but I am iffy about doing that for two reasons:

I don't know whether or not a low-pass pre-filter will interfere with the rest of the speech processing. Even though the human ear can only detect sounds lower than 20 kHz, I don't want to risk eliminating any higher order harmonics that might be necessary to process the speech (though I don't know if this is the case or not. But I don't want to take any chances).
I understand that the excitation of certain consonants (such as f, h, and s) are almost entirely white noise. I don't want to implement a noise filter that will eliminate good noise, so to speak.

Ideally, I would like to be left with only the speech of the person speaking into the microphone. If you have any ideas, or there is something that I'm missing, please let me know. Much appreciated!


texts: Discrete Fourier Transform and Opposite Convolution Theorem

texts:
I am reading the Wiki for DFT. There is a part for circular convolution theorem which sounds a bit odd saying:

$$ \mathcal{F} \left \{ \mathbf{x\cdot y} \right \}_k \ \stackrel{\mathrm{def}}{=} \sum_{n=0}^{N-1} x_n \cdot y_n \cdot e^{-\frac{2\pi i}{N} k n}
=\frac{1}{N} (\mathbf{X * Y_N})_k, \,$$   which is the circular convolution of $\mathbf{X}$ and $\mathbf{Y}$.

I feel it's incorrect because writing $y_n$ in Fourier basis, and multiplying it by $e^{-\frac{2\pi i}{N} k n}$ will result in a shift in index by $k$, and that means the result should be the cross-corrolation of two signals.

And secondly if either this or my argument is correct why this does not hold for continuous signals.

P.S. The question has been asked here as well since I was unaware of the most proper domain for my question.


texts: 3D data centering before fast fourier transform

texts:
I am using fftw to conduct a 3 dimensional fft and then apply a guassian filter. In order to shift the frequency, I need to first rotate the data matrix before fft. In 2D, I know how to do it, e.g., by multiplying (-1)^(x+y). But how to do it in 3D? Is it (-1)^(x+y+z)?
Thanks


texts: is there a single process/term to describe dissipation and amplification of a signal?

texts: I am writing a paper, and the topic is about a signal that could be either amplified or dissipated while traveling through a medium. Since I don't wanna write everywhere "amplified or dissipated ", e.g. "the plot shows the signal amplification or dissipation", is there a single word I could use for both? Maybe amplitude modulation? I don't think it is correct cause there is not a modulating low frequency signal.


texts:
I don't really understand FFTs, and filtering neither. I hope I can explain my confusion in a way that is not too confusing for you to understand.

consider the following example: I have an EEG signal, and am interested in the frequency bands from 1 to 40 Hz.

if I apply a LP filter at 40, and a HP at 1, my signal of interest would be distorted around the beginning and end because of the filter's transition band. right? why does it have a transition band at all, and how does it work? is a filter a sum of sine waves in the specified frequency range? (sorry, seems like a basic question, still I couldn't find an answer to it)
now to get the FFT of different frequency bands in the 1-40Hz range, I could either first bandpass filter the specific ranges and apply FFT to each (e.g. to 1-4,4-8,8-12,12-20). however, this would lead to a distorted result in each band, because of the transition bands, right?

if on the other hand, I do FFT over the whole 1-40 range, and then just 'extract' the results from 8-12 etc, the results would be more accurate? and how would the extraction mathematically work, is it just the mean of all values at the frequencies between 8 and 12?

I hope you can help me understand these things


texts: Possible to Derive Receive Signal Power from RSSI?

texts:
My Set-Up

I am sending packets at a specific frequency and with a specific gain from one N-210 USRP to another of the same model.

NOTE: I have the transmitting USRP connected to a series of attenuators (so that I don't overwhelm the spectrum analyzer) and then to an RF splitter - one output is connected to a spectrum analyzer and the other output is connected to the receiving USRP. This is how I am reading the receive signal's magnitude.

My Problem

I am attempting to calibrate the receive signal power measurements (in dBm) that I can visually record from the spectrum analyzer with the RSSI measurements (in dB) that are output into a file by a script (I read them from the text output of the ofdmflexrframe_rx function). My plan is to come up with a linear correction factor to relate RSSI measurements to the receive signal power measurements from the spectrum analyzer.

I read on Wikipedia:

  RSSI is the relative received signal strength in a wireless environment, in arbitrary units.

I also read this from a response to a USRP user's question about RSSI measurements:

  [The] Received Signal Strength [Indicator is] always relative to some
  signal model, incorporating considered bandwidth, assumptions on the
  modulation scheme, duration of transmission, generally: It's a
  estimation of received signal strength based on some property of
  the received signal.

Now, I think that comparing the receive signal power levels on the spectrum analyzer with the RSSI measurements that function provides is not the right way to come up with a correction factor.

What is the difference between RSSI measurements and the dBm values that I record from the spectrum analyzer?
Is it possible to relate these RSSI measurements to the receive signal power levels that I record from the spectrum analyzer?
If I cannot relate these two measurements by a correction factor, is there a way to obtain the receive signal power measurements (in dBm) programmatically (GNUradio perhaps)? (NOTE: I don't believe there is any way to obtain the receive signal power measurements through the library that contains the function I linked above)


texts:
I am trying to understand how SURF extracts features using a Hessian matrix.
What I am a bit confused about, is why a second derivative Gaussian filter is used, rather than a standard Gaussian filter? What effect does this have on the system?


texts:
SKIP to Questions below if you want :

In super resolution I can create a high resolution signal from a low resolution signal.
But it is also possible that the low resolution signal contains aliasing. The problem with aliasing is that it cannot be distinguished between the real low resolution frequencies. But in super resolution you can make multiple low resolution measurements with aliasing and calculate the shift between the signals (afaik) that can help to reconstruct a higher resolution signal. Check some papers like here :
Multi-Frame super resolution,
SR from aliased images

What I am thinking is if I can make multiple shots from one scene with aliasing and construct a higher resolution image out of it, how could I apply this to vibration from a machine.
In the first paper google shows how they use the vibrations from a hand for taking multiple photos from one scene. The vibrations from a human hand help to capture different frequencies with aliasing. Later on, all merged can be used for super-resolution.

Questions:

I have a fixed sensor on a machine which basically means I have a
camera fixed on a stand(compare photo camera wich takes multiple pictures at same time from the same scene with aliasing to construct a higher resolution). How would this affect the super-resolution
in comparison to using different angles in parallel to have multiple measurements with different offsets?
Could I apply this on
vibration with making a long measurement, dividing the measurement
in different blocks and trying to merge them as they where different
captures, instead of using different angles?


texts: Where can I get passport size face photo database for my image classification experiment?

texts: I am just trying to classify images using kNN algorithm with the aid of MATLAB. So now I am in need of some collection of images of peoples' faces. Where can I get some open database of faces (both females and males).


texts: Comparing PSD two signals of different length

texts:
I have two time signals A and B which are of different length, length(A) < length(B). Now I calculated the cumulative PSD of both signals. What is the best way to compare this result?

If I am correct you can not compare them directly because they are of different length. What is the best to compare them? I tried normalize them by dividing the PSD result by the length of the signal, but I am not sure this is actually a valid approach.

Note, I am not interested in the common frequencies I am interested in the energy content.


texts: Understanding the over sampling factor effects before a Matched Filter?

texts: I am trying to design a matched filter at the receiver side to a half sine pulse shaping function sent by the transmitter of an OQPSK transceiver. I am not sure I understand the effect of varying the over sampling factor of the ADC right before the matched filter (the minimum is twice the symbol rate). I think it may have something to do with the quality of the timing synchronization. I'm using an early-late gate synchronization algorithm, which needs at least 3 samples per symbol. I've also read in literature that the over sampling is typically a factor of 8 for QPSK systems. Does it have other effects, maybe the BER of the system?


texts:
Background: The question here is related to images in particular and not signal/waveforms.

I have been reading a lot of answers about the difference between convolution and correlation but I am stuck on my doubt as to how are they related to each other while performing operations. Most of the answers are individual answers for each of these terms.

For example:

A statement I am trying to understand says

  Here correlation operations are performed by replacing the exhausted
  convolutions with element-wise multiplications using Discrete Fourier
  Transform (DFT).

I am not able to relate the individual definitions here. Can anyone please explain this?


texts: How can I design Nyquist interpolation filters with the Parks-McClellan algorithm?

texts:
We can easily design interpolation filters that obey certain frequency-domain constraints using the Parks-McClellan algorithm.  However, it's not immediately clear how to enforce time-domain constraints; in particular, I'm interested in generating Nyquist filters.  So if I'm oversampling by a factor of N, I want the filter to have zero-crossings at kN, for non-zero integer k (this ensures that the input samples to my interpolator will appear in the output sequence).

I've seen Harris1 talk about a technique for designing half-band filters, i.e. the special case where N=2.  Is there a general solution for this?  (I know that we can easily design filters with the window method, but that doesn't give us the same control.)

[1] Multirate Signal Processing for Communication Systems, pp. 208-209


texts: What Is the Difference Between MRF and Total Variation in Noise Removal?

texts:
I have a general question in image processing. I have a noisy image. I would like to classify the noisy image into some regions. Two famous approaches I can use are:

MRF/Gibbs MRF: Model the spatial dependence between neighborhood
pixels.
Total Variation: Key idea may be based on smallest variation
of the image.

My question is: Could you tell me what is different between two approaches for noise removal? What is the key idea of TV? Which one is better? Thanks.


texts: What is the correct representation for discrete time sequence?

texts:
I am bit confused after reading Frequency-Domain Sampling and Reconstruction of Discrete-Time Signals section from Proakis where, $$x(n)=\begin{cases}
               x_p(n),\;\;\;\;0\leq n \leq N-1\\
                     \\
               0,\;\;\;\;\;\;\;\;\;\;\;     \mathrm{elsewhere}\\
            \end{cases}$$
It is clear $x(n)$ is discrete time sequence but why $n$ is defined as $0 \leq n \leq N-1$ instead of $n=0,1,\dots,N-1$


texts:
Given below is a function written in MATLAB for computing the modified periodogram, taken from Hayes, Statistical Digital Signal Processing and Modeling

function Px = mper(x, win, n1, n2)

x = x(:);
if nargin == 2
    n1 = 1; n2 = length(x);
end

N = n2 - n1 + 1;

switch win
    case 1
        w = rectwin(N);
    case 2
        w = hamming(N);
    case 3
        w = hanning(N);
    case 4
        w = bartlett(N);
    case 5
        w = blackman(N);
end;

xw = x(n1:n2).*w/norm(w);
Px = N*mypdg(xw);

Can anyone explain what is going on in second last line of the code
xw = x(n1:n2).*w/norm(w)

Is dividing by norm necessary to compute xw?


texts:
I am a mathematician trying to understand some signal processing terms. My first question regards the terms 'even' and 'odd' filters. What does this mean? Does it reference the fact that upon a Fourier transform, there is a cosine and a sine function, which in turn are even and odd? If not, a link to a reliable source would be appreciated.

Second question: What is a quadrature? I have heard the term used for the phase shift between two Fourier transforms, the method for dividing a frame into four parts, and taking a gradient-like subtraction from vertical, horizontal, and both diagonal directions. Is it all of these at once?
Also, a link to a reliable source would be great.


texts: system function $H(\omega)$ relationship to odd and even components of h[n]

texts:
What qualities of $h[n]$ are necessary for:

$$
H(e^{j\omega}) = DTFT\{h_{even}[n]\} + j\ DTFT\{h_{odd}[n]\}
$$

Do all real / causal h[n] have the property that:

$$
H(e^{j\omega}) = DTFT\{h_{even}[n]\} + j\ DTFT\{h_{odd}[n]\}
$$

where:

$$
h_{even}[n] = \frac{1}{2}(h[n] + h[-n])
$$

$$
h_{odd}[n] = \frac{1}{2}(h[n] - h[-n])
$$


texts: Help in deriving the update rule for equalizer

texts:
Considering an FIR system excited by an input $s_n$. If the input is derived from a mapping, and a +1/-1 BPSK signal is derived. Considering the Constant Modulus Algorithm (CMA). The channel is modeled as an FIR. If the received signal is fed to the equalizer, the output is : $y_n = w^T(n) x_n$ where $w^T$ is the equalizer tap vector. The update equation for weights $w_{n+1} = w  - \mu \Delta J(w)$

Will there be a modification in the derivation of the CMA and the update rule?


texts: Question about correlation of two signals

texts:
I'm working on some practice examples and in one of them I'm asked to find the convolution and the correlation of some signal pairs given. I know that correlation , given two signals x(n),g(n), is given by the following :
$$Φ_{xg}=x(n)\ast g(-n)$$ where asterisk is convolution. Now is this the same as $$g(n)\ast x(-n)$$

I think I read this somewhere $$Φ_{xy}(n)=Φ_{yx}(-n)$$ So the 2 convolutions give different results. Shouldn't the question specify which correlation I am to calculate? $$Φ_{xy}\\Φ_{yx}$$


texts: What Are The Semantics Of Wav-File Sample Values?

texts:
The title is the question.  To make this concrete, assume I have a standard 16-bit single channel wav file, with (although it doesn't matter) a 16000 Hz sampling rate.

Each sample magnitude will be a 16-bit value.

Is there a specified semantics or interpretation of these values, not relative to any physical value, but relative to the numerical systems of modern computers?  I.e., are those numbers to be interpreted as 16-bit twos-complement, or scaled into single precision floats between -1 and 1, etc?

At the highest level, I recognize this is a little pointless:  They're arbitrary scale.

But at a practical level, in the python world, librosa does one thing, scipy does another, and converting those different things into a spectrogram or analyzing with a vocoder can yield different results.

For skeptics, invite you to try the following:

import scipy
import librosa.display

np_eps        = np.finfo('f').eps
np_max        = np.finfo('f').max

x, sr = librosa.load(file, sr=16000)

plt.subplot(2, 2, 1)
D = librosa.amplitude_to_db(np.abs(librosa.stft(x)), ref=1.0)
librosa.display.specshow(D, y_axis='linear')
plt.colorbar(format='%+2.0f dB')

plt.subplot(2, 2, 2)
plt.plot(x)

sr, y = scipy.io.wavfile.read(file)
y = y.astype(float32)
plt.subplot(2, 2, 3)
D = librosa.amplitude_to_db(np.abs(librosa.stft(y)), ref=1.0)
librosa.display.specshow(D, y_axis='linear')
plt.colorbar(format='%+2.0f dB')

plt.subplot(2, 2, 4)
plt.plot(y)

Do not neglect the extra argument to librosa, forcing the sampling rate; in my case, the file I am using is from the Arctic corpus and is sampled at 16 KHz.  Do not neglect the colorbars, which will show you the very obvious shift of 90 dB from one to the other.  This should not be surprising, since 20 * log10 (32768) is approximately 90.3.

I assure you, gentle readers, this will have consequences further down the pipeline of speech vocoders, or neural networks (where normalization techniques might want you to reference to numpy epsilon to avoid negative numbers) or both.


texts: Is it possible to change sampling frequency of DSP kit if the sampling rate is fixed?

texts: I am using a DSP kit TMS320c6711. It has a fixed sampling rate of 8000 samples per second. As I have to do real time  filtering of frequencies above 5oHz with sharp cut off, I must reduce the sampling frequency to around 1000 samples per second. Suggestions are highly appreciable.


texts:
I wrote a wavetable oscillator that scans a file and extracts grains etc. The sound doesnt seem clear and I think it has aliasing.

What are ways to be able to play the audio file at any frequency without aliasing?

I know there should be interpolation, is there any specific type of interpolation that is better suited for recorded audio files?

Are there any other techniques to reduce this aliasing during playback?

I am unable to visualize how an audio file at any frequency, that is played back at different frequencies by jumping samples during indexing, aliases back into the waveform, and what depth of interpolation is needed.

Will scanning an audio file at higher frequencies cause aliasing due to index jumping, or is it rather precision loss? If the audio is already recorded, without aliasing, 44khz sample rate, bandwidth of 22khz and playback rate up to +6 octaves, is the aliasing added inside the real time signal?


texts:
I have a input voltage signal which has disturbances present in it. I am actually discretizing the input signal into say $1300$ samples and taking Hilbert transform of the signal by taking $N$ samples at a time and then taking next $N$ samples at a time leaving behind the 1st sample..
Input Signal

While the output keeps oscillating instead of being steady for $N=50$ samples, (via Scilab) my time step is $0.4\ \rm{ms}$ and input signal frequency is $20\ \rm{ms}$ and that's why i am taking $50$ samples,
Hilbert transform output for $N=50$ samples.

If I increase or decrease the samples, the response has much more ripples than obtained here.

Cant i get output something similar to this??

 this is nothing but true rms of the signal obtained in excel.

My actual objective is to obtain the instantaneous amplitude and phase of the signal using Hilbert transform and I want to know what must be the length of the window (at present i have taken $N=50$ samples) and how to move that? Blocks or point by point?

The entire signal is not available and it just getting recorded then, and within the full signal getting recorded i need to obtain instantaneous frequency and phase.


texts: How can detect the Left Ventricle (LV) Cavity into the MR dicom file?

texts:
This is MR Modality Link

MR Modality

In this file where i find the LV Cavity in above link file. Also how to draw Tumor Marker from over it


texts: Transfer function estimation using system identification

texts:
I am posting a control system problem here because of this post reference.
Question:
I have input and output datasets. These are graphically shown below.

The blue line is output and the pink line is input which was applied during specific time intervals. This dataset has been obtained from a machine that has an existing manually tuned the PID controller.
I want to make an LQI or more perfectly tuned PID controller using bode plots. However, to do that, I will have to identify the transfer function or corresponding state space model.
What is the detailed procedure to obtain the transfer function from the above data set (because this data has been obtained from a machine which already has a manually tuned PID controller controlling it)? The manually tuned values are $K_p=2.2, K_i=0.01, K_d=2$.

My attempt
I used a system identification application and obtained a state space model with a $75 \%$ fit. I tried implementing a PID controller for it and it required large $K_p, K_d, K_i$ values for the required design parameters.

The gain values used were

These cannot be implemented on a system since it has comparatively very low $K_p, K_i, K_d$ values (manually tuned).


texts: How to derive equation for second order interpolation of soft knee cutoff in a compressor?

texts:
I am reading Audio effects: theory, implementation and application by Joshua D. Reiss and Andrew McPherson.
They explain that, in the simplest case with just a hard cutoff threshold, you have a piecewise linear equation
$$y_G = \begin{cases}
x_G & x_G <T\\
\frac{x_G -T}{R} + T & x_G \geq T
\end{cases}
$$
with a corresponding diagram like

However, in many cases, you may not want this hard cutoff. You may want it to be softer like this

To do this, Reiss proposes the following formula

I found a PDF of a presentation he did where he called this "second order interpolation", but I have no idea what that means other than I think it's second order because of the quadratic terms in this version.
How does Reiss derive the middle formula in this piecewise equation?
Reference:

https://www.eecs.qmul.ac.uk/~josh/documents/2012/GiannoulisMassbergReiss-dynamicrangecompression-JAES2012.pdf
http://c4dm.eecs.qmul.ac.uk/audioengineering/compressors/documents/Reiss-Tutorialondynamicrangecompression.pdf
Differences between soft knee and hard knee in Dynamic Range Compression (DRC)?


texts: QAM constellation coder - transmitting a file

texts:
The file with size 6 MB is transmitted in 100 seconds using a constellation coder QAM with a constellation of 6 bit/symbol

$ { M = 10^6, B = 8} $

Tasks:

1. What is the signal bandwidth if the band exceedance factor is 0.5?

$W_m = (1+α)R_s$

$W_m = (1+0.5)*80 000 = 120k$

2. What is the signal sampling frequency at the modulator output if the oversampling ratio equals 3.

$ T_S = 6*0.0000020833 = 0.0000125 $

$ n = T_s / T_p $

$ 3 = 0.0000125 / T_p  = 0.00000416666 = 4.17 µs$

$ F_p = 1/T_s => F_p = 239808.15 Hz = 240 kHz  $

3. How numerous is the set of symbols of the constellation used?

$ M = 2^k $

$ k = 6 $

$ M = 2^6 = 64 $

4. What is the modulation performance?

$ R_b / W_m = 480k / 120k = 4 $

5. How long will the signal transmission last if the band exceedance factor is 1.

And this is the only task I have no idea how to do. Any suggestions?

I would also be very grateful I you check the other tasks, if they were done correctly, or if I made any mistakes. Thank you in advance!


texts:
I have created a simple plot of waveforms using matplotlib of 2 wave files on python. Here is the code:

import matplotlib.pyplot as plt
import numpy as np
from statistics import mean
import wave
import sys

plt.figure(1)

class Audio:
    def __init__(self, audio):
        self.audio = wave.open(audio,'r')
        self.signal = self.audio.readframes(-1)
        self.signal = np.fromstring(self.signal, 'Int16')
        self.fr = self.audio.getframerate()
        self.time = np.linspace(0, 100, num=(len(self.signal)))
        self.fft = np.fft.fft(self.signal)

    def plot(self):
        plt.title("Audio waveforms")
        plt.plot(self.time, self.signal, '.')

gana = Audio('gana.wav')
humm = Audio('humm.wav')

gana.plot()
humm.plot()
plt.show()

The output is this:

What I want to do is remove the noise from both the files. I have seen many examples that use FFT to remove noise. I have taken out the fft of the signals of these wave files but the problems I am having is I have no idea of what to do after taking out the FFT. Do I have to plot the FFT or do something else?

P.S The noise in the green wave form is hiss and in the blue is background sound. So is there any way that I could get rid of both of these noises? Thanks.


texts: Visualisation of bandwith of filter with sampling

texts:
I'm having some trouble visualising (and understanding) the filtering of a sampled signal.

If I understand correctly, after sampling a signal, it is desired to have it go through a lowpass to retrieve the original signal again. The literature states the following below.

$\omega_0 < B < \omega_s - \omega_0$

Where $\omega_s$ is the sampling frequency and B is, what I assume, the bandwith of the filter.

I do not understand what is meant by the above expression and the literature is kind of crap at explaining it.

I was hoping someone can visualise this for me.


texts:
I'm trying to understand why Matlab function upsample causes contraction in the spectrum, like the following figure shows.

Also showed in mathwork online document.

enter link description here

Filtering after upsampling removes the imaging artifacts, but contraction still stays, as shown below. Could anyone help explain in details why?

enter link description here


texts:
If $h(n)$ is the impulse response of the discrete LTI system, $x(n)$ is the white noise process with variance $\sigma_x^2$ and $y(n)$ is the output, prove that

  $$\mu_y = \mu_x \, \sum_{k=0}^\infty h(k)$$

  where $\mu$ represents the mean.

My approach:

$$\mu_y= \mathbb{E}[y(n)]
=\mathbb{E}[x(n) \star h(n)]$$
Where, '$ \star$' means convolution
Now how can I proceed?


texts: Nyquist sampling theorem / Nyquist - Shannon theorem evaluation over M-PSK

texts:
I am trying to simulate an M-PSK Tx-Rx System on Simulink and analyse what the effect of sampling rate reduction would have on it. More, particular I am trying to prove what Nyquist sampling theorem states that the sampling frequency should be twice the maximum frequency component.

The communication chain is as follows:

Random integer generator -> M-psk mapping (modulation) -> RRC (interpolation) -> up-conversion (carrier) -> sampling (sampling rate reduction) -> down-conversion ->Interpolation (if necessary)-> RRC (decimation) -> M-psk De-mapping ->BER

The figure bellow confirms the above arrangement.

So far I have proven that by reducing the sampling rate I simply achieve to “mix down” the signal (all the BW that sits around the carrier) within some bandwidth limits defined by the updated (reduced) sampling frequency. That stops at some point (rendering the signal unrecoverable) where the sampling rate becomes to small; probably confirming Nyquist-Shannon theorem that the sampling rate suffices to be 2 times the BW.

So, no matter how high the carrier frequency is, the signal will always return in between some frequency limits defined by the sampling frequency (-fs/2 to fs/2) following a cyclic manner (modulo). The new (aliased) frequency component will sit apart from the sampling frequency following : f_alias = abs(Fs-Fc).

For example, if Fc is exactly 2 times higher than Fs then the spectrum will return to zero (eg Fs= 8e6 Fc = 16e6 @ around 0), if is 2MHz more then it will sit 2Mhz further (eg Fs = 8e6 Fc =10e6 @ around +2Mhz), if is exactly 4 times higher then again around zero etc (eg Fs= 8e6 Fc = 32e6 @ around 0).

The above simple experiment contradicts the Nyquist theorem that states that the sampling frequency should be at least 2 times greater that the maximum frequency component.And my question is how can I evaluate that theory given a similar communication system arrangement?

If all that above is correct then the selection of and ADC's rate would merely be based on the BW of a signal to be sampled rather than the signal's modulating frequency (carrier)?


texts: Is there any research work on the possible data structures that can be used to represent an architectural design of a building?

texts: I am looking for a data structure, that can be used for description of a floor plan of a building. Its not about the format, I am asking about any other way to describe the floor plan of building, rather than in traditional map-based structure, preferably in some sort of tree or graph based structure, which can be derived from a text based description of the house.


texts:
I'm currently trying to implement some of the methods found in this paper on intelligent equalization - http://www.aes.org/e-lib/browse.cfm?elib=16792 -

The first part of the process is to build a "target equalization curve" via the means of "averaging the spectra of all songs in the dataset".

In MATLAB, if I do the following.

mag_1 = abs(fft(file));
mag_2 = abs(fft(file_2));

The two magnitude spectra have different resolution and number of frequency bins right?  Does it make sense to just mean(mag_1, mag_2)?  I can't seem to find anything that goes over averaging multiple magnitude spectra and it isn't defined in the paper how they achieved the averaged magnitude spectra for the ideal target curve.

Any help would be appreciated.


texts:
I'm working on a project where I need to generate a LUT for Hamming Window to perform a FFT.

I'm using this function (doing the table with MATLAB):

n = 0:4095;
LUT_length = length(n);
WIN_HAM_v = round(4095*(0.54-(0.46*cos(2*pi*n/(LUT_length-1)))));

I have some doubts if the function is correct


texts:
How does formant shifting in synthesizers such as Reason's Maelstrom work?

I've been pointed to the PSOLA algorithm among others, but I cannot understand how exactly does it result to formant shifting, because its base operation seems to be "formant preserving pitch shifting".


texts: How to check if a signal is mainly composed of a few impulse peaks?

texts:
I am working on a problem that produces different types signals, and I need to detect the peaks after low-pass filtering, then find_peaks from scipy. For example, the one below:

I am however faced with some signals that have an impulse shape, as in:

or,

Filtering such signals destroys the peaks' locations, and sometimes lead to merging two or more peaks into one. Thus, I am not able to extract the peaks correctly.
Is there a method to infer if the signal has impulses (a few to many spikes)? This would help me directly extract the peaks easily.
I highly appreciate any input.


texts: how to use image filters for spectrograms?

texts:
Using MATLAB, I can import an audio signal and use some FFT function to get the spectrogram into the workspace. My question from here, is how to then process this spectrogram with image filtering techniques that usually import some image. For example, here is the code that filters an image:

J = imread('example.jpg');
K = somefilter(J);

my question is the imread part; as you can see the image is imported to J and then this J is used for the processing function below. However, can I simply input my spectrogram data from the workspace directly into the filter argument? or is there some other step that needs to happen before?

Also, most image filtering matlab codes mention sometime the input is grayscale image or RGB, color, etc. My question is: how is a spectrogram in Matlab classified as? does it matter for the filter? do I need to convert the spectrogram then?


texts:
I'm doubtful about which of these is correct:

$\frac{E_b}{N_0}=\frac{SNR}{k}\frac{T_{symbol}}{2T_{sampling}}=\frac{SNR}{k}\frac{L}{2}$
$\frac{E_b}{N_0}=\frac{SNR}{k}\frac{T_{symbol}}{T_{sampling}}=\frac{SNR}{k}L$

$L$ is the interpolation factor and $k=log_2(M)$ (number of bits per symbol).

There's just a difference of a $1/2$ factor. I know where this relation comes from, and that $1/2$ depends on noise power being $\sigma_N^2=N_0/2$ (baseband) or $\sigma_N^2=N_0$ (passband). Since I want to generate $\sigma_N^2$ given a certain $E_b/N_o$, I'm not sure which one to pick.

If I have $y(t)=x(t)+n(t)$, with $n(t)$ being AWGN and $x(t)=s(t)\cos(\omega_c t)$, so $x(t)$ is passband and $s(t)$ is baseband, I'm not sure if I should pick the second one because I'm adding it to $x(t)$ or the first one.


texts: How to convert from Laplace Domain to Time Domain?

texts:
I want to convert the following equation from Laplace domain to continuous time domain:

$F(s) = \frac{-2 m k v R}{2 m R s^{2} + m k s + 2 k R}$

m, k, v, R are all constants.

If I can factor or put this into partial fractions, I can use the table here to convert: https://lpsa.swarthmore.edu/LaplaceZTable/LaplaceZFuncTable.html

eg.

$ \frac{1}{s+a} = e^{-at}$

But I don't know how to break it down into that form. Is it doable? If so, how?

I have learned how to convert Laplace into the z-domain but I have found some problems with that. In particular, I need continuous time equations to set up the [n-1] and [n-2] etc. samples for the initial run or I won't get useful outputs.


texts: Decimator effect on wide sense stationary input

texts: I've seen that the output of a decimator when a WSS process is passed through remains WSS. I am not able to immediately see why this is. What is a good explanation of why the signal maintains stationarity?


texts: Which Fourier Transform properties should be applied first when they are combined together?

texts:
A sequence $x[n]$ has Fourier Transform $X(e^{jw})$, and the Time shfiting, Frequency shifting and Time reversal properties are shown below, respectively :
$$x[n-n_d] \Rightarrow e^{-jwn_d}X(e^{jw}) $$
$$e^{jw_On}x[n] \Rightarrow X(e^{j(w-w_O)}) $$
$$x[-n] \Rightarrow X(e^{-jw}) $$

Suppose we were given that the Fourier transform of a signal $Y(e^{jw})$ is $3e^{j4w}X(e^{-{j(w-w_O)}})$, and we want to work backwards to find the original sequence.  If I use Time reversal first then:
$$x[-n] \Rightarrow X(e^{-jw})$$
followed by Frequency shifting:
$$e^{jw_On}x[-n] \Rightarrow X(e^{-j(w-w_0)}) $$
and then finally Time shfiting and Linearity:
$$3e^{jw_O(n+4)}x[-n+4] \Rightarrow 3e^{j4w}X(e^{-j(w-w_0)}) $$

So, I have obtained $y[n] = 3e^{jw_O(n+4)}x[-n+4]$,but if I had changed the order of properties which I applied,I would have gotten a slightly different answer. Should certain properties always be applied first when dealing with such problems?


texts:
How to convert grayscale intensity input to something like the following image?

It seems that no response (0, black) is converted to dark blue, then it goes through other colors to red.

How are these colors derived? Is there any standard gray-to-RGB conversion function for these?

I have tried representing intensity as hue in CIE-LCH model and then converting to RGB through CIE-Lab. But the result weren't as nice as these MRI and PET scans.

Another approach may be a rainbow, i.e. going from high frequency (violet-blue) to lower frequencies (dark red). But is there a conversion formula for this?


texts: Feature extraction from spectrum to aid categorisation

texts:
i have two sets of signals whose spectrum looks like shown in the figure. The first column of signals belongs to category 1 and the second belongs to category 2.
Looking at the spectrums I want to identify features which will differentiate the two.
Eg:mean,variance etc.


texts: frequency components of discrete Gaussian-windowed signal

texts:
When we are do a discrete Gabor transform on a 1D signal, we apply a Gaussian window on a segment of the signal. However, due to the tapering ends of the Gaussian window, the magnitude of the data points will be reduced more towards the ends, This implies that the the frequency components of the segment of the signal will be suppressed when we apply FFT to it. After performing some experiments, I found that at least 2 periods of the signals is required, so that no misjudgement of the behaviour of the signal will happen.

What is the range of frequencies in which the magnitude of the FFT is still representative of the frequencies of windowed signal?


texts: How to design all-pass filter by specifying phase requirements?

texts:
I want to delay the part of signal above (or below) specified frequency, but I don't find any books or tutorials to tell how to design is. The commonly used biquad APFs (such as the one introduced in this cookbook) produce a zigzag-shaped phase response, that is not I want.

How can I do it? Thanks for a lot!


texts: Extracting the phase of a jittering sine wave using python

texts:
The background:
I conducted an experiment where I rotated an object (with three-fold symmetry) many times and recorded the projection of the object. I should be able to tomographically reconstruct the object using an inverse Radon transform.

However, first, I need to know the exact angle of the object at each step. Unfortunately, our experiment was not perfectly stable and the object did not rotate in perfectly consistent steps. The object has a "three-leaf clover" shape, and if I plot the "asymmetry" of the projection (top-bottom/(top+bottom)), then it should look like a sine wave that has a peak every 120 degrees of phase (due to the three-fold symmetry).

Anyway, here is the real problem: In the figure below, I plotted the real asymmetry of the images in blue. In red, I plot a simple sine wave of what the asymmetry should be if there was no experimental jitter. If you look at the blue data with your eye, you can clearly see where it deviates from the model and you can infer what the actual phase of the wave should be. However, I have no idea how to convince python to assign a phase to each datapoint based on its actual asymmetry value. Any ideas??

Here is the data and the code used to make the simple plot:

text = """0.58018836  0.54719772  0.35009603 -0.07093418 -0.38495609 -0.56433715
 -0.58238798 -0.58197763 -0.79733289 -0.91822821 -0.88177671 -0.71078292
 -0.1013863   0.53490235  0.50354981  0.05747148 -0.49258035 -0.86043242
 -0.73366264 -0.11736485  0.40830017  0.75403686  0.76892275  0.74578394
  0.67759482  0.58632213  0.38615145 -0.04784894 -0.58773909 -0.80827248
 -0.70176645 -0.3018048   0.33288852  0.75568507  0.52563443 -0.09271778
 -0.74360484 -0.93751009 -0.84840137 -0.62956237 -0.18022792  0.23716647
  0.58545184  0.70021713  0.61602417  0.39086982 -0.09830602 -0.73635755
 -0.83658372 -0.36582787  0.12293982  0.51728333  0.72395964  0.77727396
  0.71986029  0.43646683  0.0821097  -0.33721973 -0.75678586 -0.59755965
  0.08786262  0.59804408  0.76638545  0.4829408  -0.01587402 -0.42493757
 -0.63032865 -0.70086101 -0.71030878 -0.55375893 -0.16074062  0.42243708
  0.78208721  0.58585221  0.14783457 -0.32854651 -0.73844163 -0.74110153
 -0.28917785 -0.03469134 -0.03392741  0.10028943  0.5273808   0.8099916
  0.60887728 -0.0270428  -0.64313701 -0.8553194  -0.82824477 -0.60014313
 -0.2471473   0.2198018   0.56404379  0.71217775  0.81460569  0.66100642 """

import numpy as np
import matplotlib.pyplot as plt

asym = np.array(text.split())
plt.plot(asym)

step = 11.76
theta = np.linspace(0,np.shape(asym)[0]*step,np.shape(asym)[0])

plt.plot(np.sin(3*theta/180*np.pi),color='r')
plt.show()


texts: Retrieve frequency value of an 2D frequency spectrum

texts:
How do you retrieve the values of the frequency component occurring in a 2D freq spectrum?

I have this, and if I wanted to know the values of the white dots.  what should i do ?

I apolgize beforehand for asking this question.


texts:
I'm looking at heart rate variability (HRV).  Frequently, in the Fourier transform  I see two peaks that are close together.   How could I determine whether this is two peaks, or just one peak and noise? Please see figure below. The units of the x axis is hertz. Thanks.
ps. I take a forty minute HRV measurement. For the figure below, I used a window of 1000 points, about 17 minutes, in the middle of the measurement. I'm not doing anything else to the data.
Below is also 100 data points from the source data.


texts:
I'm currently working on a digital audio OFDM project and found that even if I explicitly specified the same sample rate (11025 Hz), the transmitter and receiver still differ slightly in it(lesser than 1 Hz). My professor asked me to use a 'multirate interpolating filter' to synchronize the two sample rates. Could you give me some hints about it?

By the way, I'm programming with C++ and the portaudio library. Thanks.


texts: Windowing a signal by using Hanning / Hamming or... zero padding?

texts:
My question is about the use of Hanning or Hamming windows versus "square" windows.

In my understanding, I will use a Hanning or Hamming window to limit the spectral leakage. But, on the second hand, it means I will loose the begining and the end of the signal as they are faded to zero.

I read a solution which is to "ponderate" each temporal temporal signals get by the inverted DFT as I would do with "tiles" covering each other on a roof. And, that's right, if I "recover" with a 66% ratio, my "ponderated" signal looks like the original signal multiplicated by a 1.5 ratio (excluded the 2 first and 2 last frames of my whole signal on which I can't ponderate the fading).

If I'm using that solution, it means I need to wait the 2 next samples before releasing my sample "recalculated" by the inverse DFT :

first third of the sample n ponderated by samples n-1 and n-2 (I already got them so, I don't count them)
second third of the sample n ponderated by samples n-1 and n+1
last third of the sample n ponderated by samples n+1 and n+2

So, if my signal is sampled at 44.1 kHz and contains 2048 frames (to get a "good" spectral resolution), its length is approximatively 50 ms. So, my delay before releasing that signal will be 1 * 50 (=sample n) + 2 * 50 (samples n+1 and n+2) = 150 ms. It is too huge if I'm working in live conditions like in a concert.

So, my question is : why not using a square window and zero-padding BEFORE and AFTER my sample ? If I zero pad my signal before and after my original sample, won't I avoid the spectral leakage ?

Thank you for all your answers.


texts: How to remove the boundary effects arising due to zero padding in scipy/numpy fft?

texts:
I have made a python code to smoothen a given signal using the Weierstrass transform, which is basically the convolution of a normalised gaussian with a signal.

The code is as follows:

#Importing relevant libraries
from __future__ import division
from scipy.signal import fftconvolve
import numpy as np

def smooth_func(sig, x, t= 0.002):
    '''
    x is an 1-D array,
    sig is the input signal and a function of x.
    '''
    N = len(x)
    x1 = x[-1]
    x0 = x[0]

# defining a new array y which is symmetric around zero, to make the gaussian symmetric.
    y = np.linspace(-(x1-x0)/2, (x1-x0)/2, N)
    #gaussian centered around zero.
    gaus = np.exp(-y**(2)/t)

#using fftconvolve to speed up the convolution; gaus.sum() is the normalization constant.
    return fftconvolve(sig, gaus/gaus.sum(), mode='same')

If I run this code for say a step function, it smoothens the corner, but at the boundary it interprets another corner and smoothens that too, as a result giving unnecessary behaviour at the boundary. I explain this with a figure shown in the link below.
Boundary effects

This problem does not arise if we directly integrate to find convolution. Hence the problem is not in Weierstrass transform, and hence the problem is in the fftconvolve function of scipy.

To understand why this problem arises we first need to understand the working of fftconvolve in scipy.
The fftconvolve function basically uses the convolution theorem to speed up the computation.
In short it says:
convolution(int1,int2)=ifft(fft(int1)*fft(int2))
If we directly apply this theorem we dont get the desired result. To get the desired result we need to take the fft on a array double the size of max(int1,int2). But this leads to the undesired boundary effects. This is because in the fft code, if size(int) is greater than the size(over which to take fft) it zero pads the input and then takes the fft. This zero padding is exactly what is responsible for the undesired boundary effects.

Can you suggest a way to remove this boundary effects?

I have tried to remove it by a simple trick. After smoothening the function I am compairing the value of the smoothened signal with the original signal near the boundaries and if they dont match I replace the value of the smoothened func with the input signal at that point.
It is as follows:

i = 0
eps=1e-3
while abs(smooth[i]-sig[i])> eps: #compairing the signals on the left boundary
    smooth[i] = sig[i]
    i = i + 1
j = -1

while abs(smooth[j]-sig[j])> eps: # compairing on the right boundary.
    smooth[j] = sig[j]
    j = j - 1

There is a problem with this method, because of using an epsilon there are small jumps in the smoothened function, as shown below:
jumps in the smooth func

Can there be any changes made in the above method to solve this boundary problem?

Also I tried removing the zero padding in the fft source code and replaced it with a constant value, but it gave undesired results.
Can you suggest a way of removing this zero padding in the scipy fft source code?


texts: Octave: How to attenuate all frequencies except one?

texts:
I am implementing a morse code decoder in Octave. It already works, but not on very noisy signals, like the following:

I had the following idea to reduce the noise:

Perform fft to find the "beep" frequency of the morse tone
Use a bandpass filter to attenuate all other frequencies.

Unfortunately, I apparently can't do step 2. After filtering, my signal looks like this:

This shows about the first two milliseconds of the signal.

My code looks like the following. Please tell me:

whether my basic idea is promising
how to improve my code in a way that it actually does what I want

[x, f_sampling]  = wavread(filename);
t = fft(x);
l = length(t);
magnitudes = (abs(t/l))(1:l/2+1);
f = f_sampling*(0:(l/2))/l;
[peaks, locations] = findpeaks(magnitudes);
[maximum,index] = max(peaks);
f_main = f(locations(index))

f_cutoff = [0.1 0.9]*2*f_main/(f_sampling)
[b, a] = cheby1(20, 1, f_cutoff);
 y = x / max(x);

y = filter(b, a, x);


texts: Available Online Recordings for Acoustic Echo Identifiablity

texts:
I am trying to add a few samples for a publication on identification of clean signal from the one which has echo on it. For this purpose I need data set that includes playback of an already available signal in a big room/hall and the recording of the signal with a mic in this hall that includes the echo.

Does anybody know of any available online source for such data samples.


texts: What is logarithmic mapping in a filter bank such as the one in the human ear?

texts:
"The spacing of the filters in a filter bank used to make an ordinary spectrogram is linear, whereas an auditory filter bank has a spacing that corresponds to the way in which sinusoidal frequency maps onto place on the basilar membrane, a so-called tonotopic map. For frequencies of about 1kHz and above, this mapping is logarithmic."

I read this in an article about how the ear works in processing sound. Can someone explain mathematically how the mapping for frequencies of about 1 kHz and above is logarithmic? Also it would be appreciated if anyone can suggest some literature regarding this topic.


texts: What are the differences between fft and dct?

texts:
what are the differences  between Fourier transform and cosine transform.....
i try to find the typical answer
the difference ,from where ,the size of input and output form,in time complexity,the domain ,the effect on image itself .


texts: Processing signals with phantom fundamentals

texts: Is there a formula to calculate the frequency of a phantom fundamental based on the harmonics, and would this formula be applied directly to the energy bins of an FFT?


texts:
I am doing comparison between two wave files, How to calculate phase difference between those files?
Second wave is generated by making some changes in the first one .I have tried calculating phase like this:-

y, sr = librosa.load('first.wav')
D = librosa.stft(y)
magnitude, phase1 = librosa.magphase(D)

y1, sr1 = librosa.load('second.wav')
D1 = librosa.stft(y1)
magnitude1, phase2 = librosa.magphase(D1)

then calculated the difference as (phase1-phase2)

Is this the correct way??


texts:
I have to calculate SNR of an image. I have camera through which continuous images are getting generated, my requirement is to calculate the SNR in db for each image.
I have seen lot of the examples where they have added a noise to image then computed the SNR value.
I have to calculate SNR of each image independently without adding any noise.
please help.


texts: Why can the number of zeros be no more than the number of delay elements in a signal flow network?

texts:
Let $N$  be the signal graph representation network of the system function (in rationale form) of a discrete-time LTI system.  For a network $N$  with
no loops, the impulse response is no longer than the total number of delay elements
in the network. From this, we conclude that if a network has no loops, then the system
function has only zeros (except for poles at $z = 0$), and the number of zeros can be no
more than the number of delay elements in the network.
Q: Why  can the number of zeros  be no
more than the number of delay elements in the network?


texts: Frequency identification as a convex problem

texts:
Suppose I have two time series, $x(t)$ and $y(t)$: $x$ is the output of a sine wave generator, and $y$ is the generator's frequency setting. Given a training set of $(x, y)$ pairs, Is it possible to formulate a convex/quasi-convex optimization problem to build a predictor of some new $y$ given $x$?

The system is certainly non-linear, since you can double the signal amplitude without changing the output (the frequency). It seems like you need to perform an argmax in the frequency domain, which I believe is a quasi-convex operation. This limits what functions you can compose, so something like computing residuals may no longer be convex. My hunch though is there may be a clever recasting of the problem, or a duality argument, etc.

I'm sure neural nets are capable of it, but for sake of discussion assume I don't have the data/processing power.


texts: Can an unwrapped phase lag diagram have values beyond $-\pi$ to $\pi$?

texts:
I have a PLL that tracks the phase lag of a resonance system. When unwrapped, my system's phase lag has values beyond $-\pi$ to $\pi$.
After unwrapping phase lag, are there any bounds on the phase difference between signals? I assumed that the bounds were always from $-\pi$ to $\pi$ since a phase difference greater than $\pi$ could be expressed as a negative difference between $-\pi$ and $0$.
x-axes should say "Frequency"


texts:
I am trying to implement segmentation on an image. after a difference of Gaussians Process and applying weights to my edges I want to go through edge linking by local processing. I have seen this example with the van on multiple tutorials. I tried the second "real time" simplification.

what I do not understand is:

do I implement this algorithm on the gradient of the original image or the weighted edgels (I have tried both and got squat)?
If not. How can I use the gradient to implement an algorithm that will help me link the edgels?

This is) the image I am working with.

This is the weighted DOG output.

This is what I got from creating the g function from the original image and after filling gaps vertically

and horizontally I got this.

I am following the same thresholds as in the example.

This is a code snippet of my "gap filling":

        k=25
        I2 = ~imopen(~padarray(g,[0,1],1),strel('line',k,0));
        HrizontalConnection  = I2(:,2:end-1);
        I2 = ~imopen(~padarray(g.',[0,1],1),strel('line',k,0));
        VerticalConnection  = I2(:,2:end-1);
        I=HrizontalConnection|VerticalConnection.';


texts:
I am trying to downsample from $44100Hz$ to $8000Hz$, and I am not sure what I am doing wrong here.  $44100 / 8000$ is $5.5125$, so I am getting hung up a little on the fractional remainder which is causing my implementation to play back slightly slower than it is supposed to.  It seems that I need to somehow use the fractional remainder to compensate for this, but I am not sure how...

def downsample
  s = low_pass_filtered_signal_by_half_target_sample_rate
  samples_to_discard = (sample_rate.to_f / target_sample_rate).floor //   (44100.0 / 8000.0).floor

  s1_n = 0
  for (n = 0; n < N; n++) {
     if (n % samples_to_discard == 0) {
       s1[s1_n] = averaged(n, samples_to_discard)
       s1_n += 1
    }
  }
end

def averaged(start_index, samples_to_discard)
  sum = 0
    for (n = start_index, n < start_index + samples_to_discard; n++) {
      sum += s[n] || 0
    }
  sum / samples_to_discard
end


texts:
Evaluate the Z transform of $x[n] = n^3$ where the signal is two sided.

I have tried using the basic definition of the Z transform ie.,

$$X(z) \triangleq \sum_\limits{n=-\infty}^{+\infty} x[n] \, z^{-n} $$

But what I get obviously is an infinite unbounded sum of sequence. Does that mean Z transform does not exist? IS there any alternate methods to obtain the transform? Please help.


texts: How do I deal with pixel intensity overflow when performing discrete convolution on an image and the result is greater than 255?

texts:
This is my first post on this forum, and I have been working to come up with a solution for a final project in my hardware design class using system verilog on an Altera Cyclone FPGA.  I've been having trouble finding a solution for overflow when convolving a 3x3 pixel intensity matrix (RGB, values from 0-255) with a specific kernel (blurring, canny edge detection, embossing, etc.)

The question I have is this: What should I do to deal with a result that is greater than 255?  I've seen some references suggesting to "wrap the values around" from 255 back to 0 and up, and also have seen some references suggesting to multiply all pixels by a "shrinking factor" (e.g. .8) to avoid ever having to deal with pixel intensity overflow.

Sadly, I've found few references explaining how to deal with this issue in hardware realization.  Am I missing something here?  Any assistance would be greatly appreciated.

EDIT:

Specifically, what if you convolved a matrix of

[255, 255, 255; 255 255 255; 255, 255, 255] * [1, 1, 1; 1, 1, 1; 1, 1, 1]

?  Wouldn't your result contain values outside of your range of 0-255?


texts: Remove noise on edge (compression deffects)

texts:
I have the images of cartoons, that were compressed. The example:

They have such noise, that is not easy to delete.

Even though the pixels are on the gray background, the noise pixels can be of very different colors and if I look at the histogram (of grayscale image) of such image I see a lot of bins around the bin of main color. But I can't just delete the colors of value, that are not of main color, because there is a risk to delete some important features (as eye, for example). Also I've tried to posterize image (example presents 8 colors), but some pixels still remain.

Also I've tried median filter, but it is not able to vanish such heavy noise ( I use 3*3 filter).

Please, give some methods that can effectively remove the noise in this case. I would be grateful for any help!


texts:
I'm still researching this but though I'd post asking about experiences. Looks like st.com offers a variety of STM32 Discovery dev boards for ARM CPUs to, among other things, play with DSP. But searching this forum didn't find much.

Ideally I'd like something with built-in stereo ADC/DAC but I don't seem to see that. I see daughter boards that have stereo mics and stereo outputs. But mics are a non-starter for me. I need line level inputs.

Or this video recommends a Wolfson board which you splice to the STM32.

https://www.youtube.com/watch?v=TO66UN64-os

What's the difference of the M7 core over the M4? Where to order?  How do you power the board?

Thanks.

Dennis


texts:
I realize this is quite basic stuff but I am having trouble with the following question:

  Determine whether the following system with input $x(t)$ and output $y(t)$ is BIBO stable: $$y(t) = \frac{1}{x(t)}$$

My intuition says that it is not stable, just from thinking about the case where $x(t) = 0$. However, this would result in $y(t)$ being undefined, not going to infinity. Any help would be appreciated.


texts: Polyphase filter bank channelizer FFT or IFFT

texts:
I am referring to the Xilinx application note Polyphase Filter Bank Channelizer. The receiver portion should have FFT but the pdf mentions IFFT in both Tx and Rx.
The paper by F. J. Harris, C. Dick and M. Rice, "Digital receivers and transmitters using polyphase filter banks for wireless communications," in IEEE Transactions on Microwave Theory and Techniques, vol. 51, no. 4, pp. 1395-1412, April 2003 also mentions the same.
How did Xilinx get proper results by implementing IFFT in both Tx and Rx?


texts:
I´m stuck in a deduction analysis of the variance of a gaussian white noise signal in a "integrate-and-dump detector" of a baseband data transmission receiver, where $n(t)$ is white noise with double-sided power spectral density $N_0/2$ [W/Hz]

I can understand all the steps except when they deduce

How do you get to this last expression?
Thank you.


texts: Extracting original signal from overlapping area

texts:
Considering this image, where I have a spectral decomposition of a light bulb.

Due to the non point shaped light source, the spectral decomposition using a diffraction grating in front of my camera leads to overlapping spectral components. As a human though I can still perceive the structure of the main four or five spectral components because I know they should have the same shape as the lamp. Knowing the structure of the lamp, is it possible to isolate each of those spectral components with their correct amplitude, and what algorithm would be used in that case?


texts: Extracting Peak Frequencies Using FFT vs. Time Domain Peak Finding

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

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

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

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

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

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

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

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

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

Thanks!


texts: 3D wiggle plot for an analytic signal: Heyser corkscrew/spiral

texts:
Just reading The Analytic Impulse, A. Duncan, 1988, I met the name "Heyser corkscrew" for the first time in my DSP life, for a 3D display of a cisoid or complex exponential $e^{i\omega }$ (often denoted  $\text{cis}\, x$ for "cosine plus $i$ sine").

The view is mentioned as Euler formula's 3D view in the context of circular polarization:

This representation is quite common for analytic signals, but this name does not seem standard. Questions are:

What is the history of this Heyser name, with a precise reference?
What are the (earlier) history and names for this representation  (going back to the Fourier oblique view)?

After @MBaz hints, the term "Heyser spiral" seems a little more common. It is named after Richard C. Heyser, and can be interpreted as a 3D variant of Nyquist plots... Indeed, in  A Matter Of Frequency: The Nyquist Plot Explained, one can read:

  Dennis Gabor proposed the “analytic signal” which Richard Heyser
  developed into what has since been named the Heyser Spiral. The
  Nyquist plot is the “end view” shadow of the complex analytic signal


texts: Upsampling vs downsampling. Which to use when?

texts:
Downsampling reduces dimensionality of the features while losing some information. It saves computation.
Upsampling brings back the resolution to the resolution of previous layer.
My question is which is better or when to use which one?


texts: Complex output after inverse FFT of a real signal

texts:
I have a real one dimensional signal s (light absorbance in a flow cell), which has significant noise and some periodic noise after performing a deconvolution of $S$ from $S_o$. Basically fft($S$) was divided by the fft($S_o$). Let us call the output as $R$ vector. In order to remove the periodic noise, I replaced the noise region in $R$ with 0+0i in a given range. The ifft($R$) is complex. When I plot t and ifft($R$), MATLAB warns that the imaginary parts of complex R arguments ignored  . The output looks as desired, free from noise.

Forum members mention that this can arise due to numerical precision, however the imaginary part in ifft($R$) is on the order of 0.0006i.

Is there is a better window (simple one) than this rectangular window in the frequency domain?

Even in this simpler example which one can test is a simpler code, generates a complex inverse output.

t = [0:1/80:60]'; % Time
x = sin(2 * pi * 2 * t) + sin(2 * pi * 0.05 * t); % Signal
X = fft(x);
X([1:6, 4780:4801]) = 0; % Removing low frequency 0.05 Hz
z = ifft(X); % **z turns out to be complex as well**
figure(1)
plot(t, z)
hold on
figure(2)
plot (abs(X))

Thanks.


texts: Implementing a Band-Pass Filter from Scratch

texts:
I am trying to implement a band-pass filter from scratch. Here is an article which I thought explained the nuts and bolts of how to build one; the author combines a low-pass with a high-pass filter (convolving both) to create a band-pass filter. I applied it to a mixture of sine waves, but the result was not satisfactory.

import numpy as np
from numpy import sin, linspace, pi
from pylab import figure, plot, show, title, xlabel, ylabel, subplot, savefig, xlim, ylim
from scipy import fft, arange

def plotSpectrum(y,Fs):
    n = len(y) # length of the signal
    k = arange(n)
    T = n/Fs
    frq = k/T # two sides frequency range
    frq = frq[range(n/2)] # one side frequency range

    Y = fft(y)/n # fft computing and normalization
    Y = Y[range(n/2)]

    plot(frq,abs(Y),'r') # plotting the spectrum
    xlabel('Freq (Hz)')
    ylabel('|Y(freq)|')

Fs = 150.0;  # sampling rate
Ts = 1.0/Fs; # sampling interval
t = arange(0,30,Ts) # time vector

ff = 0.2;   # frequency of the signal
y = sin(2*pi*ff*t)
ff = 0.8;   # frequency of the signal
y = y + sin(2*pi*ff*t)

figure()
subplot(2,1,1)
plot(t,y)
xlabel('Time')
ylabel('Amplitude')
subplot(2,1,2)
plotSpectrum(y,Fs)
savefig('out1.png')

fL = 0.1
fH = 0.3
N = 33
n = np.arange(N)

hlpf = np.sinc(2 * fH * (n - (N - 1) / 2.))
hlpf *= np.blackman(N)
hlpf = hlpf / np.sum(hlpf)

# Compute a high-pass filter with cutoff frequency fL.
hhpf = np.sinc(2 * fL * (n - (N - 1) / 2.))
hhpf *= np.blackman(N)
hhpf = hhpf / np.sum(hhpf)
hhpf = -hhpf
hhpf[(N - 1) / 2] += 1

hp = np.convolve(hlpf, hhpf)
filtered = np.convolve(y, hp)

t2 = np.linspace(0,30.,len(filtered)) # time vector
figure()
subplot(2,1,1)
plot(t2, filtered)
xlabel('Time')
ylabel('Amplitude')
subplot(2,1,2)
plotSpectrum(filtered,Fs)
savefig('out2.png')

For the original signal:

I get an output like this:

To confirm I could filter this signal, I used a butterworth filter,

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

filtered = butter_bandpass_filter(y, fL, fH, Fs, order=2)

That worked fine (although I had to play around with the order a little bit, higher orders of this filter does not work either). What is wrong with the first filter?

Regards,


texts: Frequency domain vs time domain equalization

texts:
It is said that a frequency domain equalizer is equivalent to an infinite length time domain equalizer, why is that?

And if this is the case, why hasn't frequency domain equalization been used earlier?

Has time domain equalization any advantages when compared with frequency domain equalization?


texts: How to find phase shift and process I/Q channel

texts:
I am currently working of 24GHz FMCW radar signal processing. The Module I'm using has 2 receiving antennas and each has 2 outputs (I and Q channel). So overall I've 4 channels. I'm using combination of harware and software for the processing in frequency domain. Range and velocity is fine, if I sample any one channel out of those 4 (either I or Q) I can find the R and V values.

Now I'm interested in finding the direction of a object, when it is moving (away from sensor or towards sensor) and the position of the object (right of the sensor or left or the senor, at the best the horizontal angle).

In the datasheet of the radar sensor I read 2 lines in the features, that gives me some idea in finding the direction and location. Below are the lines:

"two receiving antenna for phase comperision operation"
"I/Q channels for direction of motion discrimination"

Now I'm thinking how to find out the position or the horizontal angle of the target using phase comperision of the 2 antennas. Below is the diagram I made for analysis.

Surely the 2 antennas receives reflection from a single target with a phase shift. Now the question is how to I find the phase shift processing the I channel (let's say) from the 2 antennas. I'm not aware of any DSP technique, so want suggestions.

Also the direction of motion, what technique can be ised with the I and Q channel of a antenna to deduce it. I'm again on poor dsp knowledge, how do I process the I and Q channel.

I appreciate any suggestion of improve my knowledge and know about more dsp techniques in signal analysis.


texts: What's the Difference between Convolution Kernel and Point Spread Function in the Context of Image Convolution?

texts:
When you use the deconvolution method to make the blurry image sharper, you will have to estimate the Point Spread Function. Is there a difference between this PSF and an image kernel?

Second question: Is the Point Spread Function in the time domain or frequency domain? Or can it be in both?